Capturing wave dispersion in heterogeneous and microstructured materials through a three-length-scale gradient elasticity formulation

1 Introduction

Dispersion of waves propagating in heterogeneous materials occurs because of nonlocal (long-range) interactions taking place in the microstructure. According to experiments, different harmonic wave components travel with different velocities [1], [2], [3], [4], [5], [6], [7], [8]. This dispersive behavior occurs because the wave length is of the same order of magnitude as the dominant heterogeneities. According to the classical elasticity theory, though, every wave number has the same phase velocity, which is reasonable for a (ideal) homogeneous medium but in contrast to the above experimental observations in (realistic) heterogeneous materials. While molecular and/or atomistic models can describe dispersive phenomena by incorporating nonlocal interactions between atoms in a discrete fashion, there exists a field of research in mechanics that has focused on how to describe these microstructural effects with continuum, rather than discrete, formulations. Evidently, there is a need for additional terms in classical continuum theory in order to capture these physical phenomena. This has given rise to the development of higher-order (or enhanced) continuum theories [9], including nonlocal elasticity [10], [11], [12] and gradient elasticity [13], [14], [15], [16], [17]. In particular, gradient elasticity theories include higher-order spatial derivatives of relevant field quantities in the equations of motions, e.g. strains, accelerations or stresses. Each higher-order term is accompanied by an internal length scale, which depends upon the underlying material microstructure, see [18] for a recent overview.

This paper presents a three-length-scale gradient elasticity formulation to capture wave dispersion observed in heterogeneous and microstructured materials. In line with previous studies from the relevant literature, in addition to the classical strain gradients [19], [20], combined (mixed) spatial-temporal derivatives are included through acceleration gradients (also called inertia gradients) [21], [22], [23], [24], [25], [26], [27]. This three-length-scale gradient formulation, recently developed by the authors [28], [29], [30], [31], is featured by an improved dispersive behavior due to the presence of two micro-inertia terms multiplying the 2nd order and the 4th order spatial derivative of the acceleration field in the governing equations of motion.

This paper complements the previous research work in the following aspects. First, the gradient model is here theoretically supported by a non-local lattice. It is well known that continualization of lattice systems with distributed microstructure leads to gradient models [31], [32]. In this regard, in [33] a two-length-scale gradient formulation with one strain gradient and only one inertia gradient has been derived by continualization of a non-local lattice model with two-neighbor interactions and with uniform distribution of micro mass. Building on this previous paper [33] but at the same time highlighting some inconsistencies noted in these earlier derivations, the lattice model is here corrected and expanded to incorporate the two inertia gradients of the present formulation. The equations of motions are obtained by applying Hamilton’s principle to the energy functionals identified by standard continualization of the discrete equations resulting from the non-local lattice model. Through these derivations, the length scale parameters can be interpreted based on geometrical and mechanical properties of the underlying microstructure. Then, the dispersive capabilities of the three-length-scale gradient elasticity model are scrutinized against experimental findings on a broad range of materials, including aluminum, bismuth, nickel, concrete and mortar.

2 Three-length-scale gradient elasticity formulation

Previous papers from the literature demonstrated that an effective class of gradient elasticity theories for capturing wave dispersion should involve mixed spatial-temporal derivatives alongside higher-order derivatives with respect to the spatial coordinates only. In physical terms, higher-order contributions should thus appear both in the stiffness and the inertia terms. The simplest (and also the most popular in this class) gradient model complying with this requirement contains two higher-order terms: one strain gradient and one inertia gradient, each accompanied by a corresponding internal length scale [13], [21]. As an extension of this model, one more inertia gradient can be added for an improved dispersive behavior, also motivated by nano-scale experimental observations on materials with microstructure [28]. The equations of motion of this three-length-scale gradient model read

where α, β, γ are three coefficients that adjust the magnitudes between three-length-scale parameters, all multiplying the length scale ℓ of the material microstructure – introduced to ensure dimensional consistency between the terms in Eq. (1). The two-length-scale gradient model discussed above is retrieved as a special case of Eq. (1) for β = 0: in this case, γℓ² may be viewed as a length scale in statics, while αℓ² represents a length scale in dynamics according to [13], [21]. Thus, βℓ⁴ is related to an additional length scale in dynamics. Through the one-dimensional format of Eq. (1), longitudinal wave dispersion for a simple thin rod of constant cross section is described by [28]

