Homogenization of the heterogeneous beam dynamics: The influence of the random Young's modulus mixing law

Abstract This paper concerns the homogenization of the dynamic response of Euler Bernoulli's beam with random Young's modulus. Considering the eigenvalue problem, special attention is dedicated to the homogenization residuals (correctors) analysis, i.e. the difference between the random heterogeneous solution and the homogenized solution. Several correlation (mixing) laws of the Young's modulus are considered and a dimensionless characteristic scale length, based on the correlation length, is introduced. The effects of the mixing law on the residuals are analyzed using numerical approaches both for sampling the random Young’ modulus and for examining the beam eigenvalue problem. Two measurements are introduced to estimate the residuals between apparent and effective solution: the normalized difference of the Young's modulus and the normalized difference of the modes' shape. The effect of the mode's order is also highlighted with reference to forced vibrations.


Introduction
It is known that the assessment of Representative Volume Element (RVE) is a fundamental stage in the homogenization [1]; in fact, by analyzing that volume, it is possible to estimate, in average, the mechanical characteristics of the heterogeneous composite material. Following [2], the RVE is the smallest material volume, sufficiently larger than micro-structure size, for which the spatially overall moduli are accurate to represent the mean response.
Except for material with periodic texture where the Periodic Unit Cell PUC can be detected, the RVE is not known a priori and then the problem is to assess the RVEs dimensions. This aspect has been taken into account in several papers [3,4,5,6,7,8] also considering the effects on material properties evaluation [9,10]; moreover, when the spatial variation of the micro-structure physical quantities cannot be ignored, the analysis of non-local interactions between heterogeneities has to be considered [11,12].
If on one side the RVE should possess statically homogeneous and ergodic proprieties [13], on the other a Statistical Volume Element (SVE), that is a mesoscale sample with finite dimensions, have to be considered in applications; the SVE's proprieties are described by the adjective apparent [14] as opposed to RVE's effective ones. In this case the errors assessment and the convergence rate of SVE to RVE become the main aspects to be considered. An equivalent approach is the periodization of random media by statically equivalent periodic unit cell [15] that, besides the classic composites, can be also applied to different material, for example masonry material [16,17,18,19].
Wanting to tackle the problem from a mathematical point of view and limiting our attention, without losing generalities, to the modelling of the mechanical problem in terms of differential equations (strong form), the homogenization can be regarded as the study of the asymptotic behavior of these equations [20]. Regarding differential operators with random coefficients the main contributions are reported in [21,22,23].
If the homogeneous solution it is only obtained in limit condition, i.e. when a characteristic scale length approaches to zero, the homogenization residuals (or correctors) analysis, that is the analysis of the difference between heterogeneous solution (apparent) and homogeneous solution (effective) is the dual problem of the convergence rate of SVE to RVE and also assumes great importance for applications. Taking into account elliptic random operators, some results have been obtained in [24,25,26,27]. Form a mechanical point of view, the problem has been firstly addressed in [28], for the case of bi-phase beam with different quasi periodic texture, and then extended, in [29], to a two-dimensional problem.
The homogenization in dynamics has been extensively studied for heterogeneous material and composites [30,31] with particular attention to wave propagation in low and high frequency range [32,33,34,35,36,37] and FGM beams [38]. Vice versa the attention the residuals analysis and size effect in dynamics field have obtained a very limited attention [39,40].
In order to the analysis in deep the homogenization residuals, the approach used in [39] is extended taking into account the correlation law (mixing law) of Young's modulus; taking into account the homogenized eigenvalue problem of the Euler-Bernoulli's beam, we shall highlight, by introducing a "scale length", the interaction among residuals, mixing law and mode's order.

