Analysis of Wave Propagation in a Discrete Chain of Bilinear Oscillators

The process of wave propagation in the discrete chain of bilinear oscillators subjected to several types of external harmonic excitation has been analysed. The phenomenon of sign inversion of the displacement is observed for tensioncompression excitation. The solution for wave propagation in a continuous 1D bimodular rod is developed and the numerical 10 results are compared.


Introduction
In this paper, we analyse the process of wave propagation in a chain of bilinear oscillators -discrete masses connected by springs having different stiffnesses in tension and compression.Due to their simplicity, discrete chains of bilinear oscillators have often been used in the problems related to non-linear vibrations of mechanical systems, such as vibrations in suspension bridges (De Freitas et al., 2004) and in the systems with the so-called fatigue cracks (Rivola and White, 1998;Ohara et al., 2007;Peng et al., 2007).Bilinear oscillators were also used in mathematical modelling of seismic isolation systems (Skinner et al., 1993;Chang et al., 2002).Layered rocks and rocks with a single set of open fractures obviously exhibit bilinear properties whereby the modulus in compression is higher than the modulus in tension due to the closure of interlayer gaps and fractures in compression.
The behaviour of the bilinear oscillators has been recently studied in Dyskin et al. (2012Dyskin et al. ( , 2014) ) and Guzek et al. (2016) for a limiting case of an infinite stiffness in compression.However, a general case of a discrete chain of bilinear oscillators has never been studied with respect to the mechanical wave propagation, which is why it has been decided to nu-merically investigate the response of the bilinear system that could represent a continuous bimodular medium.We focus on a conservative system; for the effects of damping in bilinear oscillators see Holmes (1983), Natsiavas (1990a, b), Liu et al. (2015), Dyskin et al. (2012), Klepka et al. (2015), or Guzek et al. (2016).
The purpose of the present work is to study the response of a discrete system of bilinear oscillators loaded by an external harmonic force, especially for the case of the large difference between spring stiffnesses in tension and compression.In order to compare the chain of bilinear oscillators with its homogenized counterpart, we also considered a continuous 1-D bimodular rod and developed a solution for its wave equation.In doing so, we will not restrict ourselves to small difference in stiffnesses, thus providing a more general analysis than the ones presented in Naugolnykh and Ostrovsky (1998) and Gavrilov and Herman (2012).

Mathematical formulation
We consider an infinite chain of masses and bilinear springs, where masses M are supposed to be identical, springs have the length L s and the stiffness described in the following formula: Here, U (X, T ) is the displacement, K 0 is the average stiffness of the bilinear spring, a is the stiffness ratio, and U is the difference of displacements of two adjacent masses, that is the displacement of each spring.The mass-spring chain is fixed at the right end and loaded by an external force F (T ) from the left end (Fig. 1).By introducing the Lagrangian L = E k − V , which is the difference between the total kinetic and strain energy of the system where we obtain the governing equation of the longitudinal motion of ith mass: where i = 1 corresponds to the first mass from the left end.
Since loading is applied to the left end, it follows that Here and in what follows, we consider harmonic external loading of the type F (T ) = F 0 sin ( T ), where F 0 denotes any multiplier in front of the harmonic function and denotes the external excitation frequency.We rewrite the equation of motion Eq. (3) in terms of dimensionless displacement u and dimensionless time t: where 0 is the basic frequency of the bilinear oscillator 0 = K 0 M , and c is the sound velocity in the discrete chain Here, f i and k i are the dimensionless forces and stiffnesses of the springs, respectively: Without loss of generality, we adopt that the springs are stiffer in compression and obtain higher dimensionless compressive stiffness k c = 1 + a and lower tensile stiffness k t = 1 − a for a > 0.
3 Mechanical parameters of the discrete mass-spring chain All the numerical results presented in the paper are obtained for the dimensionless parameters listed in Table 1.