where ce≡E/ρ is the one-dimensional wave velocity of classical elasticity, ω is the angular frequency, k the wave number and c = ω/k the phase velocity. Modifying the relative magnitude of the three length scales gives rise to a wide class of dispersive behaviors, as illustrated in [28], which makes the proposed model versatile in capturing wave dispersion in a range of materials. This will be demonstrated in detail later on in this paper.

3 A non-local lattice model supporting the gradient model

The model in Eq. (1), in its one-dimensional format, can be derived from the continualization of the response of a lattice. We here provide a derivation based on a similar lattice model proposed by Polyzos and Fotiadis [33] for a two-length-scale gradient formulation (without the β term, only one inertia gradient). Assuming an infinite domain for the non-local lattice model, boundary integrals will be ignored in the below derivations.

The simple one-dimensional lattice model with two-neighbor interactions is depicted in Figure 1. It consists of a periodic mass-and-spring chain model. Lumped masses M are located at uniform particle spacing ℓ, which may be interpreted as the unit cell of a non-homogeneous material having mass density ρ_M = M/Aℓ, with A denoting the cross-sectional area. The masses are connected to each other by springs of stiffness k₁ and k₂ if the closest and second closest particles are considered; note that these springs are assumed to act in parallel rather than in series. From the Hooke’s law, the stiffness of the springs is k_i = E_iA/ℓ (i = 1, 2), with E denoting the Young’s modulus. Two values of the Young’s modulus E₁ and E₂ corresponding to neighboring and non-neighboring masses, respectively, are introduced. As an example, with reference to the microstructure of a composite material like concrete, the masses M simulate the aggregates, the springs of properties k₁ reproduces the elastic characteristics of the cementitious matrix and the ratio k₂/k₁ is a measure of the non-local to local interactions within the concrete matrix due to the particle long-range interactions – cf. Figure 1.

Figure 1:

Non-local lattice model with two neighbor interactions and physical counterpart.

Continualization is performed by replacing the displacement u_n(t) of the mass particle n located at x_n with the continuous displacement u(x, t); for the neighboring and non-neighboring particles this implies that u_n±1(t) = u(x ± ℓ, t) and u_n±2(t) = u(x ± 2ℓ, t) because x_n±1 = x_n ± ℓ and x_n±2 = x_n ± 2ℓ, respectively. In order to incorporate the effect of the micro-inertia on the axial vibration of the chain, springs k₁ pertaining to neighboring mass connections are not massless as in classical lattice models, but springs with uniformly distributed mass m₁ and density ρ_m1 (micro-volume density). In particular, the micro-mass distribution is representative of a discrete set of small masses m_e, located between two neighbor particles M, connected with each other through springs of stiffness k_e, and such that ∑me=m1 and ∑1/ke=1/k1. This micro-mass distribution is responsible for the two micro-inertia terms of the proposed gradient elasticity formulation. To further extend the previous model in [33], also springs k₂ pertaining to non-neighboring mass connections are here assumed to have uniformly distributed mass m₂ and micro-density ρ_m2; whereas in the earlier paper [33] they were supposed as massless. To compute the contribution of micro-masses m_e to the kinetic energy density of the system, in Figure 1 we introduce a local coordinate system z and the velocity of each micro-mass m_e at z, by linear interpolation, is

for the micro masses pertaining to springs k₁ and, by analogy,

for the micro masses pertaining to springs k₂. Therefore, a continualization approach for the micro-inertia is performed, where the micro-mass-spring system m_e – k_e is replaced by springs of stiffness k₁ and mass m₁ for neighbor, and k₂, m₂ for non-neighboring particles.

For a unit cell ℓ corresponding to the nth mass M, the stored potential energy density (per unit macro-volume) is

By using Taylor expansions for the terms u(x ± ℓ, t) and u(x ± 2ℓ, t) up to the ℓ⁵ terms

