Continuum model of wave propagation in fragmented media : linear damping approximation

Energy dissipation during wave propagation in fragmented geomaterials can be caused by independent movement of fragments leading to energy loss on their impact. By considering a pair of impacting fragments at times much greater than the period of their oscillations we show that at large time scale, the dynamics of the pair can be described by a linear viscous model with damping coefficient expressed through the restitution coefficient representing energy loss on impact. Wave 10 propagation in fragmented geomaterials is also considered at the large time scale assuming that the wavelengths are much larger than the fragment sizes such that the attenuation associated with wave scattering on the fragment interfaces can be neglected. These assumptions lead to the Kelvin-Voigt model of damping during wave propagation, which allows the determination of dispersion relationship. As the attenuation and dispersion are not related to the rate dependence of rock deformation, but rather to the interaction of fragments the increased energy dispersion at low frequencies can be seen as an 15 indication of fragmented nature of the geomaterial and the capacity of the fragments for independent movement.


Introduction
Geomaterials are often fragmented, with the fragments covering different scales.This makes it important to understand the properties of wave propagation in such geomaterials.Fragmented materials are characterized by three major features.First is the bilinear nature of contacts when stiffness in compression is considerably higher than stiffness in tension.Bilinear oscillators feature multiple resonances, both multi-harmonic and sub-harmonic (see Dyskin et al., 2007Dyskin et al., , 2010Dyskin et al., , 2012, and the literature cited therein).Furthermore, chains of bilinear oscillators possess a rich structure of main resonances also accompanied by multi-harmonic and subharmonic ones (Shufrin et al., 2012;Dyskin et al., 2014).This structure of resonances may give an explanation for the observed spectral peaks in oscillations of blocky media (Kurlenya et al., 1996a, b, c).The effect of bilinearity on wave propagation was analysed by Kuznetsova et al. (2016).
Second is the possibility of block rotations.The bending between fragments leads to elbowing of the neighbouring fragments in the course of their mutual rotations (Pasternak et al., 2006) as well as the dependence of bending stiffness on the moments (Pasternak et al., 2012;Shufrin et al., 2014).Furthermore, non-sphericity of fragments in the presence of compression creates an effect of apparent negative stiffness (Dyskin andPasternak, 2011, 2012a, b).The resulting negative Cosserat shear modulus and its influence on wave propagation were analysed by Pasternak et al. (2016).It was shown that such a medium does not possess a critical frequency; subsequently the twist wave and both shear rotational waves of all frequencies can propagate.
Third is the energy dissipation associated with the impact of blocks characterized by low restitution.The main feature of this type of dissipation is that it acts only at the neutral position of the oscillators formed by pairs of adjacent blocks.In this paper we consider only this special type of energy dissipation and its influence on P and S wave propagation assuming that otherwise the fragmented medium is unimodular (the same stiffness in compression and tension) and neglecting the effect of rotations.We assume that the wave length is considerably larger than the fragment sizes such that energy dissipation associated with scattering over interfaces can be neglected.