Impulse harmonic excitation
In the analysis of wave propagation caused by initial excitation, simple harmonic or sinusoidal waves are of substantial interest.Due to its simplicity, let us analyse the case of a harmonic impulse first.The external loading is subjected to the left end of the chain and is described as follows: where H (t) is the Heaviside function.An explicit Runge-Kutta method with the time step t = 10 −3 is used for solving the system of N bilinear ordinary differential equations, Eq. ( 6).

Compression-tension harmonic impulse
The analysis will start with the positive sign in Eq. ( 8)in other words when compression is followed by tension.
Knowing that in the bimodular chain the compressive wave travels with a higher speed than the tensile one, one would expect the distance between compressive and tensile zones to increase with time.Figure 2 depicts the displacement field along the bilinear chain against the mass number (integer value of the coordinate) at different time moments.Since the initial load is applied from the left end of the chain, positive displacement corresponds to compression and negative displacement to tension.As expected, looking at the zones with zero deformation, i.e. horizontal regions with nearly constant positive displacement, makes it clear that the gap between compressive and tensile fronts increases with time.This phenomenon always takes place when the external excitation corresponding to a faster wave is followed by a slower one.

Tension-compression harmonic impulse
The second type of loading is described by Eq. ( 8) taken with a negative sign.This case is of considerable interest due to the fact that excitation corresponding to a slower wave speed is followed by a faster one.In this case, the faster wave front  catches up with a slower front, which leads to unusual behaviour of displacement observed in Fig. 3. Soon after the collision t * between the compressive and tensile wave fronts, the displacement gradually changes from negative to positive implying that, although a tensile impulse is applied first, the system undergoes compressive displacement after the collision.Hereafter this phenomenon is referred to as sign inversion.
The collision is defined by the time when the fast-moving wave front with negative gradient touches the wave front with the slow-moving positive gradient and is determined by the following equation: which gives t * ≈ 55 for this particular case.

Energy conservation
As an additional check on the accuracy of the numerical solution, the integral total energy E = E k +V has been calculated for the entire system of masses and bilinear springs.As seen in Fig. 4, soon after the impulse loading is applied (that is energy is added to the system), the total energy reaches its maximum and remains constant throughout the entire solution.

Continuous harmonic excitation
The second type of excitation considered here is a continuous external loading applied to the left end of the chain: As in Sect.4, two cases will be considered: compressiontension and tension-compression sequences.Obviously, in the case of continuous excitation, the difference between these two cases is the difference in the initial phase.

Compression-tension harmonic excitation
Numerical solution for displacement u at different times versus horizontal coordinate x is presented in Fig. 5.It may be observed that, due to the tensile stiffness being lower than the compressive one, the displacements close to the left end of the chain decrease with time, implying that the left part of the chain undergoes increasing tensile displacements.and 6 suggests that the numerical solution exhibits little sensitivity towards the excitation phase.This is easy to interpret given that the compression-tension and tension-compression excitations are just different phase shifts of the same continuous harmonic excitation.In this section, we want to compare the numerical results for the discrete chain of bilinear oscillators with its homogenized counterpart, a continuous 1-D bimodular rod, subjected to the same boundary conditions.In order to ensure whether the discrete chain with the given parameters can be considered as a continuum, let us estimate the dimensionless wave length λ:

Tension-compression harmonic excitation
The obtained wave length λ is much greater than the spring length L s , assumed to be equal to 1 (see Table 1), which is why the continuum approximation becomes possible.This will be done in order to check whether a numerical solution of the corresponding continuous problem can be accurate.
The wave equation for a 1-D rod made of a bimodular material reads where U , X, and T are the displacement, coordinate along the rod, and time respectively, ρ is the specific mass, E 0 is the "average" Young's modulus, and e is the difference between Young's moduli in tension E t = E 0 − e and in compression In the dimensionless form, Eq. ( 11) reads where T 0 , and T 0 is the duration of the applied external impulse.
The analytical solution for the compression-tension excitation of the frequency ω = 1 described in Sect.6.1 has been derived in Gavrilov and Herman (2012) and was later extended for the arbitrary dimensionless excitation frequency in Kuznetsova et al. (2016).
Numerical results are obtained by solving Eq. ( 12) using the explicit central difference scheme.To match the results obtained for the discrete chain of bilinear oscillators, spatial and time discretization is chosen to be the same ( x = 1 and t = 10 −3 , as in Sect.4) with all other parameters used being from Table 1.

