A dispersive homogenization model for composites and its RVE existence

An asymptotic homogenization model considering wave dispersion in composites is investigated. In this approach, the effect of the microstructure through heterogeneity-induced wave dispersion is characterised by an acceleration gradient term scaled by a “dispersion tensor”. This dispersion tensor is computed within a statistically equivalent representative volume element (RVE). One-dimensional and two-dimensional elastic wave propagation problems are studied. It is found that the dispersive multiscale model shows a considerable improvement over the non-dispersive model in capturing the dynamic response of heterogeneous materials. To test the existence of an RVE for a realistic microstructure for unidirectional ﬁber-reinforced composites, a statistics study is performed to calculate the homogenized properties with increasing microstructure size. It is found that the convergence of the dispersion tensor is sensitive to the spatial distribution pattern. A calibration study on a composite microstructure with realistic spatial distribution shows that convergence is found although only with a relatively large micromodel.


Introduction
The heterogeneity of the microstructure of composite materials causes dispersion in wave propagation associated with dynamic loading. This dispersion phenomenon, also referred to as micro-inertia, is a result of local motion of the microstructure driven by multiple wave reflections and transmissions occurring at the interfaces between the constituents. Correct evaluation or tuning of the dispersion properties of composites can be important for engineering applications, for instance impact-resistant components or devices where composites are subject to high-rate loading [1][2][3][4][5] or metamaterials with microstructures designed to show particular effective behavior [6][7][8]. Computational modelling of composites subjected to stress-wave loading typically involves three length scales, i.e. the size of the macroscopic structure, the characteristic stress wave length and the length scale of the microstructure (in fiber-reinforced composites related to the fiber diameter) [8,9]. The macroscopic length scale can be much larger than the microstructural length scale. In the case where the stress wave length is also much larger than the typical microstructural size, there is no significant transient effect within the microstructure and micro-inertia is negligible [8]. Therefore, the overall dynamic response can be solely described by averaging density and moduli. However, for a stress wave which is only few times larger than the microstructural length scale, the dispersion becomes evident and averaging properties are not sufficient to describe the dynamic response [8]. In order to account for the dispersion phenomenon, multiple models have been introduced. One type of model is the gradient elasticity model with high-order spatial derivatives of strains, stresses and/or accelerations as reviewed by Askes and Aifantis [10]. However, identification of the length-scale parameters of this method, especially in a multidimensional context, is not totally clear. Besides, Willis' elastodynamics homogenization model derives an effective constitutive relation which introduces non-classical coupling between effective stress and effective velocity and coupling between effective momentum and effective strain [11][12][13][14]. Recent development of generalized homogenization models enrich the macroscale displacement with additional generalized degrees of freedom of Bloch modes following the lines of Willis' model, see Sridhar et al. [15]. These methods are suitable for linear elastic (layered) periodic materials while an extension to materials with random microstructure is not always straightforward.
An alternative type of method which provides a more flexible framework (e.g. consideration of complex random microstructure or material nonlinearity) is representative volume element (RVE) based multiscale approaches. Initiated by Hill [16], the RVE can be defined as a characteristic sample of heterogeneous material that should be large enough to contain sufficient composite micro-heterogeneities in order to be representative, however it should be much smaller than the macroscopic structure size. Several definitions of an RVE are introduced in literature, as reviewed by Gitman et al. [17]. Generally, a micromodel can be considered representative if further increase of its size of RVE does not lead to changes in the homogenized properties. Typically, statistics studies using numerical computations are used to evaluate the homogenized physical properties (e.g. elastic properties, thermo-conductivity) or effective response (e.g. the effective stress) of a range of micromodels with increasing size, for instance in [17][18][19][20][21][22][23][24]. It should be noted that the size of an RVE depends on the specific investigated morphological (e.g. volume fraction) or physical properties (e.g. thermal or elastic) [18,20].
Multiscale methods assume multiple (at least two) spatial and (or) temporal scales [25,26]. In multiscale models finer-scale problems are considered in a (statistically equivalent) representative volume element (RVE) and information of the finer-scale is hierarchically passed into a coarser scale by bridging laws. Based on a multiscale virtual power principle [25] as a notion of a generalized Hill-Mandel lemma, Pham et al. [8] and Roca et al. [27] developed computational homogenization schemes in which the transient dynamics equations are resolved at macroscopic and microscopic scale. Asymptotic homogenization with higherorder (or first-order) expansions was proposed in [28][29][30] to capture wave dispersion and attenuation within viscoelastic composites. Fish et al. [31] introduced a general purpose asymptotic homogenization scheme in which the microinertia is resolved by introducing a eigenstrain term and is valid for nonlinear heterogeneous material. This method was further investigated by Karamnejad and Sluys [32] for impact-induced crack propagation within a heterogeneous medium. The aforementioned methods have shown certain capabilities in capturing wave dispersion for strictly periodic heterogeneous structures with simple microstructures where the RVE can be clearly defined as a unit cell. However, in a realistic composite material, the microstructure has a random nature. Therefore, the question of RVE existence needs to be answered before multiscale methods can be employed.
In this paper, a dispersive multiscale model based on asymptotic homogenization is reviewed and the existence of an RVE for this method for unidirectional fiber-reinforced composites is investigated. This paper is organized as follows: in Sect. 2, the dispersive model based on asymptotic homogenization technique is introduced. The accuracy of this numerical method for 1D and 2D elastic wave propagation is demonstrated in Sect. 3. In Sect. 4, a statistics study is performed to investigate if an RVE exists for this homogenization approach for realistic fiber reinforce composite microstructures. In Sect. 5, a batch of calibrated numerical samples based on experimentally determined spatial distribution pattern is tested for the convergence of homogenized properties.