Mechanical problem
If it be desired to find the relation between the v ε (x, t) stochastic field of the transverse displacement, with 0 ≤ x ≤ L, and the E (ω, y) random Young's modulus in the heterogeneous Euler Bernoulli's beam, this may be done by writing the following stochastic partial differential equation In the previous equation ε is the length parameter that characterizes the microstructure, y = ε −1 x.
It should be noted that, in Eq. 1, the γ linear mass density and the J transversal section inertia moment assume deterministic values constant along x and f (x, t) is the external transverse force that is deterministic. It goes without saying that it is necessary to respect the boundary conditions and the initial conditions. It is of great interest, also from the perspective of the superposition method, to determine the natural frequencies and modes of the beam; to this aim the free vibrations equation with f (x, t) = 0 have to considered. Assuming the solution, separable in space and time, v ε (x, t) = u ε (x) g ε (t), from Eq. 1 we obtain the following differential equations which correspond to a eigenproblem of a beam with random Young's modulus. Given the boundary conditions, for each choice of E (ω, y) we obtain the eigenvalues (λ ε ) 2 k and the eigenfunctions u ε k (x); (λ ε ) k are the natural frequencies and u ε k (x) are the normal modes of the beam, with k = 1, 2, 3, ... mode's order.
This corresponds to assume (Ω, , P ) a standard probability space, i.e. a set Ω : ω ∈ Ω equipped with a σ-algebra of measurable subsets and a countably additive non-negative measure P normalized by P (Ω) = 1; T x , x ∈ R d is a family of invertible measurable maps T x : Ω → Ω; moreover T x is an ergodic dynamical system [41].
In the following we assume for E(ω) Young's modulus a Gaussian law with µ E = µ(E) mean , σ E standard deviation and probability density function To complete the definition of the E (ω, y) ergodic homogeneous random field, the following correlation (mixing) model is assumed In order to analyze the convergence rate to homogeneous solution and the effects of the mixing laws, the E j Young modulus of each micro-structure elements (j = 1, 2, 3, ....N e ) has been generated using a pth-order autoregressive model; N e is the number of micro-elements with equal length that constitute the beam, in the simulation the maximum value N e,max = 200 was adopted.
Le E j be expressed by where ϕ i , i = 1, ..., p are the parameters of the model, w is a Gaussian white noise and D is a constant. Given the correlation coefficient vector ρ i , i = 1, ..., p the parameters can be obtained by the following Yule-Walker equation [42,43] Moreover, the w Gaussian white noise has zero mean and variance σ 2 w given by: and D is given by: In the application four values have been considered for the α. In the first case α = α 1 = 2 was assumed.
Let τ be the correlation length, that is the length for which the correlation loss is adequate or in other words the for ∆x > τ the Young's moduli of the microstructures can be assume uncorrelated. This can be expressed by (1/τ ) α =ρ (in the first applicationρ =ρ 1 = 0.0400 ).
For each α i i = 1, 2, 3, 4 we have considered N s = 10, 000 samples of beams with N e micro-elements and Young's moduli E js,j with j s = 1, ..., N s and j = 1, ..., N e . The numerical correlations (averaging the samples ones) are reported in Fig. 2. [ Figure 2 about here.] From the sampling of E, the random variable C = 1/E has been obtained checking, for each j s beam sample, that 0 < c 1 < E js,j < c 2 < ∞ so that 0 < d 1 < C js,j = 1/E js,j < d 2 < ∞. It is worth noting that the previous check also avoids indeterminate structural conditions when a statistically determined beam is used, as a clamped beam in the following.
As will be clarified in the next section, the Young's modulus of the homogenized beam isĒ = [µ(C)] −1 ; so that an accurate estimation is

Homogenization of the eigenvalues problem
Using a formal analogy between periodic media and statistically homogeneous ergodic media, the main results on the asymptotic behavior of random differential operators were obtained by [21] and [22].
From the results reported in these papers, it is possible to shown that in Eq. 3, when ε → 0, u ε (x) converges, for a.e. ω, to u o (x) that is the deterministic solution of the homogenized problem whereĒ is the homogeneous Young's modulus that is the reciprocal of the mean of the reciprocal of Young's modulus random variable; this results justifies the estimation given in Eq. 10.
Without loss of generality, we consider a simply supported beam in the following. It is well known that angular frequencies and modes (square root of eigenvalues and the eigenfunctions) are: So that, as ε → 0 in Eq. 3, we have λ ε k → λ o k and u ε k (x) → u o k (x) with k = 1, 2, 3, ...

