Waves in Materials with Microstructure: Numerical Simulation

Results of numerical experiments are presented in order to compare direct numerical calculations of wave propagation in a laminate with prescribed properties and corresponding results obtained for an effective medium with the microstructure modelling. These numerical experiments allowed us to analyse the advantages and weaknesses of the microstructure model.


BACKGROUND: MICROSTRUCTURE MODELLING
It is well known that the presence of a microstructure leads to the dispersion of waves propagating in the medium. The most general dispersive wave equation in one dimension is presented by Engelbrecht et al. (2005), based on the Mindlin theory of microstructure where u is the displacement, c is the elastic wave speed, ρ 0 is matter density, A, B,C, and I are coefficients defined later; subscripts denote derivatives. More particular cases of the dispersive wave equation can be found in papers by Santosa and Symes (1991), Maugin (1995Maugin ( , 1999, Wang and Sun (2002), Fish et al. (2002), Engelbrecht and Pastrone (2003), and Metrikine (2006).
As shown by Engelbrecht et al. (2005), Eq. (1) is equivalent to the system of two equations of motion (Engelbrecht et al., 1999) where the macrostress σ , the microstress η, and the interactive force τ are defined as derivatives of the free energy function and the quadratic free energy dependence holds Here c is the elastic wave speed, as before, A, B,C, and D are material parameters, ϕ and ψ are dual internal variables (Ván et al., 2008). Due to definitions (4) and (5), the equations of motion both for macroscale and for microstructure can be represented in the form which includes only a primary internal variable where I = 1/(R 2 D) and R is an appropriate constant. In terms of strain and particle velocity, Eq. (6) can be rewritten as Particle velocity and strain are related by the compatibility condition which forms the system of equations for these two variables. Similarly, introducing microvelocity w as follows: which is the compatibility condition at the microlevel, we see immediately from Eqs (7) and (10) that Integrating the latter equation over x, we arrive at (with the accuracy up to an arbitrary constant) Thus, we have two coupled systems of equations, (8), (9) and (10), (12), for the determination of four unknowns: ε, v, ϕ, and w. The goal of this paper is to examine the validity of the microstructure model. This is achieved by comparison of the results of numerical simulations of pulse propagation performed by using the microstructure model with the corresponding results of pulse propagation in a medium with known heterogeneity properties (here referred to as comparison medium).
Section 2 is devoted to the description of the high-resolution wave propagation algorithm for the solution of the coupled systems of equations (8), (9) and (10), (12). The results of numerical simulations for test problems are presented for both microstructured and comparison media. In Section 3 the conclusions are given.

Local equilibrium approximation
Waves propagate in solids with a velocity of the order of 1000 m/s. The corresponding time is of the order of hundreds or even tens of microseconds, especially in impact-induced events. It is difficult to expect that the corresponding states of material points during such fast processes are equilibrium ones. The hypothesis of local equilibrium is commonly used to avoid the trouble with non-equilibrium states.
The splitting of the body into a finite number of computational cells and averaging all the fields over the cell volumes leads to a situation known in thermodynamics as 'endoreversible system' (Hoffmann et al., 1997). This means that even if the state of each computational cell can be associated with the corresponding local equilibrium state (and, therefore, temperature and entropy can be defined as usual), the state of the whole body is a non-equilibrium one. The computational cells interact with each other, which leads to the appearance of excess quantities (Muschik and Berezovski, 2004) that are introduced here both for macroand microfields Here an overbar denotes averaged quantities, Σ is the excess stress, V is the excess velocity, Φ is the excess microstress, and Ω is the excess microvelocity.
Integrating Eqs (8) and (9) over a computational cell, we have, respectively, ∂ ∂t x n+1 Here σ = ρ 0 c 2 ε, and superscripts '+' and '-' denote values at the right and left boundaries of the cell, respectively. Determining the average quantities we can construct a first-order Godunov-type scheme for the system of Eqs (8), (9) in terms of excess quantities by the finite difference approximation of time derivatives. Here the superscript k denotes the time step, the subscript n denotes the number of the computational cell, and ∆t and ∆x are the time step and the space step, respectively.

Excess quantities at the boundaries between cells
Although excess quantities are determined formally everywhere inside computational cells, we need to know their values only at the boundaries of the cells, where they play the role of numerical fluxes.
The values of excess quantities are determined from the jump relations providing the continuity of full stresses and velocities at the boundaries between computational cells (cf. Berezovski and Maugin, 2004) where [[A]] = A + − A − denotes the jump of the enclosure at the discontinuity, A ± are the uniform limits of A in approaching the discontinuity from the ± side. The values of excess stresses and excess velocities at the boundaries between computational cells are not independent. They are related to each other at the cell boundary as follows : As shown by Berezovski et al. (2008), the excess quantities following from non-equilibrium jump relations at the boundaries between computational cells correspond to the numerical fluxes in the conservative wavepropagation algorithm (Bale et al., 2003). Therefore, the numerical scheme (17)-(22) is the reformulation and generalization of this algorithm in terms of excess quantities. The advantages of the wave-propagation algorithm are high resolution (LeVeque, 1997(LeVeque, , 2002 and the possibility of a natural extension to higher dimensions (Langseth and LeVeque, 2000).

Excess quantities for internal variables
Since the systems of equations (8), (9) and (10), (12) are coupled, we need to solve the system of equations for internal variables (10), (12) simultaneously. Representing the mentioned system of equations in the form we can construct the corresponding numerical scheme similarly to the case of macromotion where, in its turn, the values of the corresponding excess quantities are determined from the jump relations at the boundaries between computational cells accompanied with the Riemann invariants conservation where a characteristic velocity for microstructure, c 1 , is introduced by C = Ic 2 1 .

Second-order corrections
We can improve the accuracy of the algorithm by introducing second-order correction terms (LeVeque, 2002), which also can be represented in terms of excess quantities both for macromotion and for microstructure The corresponding Lax-Wendroff schemes have the form for macromotion andφ for microstructure.

Results of numerical simulations
As a test problem the one-dimensional propagation of a pulse is considered. The case of the comparison medium is analysed first. In this case the specimen is assumed to be homogeneous, except for a region of length d, where periodically alternating homogeneous layers of size l are inserted (Fig. 1a). The density and longitudinal velocity of the main specimen material are chosen as ρ = 4510 kg/m 3 and c = 5240 m/s, respectively. The corresponding parameters for the material of the inhomogeneity layers are ρ 1 = 2703 kg/m 3 and c 1 = 5020 m/s, respectively. Initially, the specimen is at rest. The shape of the pulse before crossing the inhomogeneity region is formed by an excitation of the strain at the left boundary for a limited time period (0 < t < 100∆t) u x (0,t) = 1 + cos(π(t − 50∆t)/50).
The length of the pulse is λ = 100∆x. Arrows in Fig. 1 show the direction of pulse propagation. We consider different cases of the size of inhomogeneity, namely l = 8∆x, 16∆x, 32∆x, 64∆x, 128∆x. The pulse holds its shape up to entering the inhomogeneity region. After the interaction with the periodic multilayer, the single pulse is modified because of the successive reflections at each interface between the alternating layers. Alternatively, the same pulse propagation was simulated by the microstructured model (8)-(12) (Fig. 1b) with the modified sign of the internal interaction force. This means that Eq. (12) is replaced by In these calculations of pulse propagation the value of the internal length of microstructure is kept the same as the size l of the periodic layer, as well as the density and sound velocity for inhomogeneities: I = ρ 1 ,C = Ic 2 1 . The result of the comparison of pulse propagation in laminated and microstructured media is shown in Fig. 2 for the size of inhomogeneity l = 8∆x. The coefficients in the microstructure model are chosen in such a way that the location and amplitude of resulting pulses are as close as possible, which leads to the value A = 19ρc 2 in the considered case. The length of the pulse λ = 100∆x is much larger than the size of inhomogeneity, and the influence of microstructure is rather small.
Continuing the calculations, we vary the size of inhomogeneity to l = 16∆x and 32∆x. The corresponding results are presented in Figs 3 and 4. Again, the values of coefficients in the microstructure model were adjusted to keep the location and amplitude of the leading pulse A = 97ρc 2 and A = 147ρc 2 for l = 16∆x and 32∆x, respectively. As we can see, the model is able to reproduce the leading pulse change, but fails in the description of the tail.   Up to now, the size of inhomogeneity was less than the length of the initial pulse. We performed also numerical experiments where this size is comparable with the length of the pulse: l = 64∆x and 128∆x. The corresponding results are presented in Figs 5 and 6 with the values A = 665ρc 2 and A = 2059ρc 2 , respectively.
As one can see, if the inhomogeneity size is comparable with the initial pulse length, the ability of the model to reproduce the leading pulse shape is improved. The variation of the coefficient A in the computed microstructure model may be conjectured as related to the variation of the size of inhomogeneity. However, no straightforward relation is observed. As known from theoretical considerations (Engelbrecht et al., 2005), the ratio of the length of the pulse to the scale of the microstructure plays an important role. In terms of dispersion analysis, the problem is related to the difference of dispersion curves. In our calculations the inhomogeneity size is varied according to 2 n , n = 3, ..., 7, while the length of the initial pulse remains unchanged. This means that we step over from one dispersion curve (for a small ratio of the inhomogeneity size to the length of the pulse) to another (for the ratio of the order of unity). It seems that for limiting cases of the ratio, the value of A/ρc 2 is approximately proportional to 3 n , whereas this  is not the case for intermediate values (the most illustrative case corresponds to the value l = 32∆x). The corresponding value of the coefficient B is determined by means of the shift of the location of the leading pulse: B = A 2 /(ρc 2 (1 − α 2 )), where α is the value of the Courant number used in the calculation.

CONCLUSIONS
Numerical simulations of pulse propagation in a laminated medium and in a medium with microstructure were performed to compare the results in order to check the validity of the microstructure model. Material properties and the characteristic lengths used in calculations were chosen correspondingly to match both cases. The comparison of the results of numerical simulations shows that the coefficients in the improved microstructure model can be adjusted to achieve the accordance of the amplitude and location of the leading transmitted pulses in both cases. This means that the microstructure model is able to describe the wave propagation in complex media. At the same time, it is clear that the influence of the shape and orientation of inclusions need to be investigated in the framework of two-dimensional formulation.