A pair of fragments of impact damping
Mathematically, the basic dynamic element involved in the process of wave propagation through fragmented geomaterials is a pair of neighbouring fragments, which can be modelled as a free undamped oscillator consisting of a single mass on an undamped spring complimented by a condition that velocity decreases each time the system passes through the neutral point and a second-order differential -Eq.( 1)is analysed.Herein, the velocity reduction is governed by a restitution coefficient α ∈ [0, 1]. (1) Introducing dimensionless groups X = x l , τ = ω 0 t, X 0 = x 0 l , and V 0 = v 0 lω 0 , Eq. ( 1) can be rewritten as where The solution of Eq. ( 2) prior to the first impact is a commonplace: where Analysing the second form of the solution, it is easy to determine the time of the first impact: (5) After the first impact, the velocity of oscillations reduces by the restitution coefficient α.Thus, a subsequent solution has the following form: where H (t) is the Heaviside function.
In this system, each next impact starting from the second occurs after time π from the previous impact decreasing stepwise the amplitude by α; therefore, the general solution of the system is It should be noted that the analysed problem is somewhat similar to a problem of a bilinear oscillator with an infinite stiffness in one direction (Dyskin et al., 2012(Dyskin et al., , 2013;;Guzek et al., 2016) or a ball bouncing off a solid wall (Luck and Mehta, 1993;Anagnostopoulos, 2004;Jankowski, 2006), with the same coefficient of restitution, in terms of the amplitude.The main difference is that, described by Eq. ( 1), the range of α for the former type of a problem lies in a negative region, between −1 and 0, which leads to function X being only in a positive domain.As a result, for equal and physically admissible boundary conditions and absolute values α, the odd half cycles of the solutions for those two types of problems are identical and the even half cycles are symmetrical about the axis X = 0.Both of these types of problems are demonstrated in Fig. 1.

Equivalent linear damping
Although the analysed model is relatively simple to describe, its implementation for energy dissipation during wave propagation is a challenging problem, the reason being that the boundaries of the time intervals where the system behaves linearly are not known a priori as they are influenced by incomplete restitution and hence need to be determined step by step.Therefore, an equivalent continuous damping model with effective coefficients should be selected and the relationship between it and the original model should be established.
Here, viscous damping governed by the Kelvin-Voigt model is chosen due to its simplicity and popularity.This model represents a free oscillation of a mass with damping that can be characterized in a dimensionless form by the damping coefficient ζ > 0: Equation ( 8 where Among those three solution types, only the underdamped solution is physically admissible for a comparison with the proposed model because the other two types do not intersect the axis X = 0. Consequently, the relationship between the damping coefficient and restitution coefficient is carried out using the underdamped solution. Comparing the expressions inside the sin functions in Eqs. ( 6) and ( 9), one can see that In the discrete model, the initial conditions consist of zero displacement and a given initial velocity.So, hereafter, the initial displacement X 0 is set as 0, which leads to T 1 = π .Thus, using Eq. ( 10) for the Kelvin-Voigt model, it is possible to find a relation between X (τ ) and Y (τ 2 ).For the former system, after a number of cycles N, with the passing time 2N π , the solution becomes On the other hand, for the Kelvin-Voigt model, using Eq. ( 10), one has As long as the third terms in Eqs. ( 11) and ( 12) are identical, the relationship between the damping parameters can be determined by comparison of the amplitudes of the systems: The first term of the right side in Eq. ( 13) goes rapidly to unity even for significant values of ζ .Thus, for the purpose of homogenization of the proposed model, this term can be substituted by unity (Fig. 2): It is seen from Eq. ( 14) that as ζ = 0, α = 1 (the left boundary), and when ζ → 1, α → 0 (the right boundary).The left boundary proves that the solutions for undamped oscillations (Eq.8, ζ = 0) and oscillations with full restitution (Eq.2, α = 1) are identical.The right boundary, i.e. zero restitution, can be modelled as a system where the critical damping (ζ = 1) intersects X = 0 when time tends to infinity.Also, Eq. ( 14) shows that the relationship between the damping parameters does not depend on the initial velocity of the systems or on the current time for all ζ within the range of the underdamped oscillations.Hence, taking into consideration the previous statements, Eq. ( 14) should be applied for 0 ≤ α ≤ 1.As a result, systems with both damping and incomplete restitution can be analysed by using equivalent viscous damping, reducing the complexity of the discrete problem.
When both types of energy dissipation take place, the restitution coefficient should be replaced by a damping coefficient; therefore, a reverse relationship is also important to define.It can be found relatively easily from Eq. ( 14) by selecting only positive roots.
In order to analyse the accuracy of Eq. ( 14), an example of vibration amplitudes of an oscillator representing fragmented media X(τ ) with different α and vibrations of an equivalent 3. The initial conditions for both cases are X 0 = 0 and V 0 = 1.
It is seen that between the impacts, the functions can be quite different even for high values of α.Indeed, the energy dissipation in the original system occurs at discrete times.Replacing the discrete system with a time-continuous system is equivalent to using the timescale considerably larger than the period.At this scale, the resulting damping is the same as in the original discrete system.
The dissipated energy of vibrations with the same parameters is given for both cases by the following equations, where W X is the energy dissipation function for the discrete model www.nonlin-processes-geophys.net/24/461/2017/ Nonlin.Processes Geophys., 24, 461-466, 2017 Figure 3. Vibrations for fragmented media (solid lines) and the equivalent Kelvin-Voigt model (dashed lines) for α = 0.9 (red), α = 0.6 (blue), α = 0.3 (green) and initial conditions X 0 = 0 and and W Y is the same function for the equivalent Kelvin-Voigt model.These dependencies are shown in Fig. 4.
The dissipated energy W X (τ ) cannot be approximated by between impacts, especially for small restitution coefficients, because here a continuous function is approximated by a stepwise function.Nevertheless, they approach the same values at impacts; therefore, the proposed function and the equivalent continuous function (linear damping) can be used as an approximation of the discrete one for times considerably higher than the period of free vibrations.

Wave propagation in isotropic medium with equivalent damping
Now, after establishing the large timescale equivalence of the discrete and continuous dynamics of a pair of fragments, the energy loss during wave propagation in fragmented geomaterials can be modelled by replacing the fragmented geomaterial with a visco-elastic continuum described by a Kelvin-Voigt model.The P wave velocity c p and coefficient of absorption a p are expressed by the following equations (White,  1983) and are shown in Fig. 5 for different α: where c is the P wave velocity without damping and ω 1 is the wave frequency.
It is seen that both the wave velocity and coefficient of absorption increase with frequency; however, the increase in the wave velocity becomes weaker as the restitution coefficient increases.Subsequently, the dispersion vanishes as the restitution coefficient tends to 1, i.e. the impacts are not accompanied by energy loss.It is also noteworthy that these formulae can be implemented for the S wave as well.

Conclusions
A possible mechanism of wave attenuation in fragmented geomaterials with fragment sizes much smaller than the wavelengths is the energy loss on impact of the contacting fragments with each other.The energy loss is characterized by the restitution coefficient.It is shown that the energy loss during wave propagation in such a discrete material can be modelled by an equivalent visco-elastic continuum if the characteristic times involved are considerably greater than the periods of oscillations of all neighbouring pairs of fragments.The attenuation is modelled by the Kelvin-Voigt model, its equivalent damping being expressed through the restitution coefficient and the period of oscillations of contacting fragments averaged over all pairs.For all restitution coefficients smaller than 1, the wave velocity shows a dispersion relationship, which is stronger the smaller the restitution is.The attenuation and dispersion are not related to rate-dependent rock deformation, but rather to the interaction of fragments.For that reason the effect is long-wave.Therefore, increasing damping and dispersion at low frequencies can be seen as an indication of the fragmented nature of the geomaterial and the capacity of the fragments for independent movement.

Figure 1 .
Figure 1.Vibrations for fragmented media (red solid line) and impact oscillator (blue dashed line) for absolute values of α = 0.9 and initial conditions X 0 = 0.5 and V 0 = 1.0.
Competing interests.No competing interests are present.Special issue statement.This article is part of the special issue "Waves in media with pre-existing or emerging inhomogeneities and dissipation".It is a result of the EGU General Assembly 2016, Vienna, Austria, 17-22 April 2016.