Compression-tension harmonic impulse
The displacements u(x, t) for various times are plotted in Fig. 7, which includes the analytical solution (bold dashed line) and numerical results for the discrete chain and the 1-D bimodular rod (solid and dash-dot lines, respectively).One can observe that the three approaches show good agreement at the wave front and a slight discrepancy behind the wave front, which is typical for the second-order finite difference schemes (Kukudzhanov, 2013).

Tension-compression harmonic impulse
As the analytical solution does not exist for this case, only numerical results are presented.Figure 8 shows the numerical results for the displacements u(x, t) for the discrete chain (solid line) and the 1-D bimodular rod (dash-dot line).It is interesting to note that with all parameters being equal, the discrete chain generally exhibits lower displacements throughout the entire solution.This discrepancy may be explained by the insufficiently small spatial step for the rod since it is assumed to be equal to the length of the springs in the discrete chain which equals 1.
The direct comparison with the numerical solution (with step over x equal to 1) of the partial differential equation corresponding to the continuous rod demonstrates that the solu-tions of Eqs. ( 6) and ( 12) are close.This indicates the possibility to solve the corresponding partial differential equation numerically despite the presence of discontinuous coefficients.

Conclusions
This paper analysed the response of the discrete chain of bilinear oscillators and the bimodular rod subjected to several types of external harmonic excitation.To the best of our knowledge, wave propagation in bilinear oscillators with large stiffness ratio has never been considered before.The phenomenon of sign inversion of the displacement consisting of the gradual change of displacement sign for extended times is observed for both the discrete chain and the bimodular rod under the tension-compression impulse.It suggests that the collision between the two wave fronts corresponding to compression and tension phases has a considerable effect on the dynamic behaviour of the bilinear material.
It is anticipated that this observation may play an important role in geophysical and exploration applications, making it possible to detect bilinearity and thus obtain additional information on the composition and structure of the Earth's crust.

Figure 2 .
Figure 2. Displacement u (x, t) at different time moments versus the horizontal coordinate x for the compression-tension harmonic impulse.

Figure 3 .
Figure 3. Displacement u (x, t) at different time moments versus the horizontal coordinate x for the tension-compression harmonic impulse.

Figure 6
Figure 6 represents displacement u along the discrete chain at different time moments.Comparison of Figs.5and 6suggests that the numerical solution exhibits little sensitivity towards the excitation phase.This is easy to interpret given that the compression-tension and tension-compression excitations are just different phase shifts of the same continuous harmonic excitation.

Figure 4 .
Figure 4. Integral total energy in the discrete chain with respect to time: for compression-tension (the solid line) and tension-compression (the dashed line) harmonic impulses.

Figure 5 .6
Figure 5. Displacement u (x, t) at different time moments versus the horizontal coordinate x for the compression-tension harmonic impulse.

Figure 6 .
Figure 6.Displacement u (x, t) at different time moments versus the horizontal coordinate x for the tension-compression harmonic impulse.

Figure 7 .
Figure 7. Displacement u (x, t) at different time moments versus the horizontal coordinate x for the compression-tension harmonic impulse: analytical solution (dashed line), numerical solutions for the discrete chain (solid line), and the bimodular rod (thin dash-dot line).

Figure 8 .
Figure 8. Displacement u (x, t) at different time moments versus the horizontal coordinate x for the tension-compression harmonic impulse: numerical solutions for the discrete chain (solid line), and the bimodular rod (dash-dot line).