Dispersive homogenization model
In this section, the asymptotic homogenization model originated from Fish et al. [31] and later explored by Karamnejad and Sluys [32] is described. This method allows for a decoupling of the equation of motion into a two-scale formulation. The effect of microscopic dispersion is quantified by a socalled "dispersion tensor", which is related to the acceleration influence function of the microstructure. The acceleration influence functions for microstructures with one and multiple inclusions are demonstrated as examples.

Two-scale formulation
Considering a composite structure that is in dynamic equilibrium with prescribed displacements and stress boundary conditions and given the initial conditions for displacement and velocity (see Fig. 1a), the linear momentum equation reads with boundary conditions at the boundary surface ∂ ζ = ∂ uζ ∪ ∂ tζ and initial displacement and velocity conditions where σ ζ , ρ ζ , u ζ and n ζ are stress, density, displacement and surface outward normal, respectively. The superscript ζ denotes that quantities are defined within the composite domain. For simplicity, a linear elastic material law is considered herein, namely where ε ζ is strain and S ζ is the elastic stiffness tensor which is for two-phase (i.e. inclusion and matrix) composites a piece-wise constant function of spatial coordinates. Extension to nonlinear material behavior can be done by using the instantaneous stiffness tensor as elaborated in Fish et al. [31]. Perfect bonding between inclusion and matrix is considered here while decohesion can be possibly included through the eigendeformation concept [33] or a cohesive crack formulation [34]. A relevant study on dispersive multiscale formulation with consideration of material damage is introduced by Karamnejad and Sluys [32].

Asymptotic expansions
In the asymptotic expansion approach, two spatial coordinate systems are introduced, macro-scale coordinates x defined in the macroscopic homogeneous domain and micro-scale coordinates y in an RVE domain with heterogeneous microstructures, see Fig. 1b. The y coordinate is related to the x coordinate by y = x/ζ with 0 < ζ 1. Any physical field variable, for example, the displacement field u, is a function of spatial coordinate x, y and also the physical time t. These physical fields are assumed to be y-periodic in RVE domain , namely f (x, y) = f (x, y + kl R ) in which vector l R = [l 1 , l 2 ] T represents the basic period of the microstructure (in 2D), l 1 and l 2 are the lengths of RVE along the two directions and k is a 2 × 2 diagonal matrix with integer components. The choice for periodic boundary conditions for the microscopic field is motivated by superior convergence properties that have been demonstrated by Kanit et al. [18] and Fish [35] among others. Following Fish et al. [31], we can represent this function by an asymptotic expansion around a point x in powers of ζ , namely, in which the first term on the right hand side represents a macro-scale component while the second term represents a micro-scale component. By applying the two-scale spatial ∂ y i , the strain can be expressed as where ε m i j is the micro-scale strain and (·) (i, means the symmetric gradient of a function with respect to coordinate x. The stress is expanded as where σ m i j is the micro-scale stress and the micro-scale perturbation stress σ (1) i j satisfies σ (1) i j (x, y, t) = 0 ( 1 0 ) where (·) = 1 | | (·)d is the volumetric average of (·) within the domain .
The inertia is given as the following expansion where h st i (x, y, t) is the so-called acceleration influence function which satisfies ρ (1) st and h st i (x, y, t) is a y-periodic function satisfying the normalization condition

Weak form
The weak form of Eq. (1) reads in which the test function ω Integration of the two-scale functions over the composite domain and its boundary is carried out as The test function is expanded as where the relation ω (1) By applying the two-scale integration scheme Eq. (14) to Eq. (13) and using Eqs. (9), (10), (11), (12) and (15) and neglecting higher-order terms, it can be found that where ω i = 0 on ∂ u } and the dispersion tensor D M is introduced as where (·) M denotes the macro-scale quantity. Therefore, the weak form for the two scales can be derived as The solution of these weak forms is found with separate finite element schemes for the two scales. Details about the extraction of the macro-scale stiffness S M i jkl are explained in "Appendix A" and the solution of the micro-scale problem is illustrated in the next section.