Residuals measurements and numerical results
In order to check the convergence, to estimate its rate and to analyze the effect of the mixing law, sets of N s samples of beam with random Young's modulus has been simulated considering the probabilistic law (with σ E /µ E = 0.15) and exponential laws (with α i values) as previously described.
Given α = α i i = 1, 2, 3, 4, for any sample j s = 1, ...; N s , the beam was considered with different length L = N e ∆l where ∆l is the length of the micro-structure elements (equal for each micro-element) with Young's modulus (constant) E α,js,j where j = 1 : N e , Fig 1 B).
Let N τ be the number of micro-structure elements to describe the τ correlation length; in applications N τ = 5 and it was equal for each mixing law. The approaching to zero of the dimensionless characteristic scale length ε is numerically described by the approaching to zero of the ratio τ /L = N τ /N e ; this is obtained by increasing N e .
It should be noted that in order to perform the random sampling and the structural analyses, specified numerical procedures have been developed utilizing the Matlab program. Regarding the structural analysis, a classical finite elements numerical method have been used and its reliability has checked by comparison with theoretical and numerical results reported in literature for beam composed by collinear elements [44,45,46].
For each sample beam, corresponding to (α, ε, j s ), the λ α,ε,js,k natural frequencies have been evaluated where k = 1, 2, 3, ... is the mode's order. So that the Young's apparent modulus of the equivalent beam is evaluated by For each simulation j s = 1, ..., N s the following normalized error can be introduced: whereĒ is the effective Young' modulus, see Eq. 13, that has been estimated using Eq. 10.
that measures the homogenization residuals versus the scale length (by ε), the mixing law (by α) and the mode's order (by k).
Let be µ(∆E α,ε,k ) the mean of the∆E α,ε,k random variable. The relationship µ(∆E α,ε,k ) versus k mode's order is shown in Fig. 3 for different values of ε and α.
[ Figure 3 about here.] It should be note that the converge rate∆E α,ε,k → 0 when ε → 0 is depending on k. Moreover the influence of the α exponent is evident: a lower α value corresponds to a greater correlation and larger residuals are expected.
A concise representation of this behavior is shown in Fig 4 with reference to k = 1 (A) and k = 5 mode's order (B).
[ Figure 4 about here.] The probability density functions (pdf ) of the residuals for k = 1, varying ε scale length and α mixing law exponent, are reported in Fig. 5. The Gaussian probability density functions, with first and second moments equal to the histograms ones, are reported in the same figures, it should be noted that the fitting is accurate. Similar results have been also found for the others k values.
[ Figure 5 about here.] Setting a ε value is equivalent to adopt a beam with L length, that is with N e collinear micro-elements characterized by random Young's modulus with a prefixed mixing law, given by α; so that, the obtained results, in terms of natural frequencies and modes, correspond to the statistical sample of this case. The homogenization result still holds with the following condition: the apparent coefficient is not deterministic but a measurable random variable.
An important results is obtained: the corrector, which measure the difference between the heterogeneous solution and homogeneous solution, converges in distribution to a Gaussian process also when the correlation function of the random coefficients is no longer integrable (that is an uniformly mixing condition but not strong).
If ε → 0 then the random Young's modulus becomes an ergodic field and the apparent modulus converges to the deterministic effective value: µ(∆E ε,k ) → 0 and the probability density function degenerates into a Dirac function.
An other measurement of the homogenization residuals can be obtained considering the mode's shape. For each sample beam, corresponding to (α, ε, j s ), the u α,ε,js,k ∈ N f em +1 discrete shape of the natural modes can be evaluated u α,ε,js,k = (u α,ε,js,k ) v v = 1, ..., N f em + 1 (19) where N f em is the number in finite element using in numerical analysis (to assure the solution accuracy N f em >> N e ).
Considering Eq. 15, for each simulation j s = 1, ..., N s the following normalized error can be introduced: where x v is the abscissa of the v th -node.
that measures the homogenization residuals versus scale length (by ε), the mixing law (by α) and the mode's order (by k). Let µ(∆u α,ε,k ) the mean of the previous random variable. The relationship µ(∆u α,ε,k ) versus ε and α is reported in Fig 6, A) for k = 1 and B) for k = 5.
[ Figure 6 about here.] The effects of the mixing law exponent on the homogenization residuals is highlighted by comparing the histograms of∆u α,ε,k ; for ε = 0.025, some comparisons, relative to k = 1 and k = 5, are shown in Fig. 7. It should be noted that an adequate accuracy is obtained when these histograms are fitted by the Gamma probability density function, with parameters obtained from statistical samples, Fig. 7.
[ Figure 7 about here.] Eventually, it should be noted that the obtained results highlight that modes order effects, on apparent Young modulus and modal shape, have to be considered in particular way for forced vibrations. In case of force with colored spectrum and by using the modal superposition method, greater is the modes order for which the resonant phenomenon is expected, larger must be the size of statistical volume element (SVE) to estimate the material mechanical characteristics and obtain an adequate response.

Conclusions
Considering the dynamics problem of the Euler-Bernoulli's beam with random Young's modulus the homogenization residuals are analyzed. Special attention has been dedicated to examine the effect of the mixing law of the Young's modulus; the exponential law has been adopted varying its exponent. Numerical procedures both for sampling the random Young' modulus and for examining the beam eigenvalue problem have been utilized. Adopting a correlation length, a dimensionless characteristic scale length has been introduced to estimate the convergence rate of the heterogeneous solution to the homogeneous solution. Two measurements have been introduced to analyze the residuals: the normalized difference between apparent and effective Young's modulus and the normalized difference between apparent and effective modes' shape. For both measurements, it has been highlighted that the mixing law plays an important role in the homogenization residuals behavior; a greater correlation corresponds to a lower convergence ratio that is also influenced by the mode's order. Moreover, with reference to first measurements, the obtained results show that the correctors converges in distribution to a Gaussian process even in case the correlation function is not longer integrable (that is an uniformly mixing condition but not strong). The obtained results appear interesting and encourage further studies varying both probabilistic laws and mechanical models of the beam (for example, layered beam with random elastic modulus varying along the section height and enriched beam model). Eventually, the analysis has to be extended to two-dimensional dynamic problems.

Acknowledgments
Authors gratefully acknowledge the support received from the Italian Ministry of University and Research, through the PRIN 2015 funding scheme (project 2015JW9NJT -Advanced mechanical modelling of new materials and structures for the solution of 2020 Horizon challenges).