after some algebra and truncating terms beyond 𝒪(ℓ⁵) Eq. (5) is re-written as

wherein we omitted the space and time dependence of the displacement for brevity, i.e. u = u(x, t). In (7) we have included all the terms that contribute to the strain energy up to ℓ⁴. We point out that this is the correct version of continualization as truncation of the Taylor series should only be done once all terms of the same order of ℓ are collected. This is clearly in disagreement with the model of Polyzos and Fotiadis [33] who truncated Taylor’s expansion at an earlier stage and thus did not include particular terms in the expression of 𝒲. From Eq. (7) it can be noted that the classic terms in the non-neighbor particles are scaled by a factor 4 as compared to the neighbor particles, while the ℓ² and ℓ⁴ terms are scaled by a factor 16 and 64, respectively. Therefore, Eq. (7) can be compacted in the following form

with

where E is the effective Young’s modulus and r¯2,r¯4 are two ratios. The strain energy density of a simpler lattice model where the non-neighbor particle interactions are neglected may be retrieved from (8) by setting E₂ = 0, which leads to E = E₁ and r¯2=r¯4=1 in (9) and (10).

Similarly, the kinetic energy density for the mass M (per unit macro-volume) is equal to

The kinetic energy density for the micro-mass-spring system corresponding to springs k₁ and pertaining to the unit cell ℓ is

with the coordinate z indicating the distance of each point of the spring from its left end (Figure 1), and u˙1−(z,t),u˙1+(z,t) denoting the point velocities of the springs with end velocities u˙(x−ℓ,t),u˙(x,t) and u˙(x,t),u˙(x+ℓ,t), respectively. Considering particle spacing ℓ being very small, u˙1−(z,t),u˙1+(z,t) are obtained through linear interpolation (3). Similarly to (6), Taylor expansions are applied to the resulting terms u˙(x±ℓ,t) up to the order ℓ⁶, i.e.

On the contrary, linear expansions for the velocities were considered in the earlier model by Polyzos and Fotiadis [33], but this is somehow in contrast to the quadratic expansion in space adopted for the displacements by the same authors in the same paper. By including these higher-order terms for the velocities, extra terms are obtained that were not included in [33]. After some algebra, the following expression of the kinetic energy density of the micro-mass-spring system k₁ – m₁ truncated up to ℓ⁶ is obtained

Another difference as compared to the model in [33] is the inclusion of micro masses also in the non-neighbor springs k₂, whose contribution to the kinetic energy density, by similarity with (12), can be expressed as

In (15) u˙2−(z,t),u˙2+(z,t) are the point velocities of the springs with end velocities u˙(x−2ℓ,t),u˙(x,t) and u˙(x,t),u˙(x+2ℓ,t), respectively. Similarly to the previous case, they are computed through linear interpolation (4), then Taylor expansions are considered for the terms u˙(x±2ℓ,t) [in line with (13)], which leads to

The sum of the kinetic energy densities given in (11), (14) and (16) represents the total kinetic energy density of the rod. By comparing (16) with (14), it is noted that the ℓ² terms in the non-neighbor particles are scaled by a factor of 4, the ℓ⁴ terms are scaled by a factor of 16 and the ℓ⁶ terms by a factor 64. Thus, the following compact expression for 𝒯 can be written

where ρ = ρ_M + ρ_m1 + ρ_m2 is the total mass density of the lattice system and the r^i ratios are

We point out that several terms entering the kinetic energy density (17) were not present in the formulation discussed in [33] because of the following reasons:

1. linear expansion of the velocity field was considered in [33], whereas terms up to ℓ⁶ are included in the Taylor’s expansion (13);

2. the contribution ρ_m2 was not present in [33] because the springs k₂ were assumed massless in the previous model;

3. other terms were simply ignored in [33] because the authors truncated the Taylor series at an earlier stage, before collecting all terms with equal powers of ℓ. Consequently, the higher order terms were incomplete in [33]. As an example, there is no reason for omitting terms like u˙′u˙‴ and retaining terms like u˙′′2, because they are both accompanied by ℓ⁴ powers, thus they contribute to the overall response at the same order of approximation. Similar statements also hold for the potential energy density – compare Eq. (8) with Eq. (43) of the Polyzos and Fotiadis paper [33]. Instead, the correct truncation should be done after applying the Hamilton’s variational principle in a consistent manner, as done in this paper.

