{'fraction_non_alphanumeric': 0.029682737600443335, 'fraction_numerical': 0.009905791077860904, 'mean_word_length': 4.2871268998351955, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modelling Washboard Road: from experimental measurements to linear stability analysis

When submitted to the repeated passages of vehicles unpaved roads made of sand or gravel can develop a ripply pattern known as washboard or corrugated road. We propose a stability analysis based on experimental measurements of the force acting on a blade (or plow) dragged on a circular sand track and show that a linear model is sufficient to describe the instability near onset. The relation between the trajectory of the plow and the profile of the sand bed left after its passage is studied experimentally. The various terms in the expression of the lift force created by the flow of granular material on the plow are determined up to first order by imposing a sinusoidal trajectory to the blade on an initially flat sand bed, as well as by imposing a horizontal trajectory on an initially rippled sand bed. Our model recovers all the previously observed features of washboard road and accurately predicts the most unstable wavelength near onset as well as the critical velocity for the instability.

When submitted to the repeated passages of vehicles unpaved roads made of sand or gravel can develop a ripply pattern known as washboard or corrugated road. We propose a stability analysis based on experimental measurements of the force acting on a blade (or plow) dragged on a circular sand track and show that a linear model is sufficient to describe the instability near onset. The relation between the trajectory of the plow and the profile of the sand bed left after its passage is studied experimentally. The various terms in the expression of the lift force created by the flow of granular material on the plow are determined up to first order by imposing a sinusoidal trajectory to the blade on an initially flat sand bed, as well as by imposing a horizontal trajectory on an initially rippled sand bed. Our model recovers all the previously observed features of washboard road and accurately predicts the most unstable wavelength near onset as well as the critical velocity for the instability.

I. INTRODUCTION
When submitted to the repeated passages of vehicles, gravel and sand roads can develop a ripply pattern known as washboard or corrugated road. Not only is the bumpiness of the track a disturbance to drivers, but it also causes a loss of adherence and control and is therefore a real hazard. A similar phenomenon occurs on train, tramway or metro tracks and is known as rail corrugation. Due to its obvious economic importance, rail corrugation, which is caused by wear or plastic deformation, has been extensively studied as reported in the engineering literature [1] through experiments [2], field work [3] and theoretical and numerical analysis [4][5][6][7][8][9]. The formation of a washboard road over a sand bed was first studied experimentally in [10], and theoretical models have been proposed in [11,12] but include an ad-hoc diffusion term whose physical meaning remains unclear. Recently, it has been shown that a washboard road can develop when a simple inclined blade or plow (instead of a rolling wheel) is dragged along a sand road [13,14]. The formation of an instability during a unique passage of a plow under its own weight over an initially flat surface has also been studied in viscous or viscoplastic fluids [15,16].
The aim of the present paper is to derive a model for the washboard instability caused by the repeated passages of a plow over a sand bed, from experimental mea-surements of the forces acting on the plow. Previous work has focused on force measurements in the simpler case of a blade at constant altitude steadily dragging a mound of sand over a flat sand bed [17]. Herein, we extend these results by using force measurements to probe the mechanical response to an oscillatory excitation. These measurements are then used as the basis for a linear stability analysis of the washboard road instability.
The paper is organised as follows : Section II gives a list of previous experimental observations and results that a model should recover. Section III presents the experimental methods while section IV introduces the framework and assumptions of the linear stability analysis. The relation between the trajectory of the plow and the shape left in the sand bed after its passage is discussed in section V. The expression of the lift forces acting on the plow is discussed in section VI. Finally the predictions of the linear stability analysis are compared to experimental results in section VII.

II. PREVIOUS RESULTS
Our previous work has established a number of characteristic features which a model for washboard road ought to reproduce. (i) There exists a critical velocity, v c , below which the sand bed remains flat when perturbed and above which any irregularity will develop into a regular rippled pattern. (ii) This critical velocity increases with increasing mass of the plow or wheel, following a power law. (iii) For velocities greater than the critical velocity the rippled washboard road pattern appears only gradually over several passages of the plow. While the initial wavelength, λ keeps a finite and well-defined value, the amplitude of the ripples grows continuously from zero to a saturated value.
In the present paper we focus on the onset of the instability, i.e. on velocities greater yet close to the critical velocity. The saturation and coarsening of the pattern remain to be studied in further detail.
Previous work has shown that in the case of washboard road caused by a plow, compaction of the sand bed has not been observed nor seems to play any role [14], contrary to the case of a rotating wheel. This may be due to the fact that when using a plow the granular material is constantly shuffled at every passage of the plow. Therefore compaction will be neglected throughout the present stability analysis. Note that this assumption would be highly debatable had a wheel been used.
For the range of velocities and masses of the plow used we found that the "horizontal" length of the plowed material, L 0 (typically a few centimeters), remains smaller than the wavelength of the washboard road pattern, λ typically ten centimeters. Its influence is briefly discussed in section VI C but unless otherwise mentioned it will be neglected.

III. EXPERIMENTAL SETUP
The experimental setup consists in a circular track of length 5 m (average diameter of 80 cm), width 25 cm, filled with a 20-cm high layer of granular material (figure 1). The granular material used is sand-blasting sand of size ranging from 200 µm to 400 µm. The plow consisting of a flat blade made in PVC (inclined at 45 • ) of width 14.5 cm is dragged at a constant horizontal velocity v, ranging from 0.1 m s −1 to 2 m s −1 . These values of the speed are high enough to produce a continuousflow regime and low enough to avoid a gaseous regime. The plow is either attached to an arm whose end is free to rotate (figure 2a) or mounted on a translation stage (figure 2b).
In the former case (figure 2a), the length of the arm is 50 cm while the amplitude of the ripples near onset is typically of the order of a few millimeter. Therefore, we will consider that the arm remains nearly horizontal and the pendulum motion of the arm is neglected, so that any change in either the horizontal speed or in the inclination of the plow due to a change in altitude is neglected. The vertical dynamics of the plow is therefore governed only by its own weight and the lift force caused by the plowed material, while its horizontal motion is imposed by a motor and the washboard road pattern may develop. Compared to most engineering work [1] our system is simplified in that the plow has no tire nor suspension. Movies in the supplementary material show both the steady-state washboard instability on real-time and the growth of the ripples through stroboscopic images. The vertical position of the plow is recorded using a magnetic angle sensor (ASM-PRAS1) placed on one end of the arm. Unless otherwise mentioned, prior to any experiment the sand bed is made flat by dragging (over several tens of rotations) a gradually rising vertical blade around the track.
In the second case (figure 2b) the plow is rigidly mounted on a vertical translation stage (5 µm accuracy), while its horizontal velocity remains imposed. Two force sensors (Testwell KD40S) are used to measure the vertical lift force acting on the plow. In this second case the trajectory of the plow is imposed through a NIUSB6259 card. The sand bed profiles prior to and after the passage of the plow are measured using two laser telemeters (optoNCDT 1302 from Micro-Epsilon not shown on figure 2b) of 0.02 mm accuracy. Obviously no washboard road pattern can appear using this setup used only to probe the mechanical response of the material to an imposed trajectory.
The following variables are defined : x is the horizontal position of the lower tip of the plow (set to zero at time t = 0), y(x = v(t − (n − 1)T )) its vertical position where n is the number of passages and T the duration of one passage, h n (x) is the profile of the sand bed after the n th passage of the plow (h 0 being therefore the initial profile prior to the experiment) (see figure 2).
Note that since the horizontal velocity is imposed there always exists a correspondence between the position x and the time t. Any variable can therefore be differentiated with respect to x or t, although for convenience h n will be expressed as a function of x and y as a function of t.
The sand bed being largely thicker than the amplitude of the ripples near onset it can be considered infinitely deep (having halved the thickness of the sand bed in our experiments showed no noticeable effects). The reference of the vertical variable y and h n is therefore arbitrary.

IV. FRAMEWORK OF THE LINEAR STABILITY ANALYSIS
The washboard road pattern emerges from the interaction between the plow and the sand bed. The trajectory of the plow is governed by Newton's second law where the vertical lift force, f l , acting on the plow may depend on the plowed mass M , the vertical position of the plow y, and the sand bed profile h and their derivatives. Similarly the shape of the sand bed depends on the trajectory of the plow, while the new profile left after the n th passage, h n (x), may depend on the lift, the previous profile h n−1 (x) and the trajectory of the plow y(x).
The plow position y is governed by its dynamics (weight and lift force) while an erosion/deposition law reshapes the bed profile h (compaction being neglected). Although these two processes are coupled and simultaneous a convenient way to represent the coupling is the schematics of figure 3 commonly used in the engineering literature [4,5,9]. In the following section the reshaping of the sand bed is discussed whereas the dynamics of the plow is studied in section VI. Figure 3. Diagram representing the coupling between the dynamics of the plow and the reshaping of the sand bed due to erosion and deposition.
As mentioned in the introduction the present paper proposes a linear stability analysis near the onset of the washboard road instability. It will be shown that although nonlinearities may be needed to account for the saturation and coarsening of the ripple pattern, a linear model accurately reproduces the experimental data. The linear response is probed using sinusoidal excitations (of either the vertical position of the plow, or of the initial bed thickness) of wave-number k, or corresponding angular velocity ω = vk. As mentioned above, the origin of the vertical axes is arbitrary and does not play any role. Since compaction is here neglected the average height and thickness profile is constant and can be set to zero. Therefore, the system can be described using complex variables indicated by an under bar.
where M is the plowed mass on time-average, and similarly y = A y e iωt and h = A h e ikx . As a reminder, based on previous experimental observations the four following assumptions are made. (i) A 2-dimensional model is suitable, (ii) compaction is neglected and (iii) the length of the plowed material (L 0 on figure 2) is smaller than the wavelength of the pattern and (iv) near the onset of the instability studied here the plow remains in contact with the sand bed.

V. RESHAPING OF THE SAND BED
In rail corrugation, the time evolution of the bed profile can involve wear or plastic deformation [1] and is governed by the normal and tangential forces acting between the wheel and the track, as well as by the trajectory of the wheel y(t) and the previous shape of the track, h n−1 (x), these three variables being coupled. Here in the case of a plow running over a sand bed the sand is simply eroded and redeposited (compaction being neglected).
In order to investigate the relation between h n (x) and the other variables the plow is dragged at constant horizontal velocity v while a vertical sinusoidal trajectory of angular velocity, ω, (and corresponding wavenumber k = ω/v) is imposed to the plow mounted on the translation stage : y(t) = A y e iωt . The profile, h n (x) is then recorded. Figure 4 shows a typical example of the outcome for v = 0.5 m s −1 , A y = 0.6 mm and k = 21 m −1 (corresponding wavelength λ = 30 cm).
x (m) y, h (mm) The first result is that there is no phase-lag between the imposed vertical position, y, and the profile left after the passage of the plow, h n (see Figure 4a). Our experiments have shown that this remains true for the range of values of v, A y and k of interest in this paper (near the onset of the instability, v is ranging from 0.2 to 2.0 m s −1 and k is ranging from 10 to 125 m −1 ).
The major result is that the amplitude of h n is slightly smaller than that of y, meaning that the plow leaves a footprint in the sand which is similar yet not identical to its trajectory. The ratio between the two amplitudes, χ, is plotted as a function of wavenumber and velocity (in figure 4 b and c). The ratio χ is clearly independent of the wavenumber and depends very weakly on the velocity. The ratio χ can be seen as an imprinting efficiency ratio. It seems intuitive that it should tend toward one for low velocities since when trying to carve a shape in a sand bed, one should be as delicate (and slow) as possible. Note however that over the range of velocities of interest χ remains rather constant at χ ≈ 0.9 and any dependence on v will be neglected in the following. Moreover, the ratio χ is also found to be independent of the amplitude of the imposed trajectory, A y , for the range of small amplitudes studied (A y < 5 mm). In the present case where the trajectory is imposed by the translation stage the profile left after the n th passage of the plow, h n , is found to be independent of that prior to this passage, h n−1 . In the following, we will therefore consider that the profile left after the n th passage of the plow, h n (x) is simply given by : where χ = 0.9±0.05. Although this value is close to unity it will emerge that the difference between the trajectory of the plow and the shape left in the sand bed is a key element of the washboard road instability.

A. Jerk equation
In a situation where the plow is free to move vertically its dynamics is simply governed by the lift force acting on the plow f l (its horizontal position being imposed) : Our recent study [17] has shown that in a uniform motion over a flat bed (constant horizontal velocity and no vertical motion), the lift force is simply proportional to the plowed mass M and is independent of the velocity : f l = µM g, where µ is an effective solid friction coefficient which depends solely on the inclination of the plow (or angle of attack), and µ = 0.56 for an inclination of 45 • . The goal of this section is to extend our previous analysis to non-uniform motion over a rippled sand bed.
From mass conservation (and having neglected the length L 0 , see figure 2) the rate at which mass is gained or lost is simply given by : where ρ is the 2D density of the sand (ρ = bulk density of sand × packing fraction × plow width). In a steady state (horizontal or sinusoidal motion) equations (2) and (3) simply give a relation between the average plowed mass and the mass of the plow : On an initially flat sand bed (h n−1 (x) = 0) differentiating equation (2) and using equations (1) and (3) yields : where K = (µgρχ/(mv 2 )) 1/3 . This equation known as a jerk equation governs the dynamics of third-order oscillators and charged particules in motion in their own electromagnetic field. Equation (5) has oscillatory solutions (of wavenumber √ 3K/2) but with an exponentially growing amplitude along the track during one unique passage over a flat bed, meaning that the road is always unstable. This is in total contradiction with the experimental observation of a critical velocity, at least in the range of masses and velocities under study. Therefore additional dissipative terms are needed in equation (5). Note however that the wavelength predicted by the jerk equation (5) is rather close to that observed near the onset of the instability (i.e. the most unstable mode for a velocity close to the critical velocity). Indeed, for v = 0.8 m s −1 and m = 0.25 kg, equation (5) predicts a wavelength λ = 4π √ 3K = 0.41 m close to the experimental wavelength of the pattern for similar values of the parameters [13,14].
To the first order described in this section, the lift force depends on the plowed mass M , i.e. on h n and h n−1 as seen from equation 3. The simplest mathematical form is to include the derivatives of h n and h n−1 , whose physical meaning will be discussed below : . Under the linearity assumption used here it will be assumed that the effects of all further terms are additive. Section VI B is devoted to the study of the case h n−1 = 0 whereas section VI C will focus on the case where h n = 0 and h n−1 (x) = A h e ikx .

B. Lift force over an initially flat bed
This section aims at determining the dependence of the lift force on h n and h ′ n in the case of an initially flat and horizontal bed, h n−1 = 0. The plow is mounted on the translation stage and the vertical force is recorded as the vertical position is imposed.
An initial mass of sand M is plowed by lowering the plow in the sand bed. From this position (at t = 0) a sinusoidal motion is imposed to the plow : y = A y e iωt = A y e ikx . The sand bed is flattened at every rotation using a vertical blade attached to another arm which erases any pre-existing profile. This allows the data to be averaged over typically 30 oscillation periods. The lift force reads : where b and c are real and positive coefficients whose dependence on the average plowed mass M and on the velocity is key to the stability analysis. Note that in the steady regime studied here equation (6) could be expressed as a function of h n and h ′ n using equation (1). Although the two additional terms by and cẏ are introduced as the simplest linear extension of the lift force, their physical meaning is easily understood. The first additional term can be seen as a restoring force acting on an intruder, adding an extra force which is proportional to the penetration depth y. Clearly when pushing an intruder into the sand bed (y < 0), the additional force is positive, hence the negative sign in −by. Such linear restoring forces have been reported in the literature [18][19][20][21].
The second term may be seen as a contribution of the vertical penetration speed (or equivalently of the local slope of the profile when h n−1 = 0). Again the additional force should be positive when penetrating into the sand bed (ẏ < 0), hence the negative sign in −cẏ.
In order to measure the coefficients b and c and to study their dependence on the average mass M and velocity v, the transfer function H is defined : H = f l /y. Equation (6) yields the following band-stop filter expression : Equation (7) predicts that there exists a minimum in |H| when H is real and negative (ψ = arg(H) = π), i.e. for ω = ω c ≡ (µgρχv/c) 1/2 . Therefore the coefficient c can be computed from the value of ω c while the coefficient b is simply the minimum reached by |H|. In the framework of equation (6) measuring min(|H|) and ω c is sufficient to experimentally determine the values of the two coefficients b and c. The experimental data of the absolute value has a higher noise level (about 15%) than the argument (about 5%) and hence, the value of ω c was measured from the argument ψ(ω c ) = π (dashed line on figure 5).
The transfer function H was measured for various plowed masses, M , and velocities, v, and a typical example is shown on Figure 5 (v = 0.3 m s −1 and M = 0.58 kg). The experimental data shows the features predicted by equation (7) (minimum in |H|, reached for H real (ψ = π), ψ ranging from π/2 to 3π/2) and both the absolute value and the argument are well fitted by equation (7) (solid lines). For the small amplitudes imposed here (typically A y < 5 mm) we found no dependence on the excitation amplitude A y . These results validate the linear form of equation (6).
The values of ω c computed from the transfer function are plotted in figure 5c as a function of v for an averaged plowed mass M = 0.58 kg. It appears that ω c is a linear function of the velocity v : The physical meaning of this dependence remains to be explained. Still this empirical expression along with measurements of min(|H|) for various plowed masses and velocities allows to plot the two coefficients b and c as functions of M and v with a typical uncertainty of 10% ( figure 6). The following observations emerge : over the range of parameters of interest, b does not show any systematic dependence on M and is a linear function of v ( Fig. 6a and 6b) and the following expression will be used b = B 0 v, with B 0 = 560 kg m −1 .s −1 . The coefficient c appears to be proportional to the mass M (solid line of slope one in figure 5c) whereas its dependence on v has already been determined from the empirical expression of ω c (solid line on figure 5). Overall the following expression emerges :

C. Lift force over a rippled sand bed
The previous section has established empirical expressions for the contribution of the derivatives of y(t) (or equivalently h n (x)) to the lift force, f l . In this section the contribution to f l of a pre-existing profile prior to the n th passage of the plow is studied by shaping a sinusoidal profile (h n−1 = A h e ikx = A h e iωt ) over which a horizontal trajectory is imposed (y = 0, and hence h n = 0) with a blade plowing an average mass of sand M .
Similar additional terms to the lift forces are expected although with an opposite sign. Indeed a bump in the sand bed (h n−1 > 0) will create an extra positive force as will a positive slope h ′ n−1 : whereb andc and coefficients whose dependence on the averaged plowed mass M and on the velocity v has to be determined. It is expected that the role of y in equation (6) and h n−1 in equation (8) are symmetrical and we propose thatb = b andc = c. A new transfer function can be defined as G = f l /h n−1 and is plotted as a function of ω on figure 7 (for M = 1 kg and v = 1.5 m s −1 , one can note that v is high compared to the previous experiment, this is to prevent the effect of the length L 0 and make sure that λ ≫ L 0 ). The solid lines correspond to the predictions of equation (8)

A. Exponential growth rate
Having empirically established the mass and velocity dependence of additional terms in the lift force in section VI using a plow whose trajectory is imposed, the equation of motion of a plow free to move vertically can be given by combining equations (2), (6) and (8) under the assumption of additive forces made above : Differentiating equation (9) and using equation (1) the growth rate of the instablitiy, σ = ln |hn| |hn−1| , can be expressed : where α = µgρ/(mv 2 ), β = B 0 /(mv) and Experimentally this growth rate is measured by preparing an initially rippled pattern (similarly to section VI C) and instead of imposing a vertical trajectory to the plow it is left free to move vertically. After several rotations of the plow the ripples either vanish or increase depending on the wavelength of the pattern and on the velocity of the plow. Figure 8 shows two typical examples of these behaviors. During the first rotations (at least 10), the amplitude of the ripples can be fitted to an exponential curve and the growth rate is calculated from the exponental fit.
Using the values of the coefficients B 0 and C 0 determined experimentally in the previous sections, the growth rate σ can be computed from equation 10 as a function of the wavenumber k (or corresponding wavelength λ) for various values of the average mass M and velocity v. The predictions of the model are in good qualitative and quantitative agreement with experimental results and show an important feature : for low velocities, the growth rate remains negative for all wavelengths, meaning that any perturbation will be gradually eroded, whereas for high velocities there exists a positive maximum corresponding to the fastest growing mode or the most unstable wavelength. of a critical velocity but also predicts the fastest growing mode : the theoretical and experimental most unstable wavelengths are in excellent agreement as shown in figure 9.
It should be noted that there is some degree of uncertainty as to the asymptotic behavior of σ since the expressions for the additional terms in equation (6)  The critical velocity v c for which the first positive value of σ is observed depends only on one parameter : the average plowed mass M , i.e. the mass of the plow m = µ M . Figure 10 (top) shows the values of v c inferred from our model for a plow mass ranging from 100 g to 2 kg. The critical velocity follows a power-law as a function of this mass : v c ∝ m 0.35 , as found in previous experiments [14], yet with a lower exponent (v c ∝ m 0.25 ). Although it is satisfactory to recover a power-law, the origin of the discrepancy in the exponent remains unclear. However the relevance of a power-law over only one decade is debatable. Again expanding the predictions of equation (10) to a wider range of masses, m, would be questionable since the expression of the lift force was determined in a somewhat narrow range of parameters near the onset of the instability.
Finally the value of the predicted most unstable wavenumber, k th , can be compared to that measured experimentally, k exp . The latter was not measured from the plots of σ but instead was simply deduced from the pro-file of the sand bed after a few rotations. These values of k th and k exp depend on the mass of the plow m as well as on the velocity, v (clearly for v > v c (m)). For simplicity, they are plotted not as a function of m and v but instead k exp is plotted as a function of k th on figure 10 (bottom). The agreement between the predicted and experimentally measured most unstable wavenumber is excellent for all values of the mass and velocity as shown by the solid line of slope one.

VIII. CONCLUSION
We have presented a linear stability analysis based on experimental measurements of the lift force acting on a blade plowing a mound of sand on a sand bed. From previous experimental work [17] which focused on the case of a steady plow over a flat surface an equation for the dynamics of the plow was derived. However we found that further terms were needed to recover the observed features of the instability. The empirical expressions of these additional terms were obtained by probing the mechanical response of the system to sinusoidal excitations, both in the trajectory of the plow and in the initial profile of the sand bed. From this the growth rate of a sinusoidal perturbation was calculated and we showed that above a critical velocity, which depends on the mass of the plow, it displays a positive maximum corresponding to the most unstable or fastest growing mode. The critical velocity as well as the most unstable wavenumber were computed and showed excellent agreement with experimental measurements.
The expression of the lift force was obtained for a range of parameters near the onset of the washboard road instability. Hence the model recovers all the experimental observations near the threshold of the instability but fails to include the saturation of the pattern observed in experiments. Further measurements would be necessary and a thorough modelling will require to include nonlinearities.
Finally it would be interesting to extend the present study to the case of a rolling wheel. However this remains challenging since on the one hand the transport mechanism is clearly different and may involve plastic deformation and on the other hand compaction is expected to play a crucial role.