Micro-scale problem
It can be observed from Eq. (19) that the micro-scale problem is treated as a "quasi-dynamic" formulation following Fish et al. [31]. Considering that S ζ i jkl (x) = S m i jkl ( y) and leaving out the first order remnants in Eqs. (8) and (9), the constitutive relation Eq. (6) can be rewritten as The micro-scale strain ε m kl (x, y, t) is divided into two parts: one is caused by the macro-scale strain u (0) (i,x j ) (x, t); the other one is driven by the macro-scale acceleration gradient. The micro-scale strain, from Eq. (8), is expressed as macro-scale strain related terms macro-scale acceleration gradient related terms (21) where H st k is a displacement influence function satisfying y-periodicity and H st k = 0, and f kl (x, t) is a macroscopic quantity proportional to the macro-scale acceleration gradient so that f (0) kl ≡ 0 in the absence of macro-scale acceleration gradient. It is to be noted that in Eq. (21) the first of the macro-scale strain related terms corresponds to a uniform strain on the microstructure while the second term is a fluctuation contribution to keep the microstructure to be in self-balanced state. Similarly, the first term of macro-scale acceleration gradient related parts in Eq. (21) represents a uniform acceleration gradient over the microstructure while the second terms refers to a fluctuation contribution to keep the microstructure in equilibrium under the excitation of a macroscopic acceleration gradient.
Considering the arbitrariness of u (23) and Of these, the second one is solved for the acceleration influence function h st k with finite elements over the micro-scale domain. Details about the solution procedures of Eq. (24) are given in "Appendix B". By use of Eq. (17), h st k gives the dispersion tensor which is used in the macro-scale weak form, i.e. Eq. (18).

Acceleration influence functions
The effect of inertia in the micro-scale is considered by an eigenstrain [36] and therefore only the solution of a typical quasi-static balance problem is needed. As it is pointed out in Fish et al. [31], it is possible to get a closed form solution of the micro-scale balance Eq. (19) in the case of a 1D composite bar with two different phases of material with elastic material constants of E 1 , ρ 1 and E 2 , For a microstructure with one circular inclusion, the six acceleration influence functions are plotted in figure 2. The elastic modulus, Poisson's ratio, mass density of the inclusion and matrix are E i = 200 GPa, ν i = 0.2, ρ i = 10,000 kg/m 3 and E m = 2 GPa, ν m = 0.2, ρ m = 4000 kg/m 3 . The volume fraction of the inclusion is 60% and its diameter is 5 µm. From Eq. (21), it is known that the acceleration influence functions h mn k can be interpreted as the first-order microstructural correction to the eigendeformation field as triggered by the macro-scale acceleration gradient. For instance, the negative gradient of h 11 1 along y 1 -direction in the domain occupied by the stiff inclusion implies that the true acceleration induced strain in the inclusion is smaller than that in the matrix domain. For a microstructure with 25 × 25 randomly positioned inclusions while keeping the other properties unchanged, the influence function h 11 k (k = 1, 2) is plotted in Fig. 3. Some observations can be made: (1) the gradient h 11 (1,y 1 ) in the inclusions phase is still negative and has similar magnitude among all the inclusions, and (2) both h 11 1 and h 11 2 show regions where peaks are higher and regions where peaks are lower as a consequence of the mesostructure.