Once the energy functionals (8) and (17) are determined, the derivation of the equations of motion follows the usual procedures underlying the Hamilton’s variational principle, which yields

which is the sought equation of motion of the non-local lattice model. The micro-inertia properties of the material are characterized by the ratio ρ_m1/ρ (micro-to-macro density ratio) and by the r^i ratios that control the ρ_m2/ρ_m1 parameter, while the micro-stiffness properties are governed by the r¯i ratios that reflect the E₂/E₁ parameter. It is observed that this continualization approach provides results that are consistent with continualizing the equations of motions. For simplicity, we ignore for the moment the influence of the spring-mass system k₂ – m₂ by setting E₂ = 0 and ρ_m2 = 0 so that all the r^i and r¯i ratios are equal to one and Eq. (19) can be factorized as

which is different to the final equation presented by Polyzos and Fotiadis [33] for the reasons outlined above. The destabilizing terms entering Eq. (20) can be resolved using Padé approximations of the differential operator which, after some algebra, leads to

Comparing Eq. (1) and Eq. (21) leads to a clear, one-to-one identification of the α, β, γ constants of the proposed three-length-scale gradient elasticity formulation

which are linked to geometrical and mechanical properties of the material microstructure, represented by the non-local lattice model.

4 Comparison with experiments

Here, the ability of the three-length-scale gradient model to capture wave dispersion for a broad range of materials is scrutinized. Comparison against experimental findings on metals and alloys (aluminum, bismuth, nickel), as well as concrete (in the fresh and hardened state) and mortar is displayed in Figure 2 and Figure 3, respectively. The length scale parameters in Eq. (2) result from a non-linear least square minimization procedure. It is seen that the proposed gradient elasticity formulations are versatile, as they are able to describe a range of dispersive characteristics arising from experiments. These include: i) inflection points for metals and alloys beyond the first Brillouin zone – this experimental observation cannot be captured without the additional β term, as already demonstrated in [28]; ii) a phase velocities curve that decreases with increasing wave numbers in fresh concrete; iii) a phase velocities curve that increases with increasing wave numbers in hardened concrete and hardened mortar. Overall, the proposed gradient model captures all qualitative aspects of the experimental results.

Figure 2:

Experimental dispersion curves for aluminum, bismuth [3], [4] and nickel [5] vs. simulations with proposed gradient elasticity model.

Figure 3:

Experimental dispersion curves for fresh concrete [7], hardened concrete and hardened mortar [8] vs. simulations with proposed gradient elasticity model.

5 Conclusions

Wave dispersion occurring in heterogeneous and microstructured materials can be described through generalized continuum theories with higher-order terms. To this end, a three-length-scale gradient formulation has been used to capture experimental results concerning dispersive wave propagation in a range of materials. This formulation is featured by one strain gradient and two inertia gradients – the latter ensuring an improved dispersive behavior in comparison with a simpler model with just one micro-inertia term, widely employed in the literature. Despite being introduced by a simple phenomenological deduction based on observations of the dispersion curve of microstructured materials, this three-length-scale gradient model also has some physical justifications. This has been demonstrated here by the authors through the introduction of a non-local lattice model with two neighbor interactions, which can be viewed as the idealization of matrix stiffness and particle long-range interactions in composite materials and heterogeneous media. A similar non-local lattice model was earlier presented in the literature for a simpler gradient formulation, but the authors have pointed out in the present paper that the previous model did not contain all the terms that contributed to the response at the same order – thus, the higher order terms were incomplete. By correcting the above inconsistencies, and by properly extending the Taylor series the three length scale parameters are directly linked to geometrical and mechanical properties of the materials microstructure. Finally, the ability of the gradient formulation to capture wave dispersion characteristics arising from real experiments on a broad range of materials has been scrutinized. This has revealed a great versatility of the model to fit a variety of dispersive behaviors for different materials, spanning from metals and alloys to concrete and cementitious materials at different states.