Properties of dispersion tensor
The magnitude of the dispersion tensor depends on several characteristics of the microstructure. (i) The contrast of material properties between the constituents. In terms of a one-dimensional two-phase composite bar, if the wave impedance is the same for the two phases, the dispersion tensor is null [31]. This is consistent with the fact that when the wave impedance is the same, wave dispersion due to material heterogeneity is not present. To investigate the influence of material contrast on the dispersion tensor, the dispersion tensor is computed for a single microstructure with 15 × 15 fibers and a fiber volume fraction of 60%. The Young's modulus and density of the inclusion are varied proportionally as E i = cE m and ρ i = cρ m for c ∈ [2,4,9,16,30], respectively, with the Poisson's ratio unchanged. The six independent components of the dispersion tensor D M i jkl , normalized with respect to the values of the sample with the lowest property contrast, D M i jkl (c = 2), for the five cases are shown in Fig. 4a. It can be seen that the magnitude of all six components are increasing with larger properties contrast, which means that the expected dispersion effect should be larger for a higher contrast composite. (ii) The diameter of inclusions. Similarly, the dispersion tensor is computed for five cases with the same material properties and microstructure with 15 × 15 inclusions but with different inclusion diameter, d i ∈ [0.005, 0.025, 0.1, 0.25, 0.5] mm. The variation of D M 1111 as function of d i is plotted in Fig. 4b along with a quadratic fit. It can be observed that the relation between D M 1111 and d i can be described very well as a quadratic relation and the same was found for the other components of the dispersion tensor. This quadratic relation can also be found in the closed-form solution of the dispersion tensor for a one-dimensional two-phase composite bar given in [31]. This means that for larger inclusions, the influence of dispersion is larger. To ensure a meshindependent solution of D M 1111 , five different mesh sizes of the case with the fiber diameter d i = 0.25 mm are solved, with the number of elements ranging from 22,624 to 288,800. It is shown in Fig. 4c that after the mesh size has reached the medium size, the dispersion tensor component D M 1111 has reached converged values. Therefore, this geometry y 1 ( y 2 ( mesh size is adopted for the study described in Sects. 4 and 5.

Comparison with direct numerical simulation
In this section, elastic wave propagation in a periodic composite medium is simulated to investigate the performance of the introduced dispersive homogenization model in comparison with direct numerical simulation (DNS). Two cases are considered, one with a one-dimensional microstructure and one with a two-dimensional microstructure.

1D elastic wave propagation
Elastic wave propagation in a periodic two-phase composite bar is studied similar to Karamnejad and Sluys [32]. Geometry and properties are such that wave propagation is purely one-dimensional. The material properties of the two phases are E 1 = 200 GPa, ρ 1 = 10,000 kg/m 3 , ν 1 = 0 and E 2 = 2 GPa, ρ 1 = 4000 kg/m 3 , ν 2 = 0. The wave impedance con- trast between the two phases is E 1 ρ 1 /E 2 ρ 2 = 250. The left end of the bar is fixed. Two types of loading are considered. Firstly, a incoming harmonic wave is imposed on the right end of the bar so that the horizontal displacement u(t) on the right edge satisfies in which A 0 = 0.025 mm represents the magnitude, f = 50,000 Hz is the wave frequency and H (·) is the Heaviside step function. The bar consists of 100 repeating unit cells for a total length of 500 mm as shown in Fig. 5. Four numerical models are considered: two single-scale models, namely the fine heterogeneous model (DNS) and the fine model with homogenized properties, and two multiscale models, namely the coarse non-dispersive model and the coarse dispersive model. The two fine models consist of 4500 quadrilateral elements. For the one with homogenized properties, the elastic modulus and density of the two phases are prescribed to be the homogenized values, i.e. E M = 4.926 GPa and ρ M = 7600 kg/m 3 , respectively. This model therefore neglects wave reflec-tion and transmission at material interfaces. By contrast, the heterogeneous DNS model considers the heterogeneity of the composite bar and is therefore considered as the reference exact solution of this problem. The two coarse multiscale models use 500 quadrilateral elements with one Gauss-integration point corresponding to a unit cell domain. Homogenized elastic properties are used and the analysis is performed with and without dispersion tensor. The main dispersion tensor component D M 1111 for this simple microstructure can be obtained by a closed-form solution following Fish et al. [31]. The dispersion tensor is evaluated numerically as It is shown in Karamnejad and Sluys [32] and Fish et al. [31] that the ratio between macroscopic wave length l M and the unit cell size l m determines the extent of dispersion. If l M is more than five times larger than l m , the dispersive effect can be neglected. The dispersive model shows considerably better accuracy than the non-dispersive model within the first pass frequency band 1 . In this case, the wave length is calculated by l M = E M (ρ M ) −1 / f = 16.1 mm, while the unit cell size l m = 5 mm. Therefore, the ratio l m /l M = 0.3106 in this case, which corresponds to a shorter wave length than Karamanejad and Sluys [32]. The displacement field along the bar is plotted for four typical time instants in Fig. 6. It can be seen in the DNS reference solution that the input sinusoidal pulse is not maintained during the propagation and there exists significant amplitude decay. The resultant displacement field shows sharp kinks representing high velocity gradients. This flutter characteristic is mainly due to the high wave impedance contrast between the two phases of the structure, which causes wave reflection. The result of the fine homogenized model shows a well maintained profile although with a small amount of leading oscillation related to the discretization. A comparison between the DNS model and fine homogenized model shows that dispersion is indeed 1 According to Andrianov et al. [37], a periodic elastic composite behaves like a discrete wave filter. A discrete pass frequency band and stop frequency band structure is formed. Whenever the wave frequency is within the stop frequency band, its magnitude is exponentially attenuated such that the wave is not able to propagate. Only when the wave frequency is within the pass frequency bands, propagation is admissible.
significant for the studied wave length. The solution from the coarse non-dispersive model shows very strong oscillations due to numerical dispersion effects caused by a coarse mesh and the match with the DNS model solution is poor. The coarse dispersive model shows a relatively smooth response. This is because for the dispersive model there is no physical interface in the macro-scale model, so the wave can not "feel" the interface. Nevertheless, there exists reasonably good agreement between the dispersive solution and the DNS model solution before wave reflection at the left edge (a-c) and after reflection (d) although the dispersive model predicts a stronger decay of magnitude for the oscillation at the rear of the wave. A higher order homogenization scheme could result in a higher accuracy but with more computational costs, see for instance [38]. Secondly, loading which mimics an impact-induced loading pulse is considered, to demonstrate the capability of the introduced dispersive numerical model for impact problems. The problem setting is the same as described in Fig. 5 except that the horizontal displacement on the right edge is prescribed as capturing dispersion with a discretization at macroscale that is coarser than the microstructural resolution.

Two-dimensional wave propagation
Next, elastic wave propagation in a material with a two dimensional microstructure subjected to an incoming sinusoidal wave is considered. The geometry consists of 100 repeating microstructures with a total length of 57.04 mm.
There are two phases of materials, circular inclusions with a diameter of 0.5 mm and a surrounding matrix. The top  Fig. 10). For the homogenized one, the material properties of the complete domain are prescribed to be the homogenized values of each phase. By contrast, the DNS model keeps the different properties of each phase and therefore the dispersion caused by material heterogeneity is naturally included. The coarse models are discretized with 500 linear quadrilateral elements at the macro-scale while each Gauss integration point corresponding to a micro-scale problem solved within a unit cell domain as shown in Fig. 11. The coarse non-dispersive model neglects the contribution of dispersion at the macro-scale while the dispersive model adds the dispersion tensor contribution to capture the dispersion effect. The RVE for this structure is clearly identified as the unit cell. The dispersion tensor is found using the procedures illustrated in "Appendix B" as: The averaged horizontal displacement u(x), i.e. the volumetric average of horizontal displacement within a constant distance l s , along the bar for two time instants is plotted in Fig. 12. In this study, l s is equal to 1/500 of the total length of the bar. According to the reference DNS solution, the input sinusoidal wave does not maintain its profile during the propagation and it breaks into several pulses with obvious magnitude decay. The profile also shows significant oscillations caused by reflections at material interfaces. The magnitude decay of the DNS model exhibits a more gradual process than the 1D wave propagation problem considered in Sect. 3.1. This is due to the two-dimensional nature of the microstructure. The coarse dispersive model captures the magnitude of the wave well although it predicts a more smooth wave profile. The computational time per time step for the dispersive multiscale model is around 0.04 seconds while for the DNS model around 2.5 seconds is needed on the same system (a PC with a 16 GB of memory and 3.5 GHz Intel Xeon CPU).

RVE existence study
In the previous section, the accuracy of the introduced dispersive model has been validated for periodic composite structures in which an RVE can be clearly defined as the building unit cell. However, for engineering materials such as fiber reinforced composites, the assumption of a periodic microstructure is not representative. Spatial variations in the fiber distribution can affect the material response. The existence of an RVE for this multiscale approach needs to be assured before the dispersive multiscale model can be applied to simulate realistic composites with random microstructure. The definition of an RVE prescribes that the size should be sufficiently large such that further increase of the size does not lead to change in the homogenized properties (or effective response). Therefore, the homogenized quantities, i.e. the dispersion tensor and the stiffness tensor, of random microstructures with different sizes are calculated in this study. It is noted that with reference to undirectional fiber-reinforced composites circular inclusions are considered here while non-circular inclusions, for instance, ellipses with different aspect ratio, can be considered for other materials, following Li et al. [40]. In this section, the numerical tool used to generate random microstructures is firstly described, followed by a study of the spatial distribution features of generated numerical samples. Finally, convergence of the dispersion tensor and stiffness tensor for realistic composite microstructure with respect to micromodel dimensions are presented.

Generation of random microstructure with DEM
Following van der Meer [41], a Discrete Element Method (DEM) solver called HADES, developed by Stroeven [42], is used to generate numerical samples of random microstructures. DEM allows for generation of high packing density of granular samples, as the final configuration is a result of stochastic initial conditions and particle-to-particle collisions under external environment force and boundary conditions. The process for generating a two-dimensional numerical sample with desired packing density V i is described as follows (see Fig. 13): A large box-shaped reference body is created initially, which contains a predefined number of N i circular particles with either a given constant diameter D i or a given diameter distribution curve f (D i ). The initial posi-tions of the particles are random perturbations of horizontally and vertically aligned locations. The initial velocity of any particle is prescribed to be the same in the horizontal and vertical direction but with a randomly assigned magnitude and sign for each particle. The reference body is gradually diminishing with the same velocity in the horizontal and vertical direction (see Fig. 13a). Periodic boundary conditions are applied to the reference body to avoid any wall effect [17]. This means that during the simulation whenever a moving particle is crossing the edge (or corner) of the diminishing reference body, it reappears in the corresponding edge (or corner) during the DEM simulation (see Fig. 13b). The collision between particles are treated with the Hertz contact law [42]. The inter-particle contact can introduce energy dissipation characterized by contact damping. In the case that damping coefficients are non-zero, continuous energy loss due to contact results in a clustering effect on the resulting particle system. It will be demonstrated in Sect. 4.2 that this phenomenon can be clearly demonstrated with the spatial distribution pattern of the particle system. A minimum contact distance parameter d min is enforced so that any two particles can not be closer than this value, following the same hardcore model concept as in [43][44][45][46]. It can be expected that decreasing the size of the reference body has the effect of increasing the volume fraction of particles. The simulation stops when the reference body has decreased such that the desired packing density, i.e. volume fraction, is reached (see Fig. 13c). After this, a mesh is generated with GMSH [47] for the inclusions and the matrix.

Spatial point distribution analysis
For the random microstructure of composites reinforced with long fibers, the inclusions can be considered as discrete circular objects dispersed into a continuum matrix in a two-dimensional domain. The spatial positioning of these inclusions is stochastic for realistic composites due to the manufacturing processes and conditions. By treating the center of each inclusion as one point, statistical spatial point pattern analysis can be applied to characterize the microstructure. As it is commented in Bailey and Gatrell [48], the basic interest in analysing the spatial point process is in whether the observed events exhibit any systematic pattern, as opposed to being distributed randomly in space. Possibilities for these patterns are classified as first order or second order effects. First order effects relate to variation in the mean value of the process in space representing a global trend in the distribution of inclusions. Second order effects result from the spatial dependence in the process, representing a local effect. This can be described by the probability density function of the nearest neighbour distances and Ripley's K function [48].
In this study, numerical samples of composite microstructures generated by HADES are first evaluated for their (a) (b) Fig. 14 A representative sample of a batch A and b batch B spatial distribution pattern. By introducing the contact damping within the DEM solver, the inter-particle (inclusion) distances become smaller, representing local effects. To demonstrate the influence of contact damping, two batches of numerical samples are considered, one with nonzero contact damping (batch A) and the other with zero contact damping (batch B). The probability density function of the nearest neighbour distances and Ripley's K function are evaluated for these two batches of samples.

Nearest neighbour distances
The nearest neighbour distance (NND) measures the shortest inter-inclusion interaction. The N th (N = 1, 2, . . .) nearest neighbour distance is the distance between a randomly chosen inclusion and its N th closest neighbour in the studied domain [48]. It is assumed that the studied numerical sample is in some sense representative of any region from a realistic composite microstructures. Therefore, it is needed to determine at which size the NND distribution function has converged. Eight different sizes are considered, by varying the number of inclusions per sample N i with values among The diameter of the inclusions D i is kept as a fixed value of 5 µm and their volume fraction V i is 60%. The minimum contact distance d min is set to be 0.18 µm. In Fig. 15a, the probability density of the 1st NND at different distances corresponding to different sample sizes for batch A is shown. The value of the probability density at a distance h represents the "chance" of the 1st NND for an inclusion to be h. A detailed view of the probability density of the 1st NND at two typical distances for different sample sizes is plotted in Fig. 15b. It is noted that in order to eliminate the influence of the sample boundary on the probability density function of the 1st NND for a given numerical sample (e.g. Fig. 13c), a toroidal correction is applied. It assumes that the top and left of the sample domain is connected to the bottom and right, respectively, as if the sample domain is a torus [49,50]. It can be found from Fig. 15 that as the size of sample increases, the mean value of the 1st NND density function gradually converges and the standard deviation decreases. When the sample has 35 × 35 inclusions, the 1st NND function represents a size-independent function for this batch of samples. The same holds for the numerical samples of batch B. With this converged size, the probability density function of the 1st and 2nd NND for batch A and batch B are evaluated, with the mean value and standard deviation shown in Fig. 16. It can be found that the largest probability of 1st NND occurs around 5.21 µm for both batch A and batch B, however, batch A has a larger peak value than batch B. The 2nd NND shows the largest probability around 5.24 µm for batch A while the largest probability for batch B occurs at around 5.38 µm. Again, the peak value of batch A is larger than batch B. These observations show that the spatial distribution of inclusions of batch A is more clustered than that of batch B, which is a result of the contact damping included for batch A.

Ripley's K function
The nearest neighbour distance considers the point pattern on the smallest scale. Information on larger scales is considered by Ripley's K function [48]. The K function is the ratio between the expected number of inclusions within a circle of radius h of an arbitrary inclusion and the mean number of inclusions per unit area. A direct estimate of K (h) from a numerical sample is given by [48] (28) in which A is the area of the sample, N i is the total number of inclusions, d i j is the distance between inclusion i and j, I h (x) is an indicator function which has a value of 1 if the condition x is true otherwise 0, and w i j is a weighting function for edge correction following Zangenberg and Brøndsted [50]. A homogeneous point process with no spatial dependence is the Poisson process with K = π h 2 [48]. The mean value and standard deviation of the K function for batch A and B are plotted in Fig. 17 along with the Poisson solution. The standard deviations for the two batches are both very small. The

Convergence of dispersion tensor
The influence of micromodel size on the dispersion tensor for the two batches is investigated to examine the convergence of the dispersion tensor. With the methodology mentioned in Sect. 2, the dispersion tensor can be solved using finite element models. The material properties for inclusions are prescribed as the Young's modulus E i = 74 GPa, Poisson's ratio μ i = 0.2, density ρ i = 2500 kg/m 3 and the matrix properties are Young's modulus E m = 3.76 GPa, Poisson's ratio μ m = 0.3, density ρ m = 1200 kg/m 3 . There are six independent components for the dispersion tensor due to symmetry. 2212 are calculated to be very close to zero, as expected for a transversely isotropic material. However, certain different findings can be made for the two batches. It is seen that the mean value of the six independent components of the dispersion tensor for batch A is gradually approaching a constant value after N i becomes larger than 20 2 with decreas-ing standard deviations. For batch B, the mean value of the dispersion tensor components shows the trend of converging to a constant value after N i becomes larger than 30 2 with the standard deviations tend to decrease as well. The converged values of the dispersion tensor for these two batches seem to be different, which can be possibly due to different spatial distribution patterns. The dispersion tensor converges with a very large size for both these two batches, although for batch A the convergence rate is slightly faster.

Convergence of stiffness tensor
The influence of the numerical micromodel size on the stiffness tensor for the two batches is investigated to also examine convergence of the homogenized stiffness tensor. The stiffness matrix is obtained by using the numerical scheme explained in "Appendix A". Plane-strain conditions are assumed herein. The calculated mean and standard deviation of the stiffness tensor is shown in Fig. 19. It can be found that for both these two batches, the standard deviation values for the stiffness tensor are small compared with the mean values. The mean values already converge for relatively small micromodel size. The calculated mean values of the stiffness tensor for batch A and batch B in matrix notation are: To verify the results, the isotropy of S M is checked. It is known that for isotropic material under plane strain condition, the stiffness matrix is where E is the Young's modulus and ν is the Poisson's ratio. Therefore, an isotropic law should satisfy that 2S M 1212 /(S M 1111 − S M 1122 ) is equal to one. An error measurement S is introduced as According to Eqs. (29) and (31), the errors for batch A and batch B are calculated to be 0.0486% and 0.0496%, respectively. This shows that the calculated stiffness is very close to theoretical values. As it is seen in Eq. (29), the difference between any component of the calculated stiffness tensor for these two batches is less than 2%, which means that the stiffness tensor is not sensitive to the considered differences in spatial distribution of the microstructure.

Experimental calibration
It is demonstrated in Sect. 4 that the convergence performance of the dispersion tensor depends on the spatial distribution pattern of the inclusions. Therefore, realistic microstructures should be studied to evaluate the appropriateness of this method to be applied in real composite structures.
With image-analysis techniques, spatial distribution patterns can be extracted from snapshots of real composites (see e.g. Czabaj et al. [51]). Computational approaches can thereafter be used to generate random numerical samples according to the experimentally determined spatial patterns. In this study, numerical samples are generated with the DEM solver HADES for the microstructure of a CFRP composite, HTA/6376, with a fiber volume fraction of 59.2% [52,53]. First, the DEM settings are calibrated to match experimentally observed NNDs. Then, the RVE existence study from the previous section is repeated with these settings. For the calibration the number of fibers is fixed at 1296 (= 36 2 ) in line with the experimental observations on 1300 fibers. The adopted RVE size is large enough to consider the spatial distribution in the numerical samples as representative based on the observation in Fig. 15. A fiber diameter distribution function is predefined according to the experimentally determined size distribution (see Fig. 20). The contact parameters, i.e. damping and minimum distance d min of the DEM simulation are tuned such that the spatial distribution pattern of generated numerical samples can match the experimentally-determined patterns. Herein, 100 numerical samples are generated with the DEM solver. The mean value and standard deviation of the actual fiber diam- Distance h (μm) Probability density (µm −1 ) Yang et al. [43] Ismail et al. [54] This study Experiment [52] Fig. 21 The probability density function of 1st NND eter distributions used in the DEM simulations are shown in Fig. 20 along with the experimental quantification to verify that the input distribution is recovered. The calculated mean and standard deviation of the probability density function of the 1st NND for the calibrated numerical samples are demonstrated in Fig. 21 along with the experimental result and two reference solutions from other numerically generated fiber distributions by Yang et al. [43] and Ismail et al. [54]. It is seen that the probability density function of the 1st NND for the generated numerical samples has a good match with the experimental measurements. A considerably better agreement is found for the current study than for the two reference solutions. The mean and standard deviation of the probability density function of the 2nd NND for the numerical samples and the reference solutions are shown in Fig. 22. Again, the result for the generated numerical samples matches well with the experiment measurements. A higher agreement is found for this study than Yang et al. [43] while the solution from Ismail et al. [54] also shows a good match. It is therefore validated that the generated numerical samples Distance h (μm) Probability density (μm −1 ) Yang et al. [43] Ismail et al. [54] This study Experiment [52] Fig. 22 The probability density function of 2nd NND  Fig. 23 Convergence study of calibrated sample are sufficiently representative for the real microstructure of the CFRP composite, following the criterion proposed by Liu and Ghoshal [55]. A convergence study on the dispersion tensor is performed next by increasing the size of samples with the same DEM solver settings as the calibrated numerical samples. Again, eight different sizes are considered, and for every size 100 random numerical samples are generated and solved for the dispersion tensor with the FEM model. The influence of the micromodel size on the six independent components of the dispersion tensor is plotted in Fig. 23. The standard deviation of all six components shows a decreasing trend as the number of fibers N f is larger than 400 and the mean values converge to representative values. This means that for this specific composite material the dispersion tensor components converge although convergence is still relatively slow when compared with the stiffness tensor components as shown in Fig. 19 or the elasto/plastic response of composites as reported in [41].

Conclusion
In this paper, a multiscale model is introduced to capture wave dispersion in composites. By using asymptotic homogenization, it is found that the dispersion effect can be characterized by a dispersion tensor, the magnitude of which is dependent on the material property contrast of inclusion/matrix and the size of the inclusion. The dispersive multiscale model is applied to simulate elastic wave propagation problems in a bar with one-dimensional and two-dimension microstructures. Comparison with a DNS model shows that the dispersive multiscale model has a significantly improved accuracy, compared with non-dispersive homogenized models. To test if an RVE can be defined for realistic composites with random microstructures, the dispersion tensor and stiffness tensor are computed for random numerical samples at different sizes. It is found that the convergence performance of the dispersion tensor is considerably slower than that of the stiffness tensor and that the convergence depends on the spatial distribution pattern. Finally, a batch of calibrated numerical samples of CFRP composites is tested for the convergence of the dispersion tensor. It is argued that careful definition of microstructure geometries is required to achieve representativeness. After solving the incremental form of system equation K δu = δ f , the stiffness matrix S M is obtained by in which subscript p denotes the degrees of freedom of the three controlling nodes and subscript f represents the other free nodes.

Appendix B
To solve Eq. (24) (24) is further simplified by using the FEM formulation as where (1,1) ω (2,2) 2ω ( in which the cofficient matrix c = N . Equation (36) can then be treated as a constraint for the linear system of equations given by Eq. (35). The y−periodicity condition of h mn