1 Introduction

So far, many models, such as the cellular automaton (CA) models and the car-following models, have been applied to analyze the complex traffic phenomena. These models are used to describe the micro characteristics of car motion and explore the microscopic traffic flow phenomena. The car-following models show the microscopic characteristics of traffic flow by using the interaction between front and rear cars. Herman and Chandler et al. proposed the linear models^[1]. Subsequently, Gazis et al. proposed the nonlinear models^[2]. Based on the above models, remarkable improvements were achieved^[3-7]. One can deduce from Ref.[3] that Bando and others proposed OV (Optimal Velocity) model according to traffic analysis in 1995 which was then widely applied^[4]. And many scholars proposed a series of OV extended models, which considered many cars' relative position or the relative speed of the near vehicle, to enhance the stability of the traffic flow.

Most of the existing car-following models focus on the normal traffic condition^[8-12]. In fact, under ice and snowfall conditions, the visibility and friction factor decline significantly, and the traffic property and drivers' behaviors are crucially distinct from those under the normal state. Drivers will keep bigger headway and space headway, and make slighter acceleration or deceleration under ice and snowfall conditions. As a result, the car-following model is not the same as the normal one. However, there are very few car-following models which focus on ice and snowfall conditions. Hence, in this paper, we directly focus on car-following model to describe the characteristics of car motion under ice and snowfall conditions, then grasp the stability conditions of traffic flow by linear analysis, and finally verify the rationality of the proposed model by nonlinear analysis and numerical simulation.

2 Model Establishment

Bandoet al.^[3] proposed this model:

$ \frac{{{\rm{d}}{v_n}}}{{{\rm{d}}t}} = {a_1}\left[ {V\left( {\Delta {x_n}} \right) - {v_n}} \right] $

(1)

where, a₁ is the sensitivity coefficient to the vehicle in front, Δx_n=x_n+1-x_n, is the space headway, x_n and x_n+1 are the locations of the flowing car and the first car respectively, v_n is the velocity of the following car, and V (Δx_n) is the optimal velocity.

According to the observed data, Helbing and Tilch identified the parameters of OV model^[13]. It turns out that there are excessive acceleration and deceleration issues among OV model. Therefore they proposed the following GF (General Force) model:

$ \frac{{{\rm{d}}{v_n}}}{{{\rm{d}}t}} = {a_1}\left[ {V\left( {\Delta {x_n}} \right) - {v_n}} \right] + {\lambda _1}\Delta {v_n}H\left( { - \Delta {v_n}} \right) $

(2)

where, Δv_n=v_n+1-v_n is the relative velocity, H(-Δv_n) is the Heaviside function, and λ₁ is the reaction coefficient of the relative velocity to the vehicle in front.

In the real traffic flow, when the current car is much faster than the following one, the following car will not slow down although the space headway is shorter than the safe distance. Neither the OV model nor the GF model can explain this phenomenon, so Jiang et al.^[7] proposed the FVD (Full Velocity Difference) model:

$ \frac{{{\rm{d}}{v_n}}}{{{\rm{d}}t}} = {a_1}\left[ {V\left( {\Delta {x_n}} \right) - {v_n}} \right] + {\lambda _1}\Delta {v_n} $

(3)

All the above models are proposed based on a construction of normal road traffic conditions. However, the visibility and friction factor decline significantly under ice and snowfall conditions, where drivers will keep bigger headway and space headway as well as make slighter acceleration and deceleration. Therefore, considering the impact of ice and snowfall conditions on car-following behavior, this paper proposes the car-following model under ice and snowfall conditions based on Jiang's model^[7].

$ {{x''}_n}\left( t \right) = a\left[ {V\left( {\Delta {x_n}\left( t \right)} \right) - {v_n}\left( t \right)} \right] + {\mu _0}\left( {{f_r}/f_r^0} \right)\Delta {v_n} $

(4)

where, μ₀ represents the reaction coefficient of the relative velocity of the vehicle in front and the road surface friction coefficient under ice and snowfall conditions, it may be gotten by field observation; f_r⁰ is the road surface friction coefficient under normal traffic condition; f_r is the road surface friction coefficient under ice and snowfall conditions. f_r⁰ and f_r may be gotten by field observation, under very smooth ice film f_r=0.1, very smooth compacted snow f_r=0.15, ice sheet f_r=0.175, ice film f_r=0.225, ice sheet under snow f_r=0.25, mild compacted snow f_r=0.3 as well as f_r⁰=0.6.

3 Linear Stability Analysis

In this section, according to Refs.[7, 14], we study the linear stability of the model established in this paper. Setting λ=μ₀(f_r/f_r⁰), assuming that the initial state is steady. The space headway is b. The corresponding speed is V(b). Now the solution of vehicle location of the steady traffic flow is:

$ x_n^0\left( t \right) = bn + V\left( b \right)t $

(5)

For uniform flow solution (5), plus a perturbation y_n(t), we get:

$ \begin{array}{l} {x_n}\left( t \right) = x_n^0\left( t \right) + {y_n}\left( t \right)\\ {y_n}\left( t \right) = {x_n}\left( t \right) - x_n^0\left( t \right) \end{array} $

(6)

Seeking the first and the second derivatives, we get:

$\begin{array}{l} y{'_n}\left( t \right) = x{'_n}\left( t \right) - x{'_n}^0\left( t \right)\\ \;\;\;\;\;y'{'_n}\left( t \right) = x'{'_n}\left( t \right) \end{array}$

where $x{'_n}(t) = {v_n},x{'_n}^0(t) = v(b) = {\left. {\frac{{{\rm{d}}V(\Delta {x_n})}}{{{\rm{d}}\Delta {x_n}}}} \right|_{\Delta {x_n} = b}}$.

Then

$ \begin{array}{*{20}{c}} {\Delta {x_n} = {x_{n + 1}} - {x_n} = b + {y_{n + 1}} - {y_n} = b + \Delta {y_n}}\\ {\Delta {v_n} = {{x'}_{n + 1}} - {{x'}_n} = {{y'}_{n + 1}} - {{y'}_n} = \Delta {{y'}_n}} \end{array} $

Substituting the above equations into Eq.(1), linearizing the equation and getting:

$ \begin{array}{l} {{y''}_n} = a\left[ {V\left( {\Delta {x_n}} \right) - {v_n}} \right] + f\Delta {v_n} = a\left[ {V\left( {b + \Delta {y_n}} \right) - } \right.\\ \;\;\;\;\;\;\;\left. {v\left( b \right) - {{y'}_n}} \right] + f\Delta {v_n} \end{array} $

(7)

Setting y_n(t)=e^(zt+α_ni), by Fourier series expansion for y_n(t), and the results substituted into Eq.(7), and then we can get:

$ {z^2} + z\left[ {a - f\left( {{{\rm{e}}^{{a_k}{\rm{i}}}} - 1} \right)} \right] - aV'\left( b \right)\left( {{{\rm{e}}^{{a_k}{\rm{i}}}} - 1} \right) = 0 $

(8)

Setting z=ωi, e^a_ki=cos α_k+isin α_k, we get:

$ \begin{array}{l} - {\omega ^2} + wi\left[ {a - f\left( {\cos {\alpha _k} - 1 + {\rm{isin}}{\alpha _k}} \right)} \right] - \\ \;\;aV'\left( b \right)\left( {\cos {\alpha _k} - 1 + {\rm{isin}}{\alpha _k}} \right) = 0 \end{array} $

(9)

Making the real and the imaginary part zero, we get:

$ \begin{array}{l} V'\left( b \right) = {\left\{ {\left( {1 - \cos {\alpha _k}} \right)\left[ {a + f\left( {1 - \cos {\alpha _k}} \right)} \right]} \right.^2} + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;f{{\sin }^2}{\alpha _k}\left[ {a + f\left( {1 - \cos {\alpha _k}} \right)} \right]} \right\}/a{\sin ^2}{\alpha _k} \end{array} $

(10)

When a_k→0, apply the L'Hospital law to obtain the critical stable curve:

$ V'\left( b \right) = \lambda + \frac{a}{2} = {\mu _0}\left( {{f_r}/f_r^0} \right) + \frac{a}{2} $

When the space headway satisfies the following relationship, the system will be in a steady state

$ V'\left( b \right) < {\mu _0}\left( {{f_r}/f_r^0} \right) + \frac{a}{2} $

(11)

When f_r/f_r⁰=1, Eq.(11) is consistent with the model stability in Refs.[7, 14]. Meeting Eq.(11), the disturbance is small, the traffic flow will be stable; on the contrary, the traffic flow will be unstable. With the introduction of road surface friction coefficient in the model, the stability of traffic flow under different road conditions can be explored.

Parameter f_r is the road surface friction coefficient under ice and snowfall conditions. It may be gotten by field observation, under very smooth ice film f_r=0.1, very smooth compacted snow f_r=0.15, ice sheet f_r=0.175, ice film f_r=0.225, ice sheet under snow f_r=0.25 as well as mild compacted snow f_r=0.3. The simulation rendering critical stable curves under the seven different road conditions are shown in Fig. 1, where the region above the curve is stable, otherwise it is unstable. It shows that when road conditions become harsher, especially the friction coefficient becomes smaller, the stable area of traffic flow gets much smaller.

4 Nonlinear Analysis

In this section, with Refs.[7, 14-15], we study the nonlinear analysis of the model established in this paper. When traffic flow is in the unstable region or near the critical point, the slow variables will change. Density wave is the form of kink-anti-kink solitary wave^[15]. To refer to Ref.[15] for more details. When m=1, the model is equivalent to this paper's situation.

$ {{x''}_n}\left( t \right) = a\left[ {V\left( {\Delta {x_n}\left( t \right)} \right) - {v_n}\left( t \right)} \right] + \lambda \Delta {v_n} $

(12)

Using perturbation method kinetic for nonlinear analysis to obtain kink-anti-kink solitary wave solutions. Eq.(12) can be rewritten as:

$ \begin{array}{l} \frac{{{{\rm{d}}^2}\left( {\Delta {x_n}\left( t \right)} \right)}}{{{\rm{d}}t}} = a\left[ {V\left( {\Delta {x_{n + 1}}\left( t \right)} \right) - V\left( {\Delta {x_n}\left( t \right)} \right) - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{{{\rm{d}}\Delta {x_n}}}{{{\rm{d}}t}}} \right] + \lambda \left( {\frac{{{\rm{d}}\left( {\Delta {x_n} - 1} \right)}}{{{\rm{d}}t}} - \frac{{{\rm{d}}\left( {\Delta {x_n}} \right)}}{{{\rm{d}}t}}} \right) \end{array} $

(13)

Considering changes in the slow variable behavior near the critical point, a=a_c, $\varepsilon = \sqrt {\frac{{{a_{\rm{c}}}}}{{a-1}}} $, 0 < ε≤1.

And defining the variables X, T:

$ \begin{array}{l} X = \varepsilon \left( {N + bt} \right)\\ \;\;\;\;\;\;\;T = {\varepsilon ^3}t \end{array} $

where b is the unknown parameter.

Making space headway as Δx_n(t)=h_c+εR(X, T).

Eq.(13) expands to ε⁵ order, and then the nonlinear partial differential equation is obtained as follows.

$ \begin{array}{l} {\varepsilon ^2}\left( {b - V'} \right){\partial _X}R + {\varepsilon ^3}\left( {\frac{{{b^2}}}{a} - \frac{{V'}}{2} - \lambda \frac{b}{a}} \right)\partial _X^2R + \\ \;\;\;\;{\varepsilon ^4}\left[ {\partial _T^2R + \left( { - \frac{{V'}}{6} - \lambda \frac{b}{{2a}}} \right)\partial _X^3R - \frac{{V'''}}{6}\partial _X^3{R^3}} \right] + \\ \;\;\;\;{\varepsilon ^5}\left[ {\left( {\frac{{2b}}{a} - \frac{\lambda }{a}} \right){\partial _T}{\partial _X}R + \left( { - \frac{{V'}}{{24}} - \frac{{\lambda b}}{{24a}}} \right)} \right.\partial _X^4R - \\ \;\;\;\;\;\left. {\frac{{V'''}}{{12}}\partial _X^2{R^3}} \right] = 0 \end{array} $

(14)

where

$ V' = \frac{{{\rm{d}}\mathit{V}\left( {\Delta x} \right)}}{{{\rm{d}}\Delta x}}\left| {_{\Delta x = {h_c}}} \right.,V''' = \frac{{{{\rm{d}}^3}\mathit{V}\left( {\Delta x} \right)}}{{{\rm{d}}\Delta {x^3}}}\left| {_{\Delta x = {h_c}}} \right. $

Taking b=V′, $\frac{{{a_{\rm{c}}}}}{a} = 1 + {\varepsilon ^2}$, at the critical point near (a_c, h_c), Eq.(14) can be simplified as:

$ \begin{array}{l} {\varepsilon ^4}\left[ {{\partial _T}R + {g_1}V'\partial _X^3R + {g_2}{\partial _X}{R^3}} \right] + {\varepsilon ^5}\left[ {{g_3}\partial _X^2R + } \right.\\ \;\;\;\;\;\;\left. {{g_4}\partial _X^2{R^3} - {g_5}\partial _X^4R} \right] = 0 \end{array} $

(15)

where

$ \begin{array}{*{20}{c}} {{g_1} = \frac{{V'}}{6} + \frac{1}{2}\frac{{\lambda V'}}{{{a_c}}},{g_2} = - \frac{{V'''}}{6},{g_3} = \frac{{V'}}{2},}\\ {{g_4} = \left( {\frac{{2V'}}{{{a_c}}} - \frac{\lambda }{{{a_c}}}} \right)\left( {\frac{{V'}}{6} + \frac{1}{2}\frac{{\lambda V'}}{{{a_c}}}} \right) + \left( { - \frac{{V'}}{{24}} + \frac{1}{6}\frac{{\lambda V'}}{{{a_c}}}} \right)}\\ {{g_5} = - \frac{{V'}}{{12}}\left( {1 - \frac{{4V'}}{{{a_c}}} - 2\lambda } \right)} \end{array} $

Making the following transformation to Eq.(15), setting T′=g₁V′T, $R = \sqrt {\frac{{{g_1}}}{{{g_2}}}} R'$, we can get the MKdV equation which contains o(ε) correction term:

$ {\partial _T}R' - \partial _X^3R' + {\partial _X}{{R'}^3} + \varepsilon M\left[ {R'} \right] = 0 $

(16)

$ M\left[ {R'} \right] = \sqrt {\frac{1}{{{g_1}}}} \left[ {{g_3}\partial _X^2R' + \frac{{{g_1}{g_5}}}{{{g_2}}}\partial _X^2{{R'}^3} + {g_4}\partial _X^4R'} \right] $

(17)

Omitted o(ε) item, the kink-anti-kink wave solution is:

$ {{R'}_0}\left( {X,T'} \right) = \sqrt c \tanh \sqrt {c/2} \left( {X - cT'} \right) $

(18)

In order to obtain the propagation velocity c to Eq.(17), solvability condition must be met:

$ \left( {{{R'}_0},M\left[ {{{R'}_0}} \right]} \right) = \int_{ - \infty }^{ + \infty } {{{R'}_0}\left( {X,T'} \right)M\left[ {{{R'}_0}\left( {X,T'} \right)} \right]{\rm{d}}X} = 0 $

(19)

By integrating, we can get the speed of kink-anti-kink wave propagation

$ c = 5{g_2}{g_3}/\left( {2{g_2}{g_4} - 3{g_1}{g_5}} \right) $

so the solution to above equation is:

$ \begin{array}{l} R\left( {X,T} \right) = \sqrt {\left( { - V' + \frac{{3\lambda V'}}{{{a_c}}}} \right)\left( {\frac{{{a_c}}}{a} - 1} \right)\frac{c}{{V'''}}} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\tanh \sqrt {\frac{c}{2}} \left( {X - \frac{{cTV'}}{6} + \frac{{\lambda V'}}{{2{a_c}}}} \right) \end{array} $

(20)

so the solution to kink-anti-kink density wave of the space headway is:

$ \begin{array}{l} \Delta {x_n}\left( t \right) = {h_c} + \sqrt {\left( { - V' + \frac{{3\lambda V'}}{{{a_c}}}} \right)\left( {\frac{{{a_c}}}{a} - 1} \right)\frac{c}{{V'''}}} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\tanh \sqrt {\frac{c}{2}} \left( {X - \frac{{cTV'}}{6} + \frac{{\lambda V'}}{{2{a_c}}}} \right) \end{array} $

(21)

so the amplitude of the solution to kink-anti-kink density wave is:

$ B = \sqrt {\left( { - V' + \frac{{3\lambda V'}}{{{a_c}}}} \right)\left( {\frac{{{a_c}}}{a} - 1} \right)\frac{c}{{V'''}}} $

Kink wave solutions stand for the coexisting phases. And the coexistence curve can be plotted by the following formula:

$ \Delta {x_n}\left( t \right) = {h_c} \pm B $

(22)

5 Numerical Simulation

To verify the application of the model proposed in this paper, a numerical simulation was conducted. In the simulation, the OV function proposed by Helbing et al. was adopted^[8].

$ V\left( {\Delta x} \right) = {v_1} + {v_2} + \tanh \left[ {{C_1}\left( {\Delta x - {l_c}} \right) - {C_2}} \right] $

(23)

where, l_c is a car length, taking 5 m, v₁=6.75 m/s, v₂=71.91 m/s, C₁=0.13 m^-1, C₂=1.57 m^-1. In the periodic boundary conditions, take the loop length L=1 500 m, the number of vehicles N=100, the initial disturbance is:

$ \begin{array}{*{20}{c}} {{x_1}\left( 0 \right) = 10{\rm{m}},{x_n}\left( 0 \right) = \left( {n - 1} \right)L/N,\left( {n = 2,3, \cdots ,N} \right)}\\ {{v_n}\left( 0 \right) = V\left( {L/N} \right)\left( {n = 1,2,3, \cdots ,N} \right)} \end{array} $

Taking a=1.85 s^-1, μ₀=0.2 s^-1, from the critical stability curve we know all the models under various road conditions are stable. Where there is initial disturbance, taking t=100, 300, 1 000, 1 500 s, the numerical simulation of vehicle speed under different road conditions are shown in Fig. 2.

From the simulation results in Fig. 2, it is found that the speed of vehicles under different road conditions fluctuates around the initial velocity v₀=4.664 7 m/s with the initial disturbance for a long time. As time goes by, the speed of the vehicle will tend to be stable. At the same time, it is obvious that the harsher road conditions are, the larger vehicle speed volatility is, which indicates traffic flow more unstable. If the sensitivity coefficient takes a smaller value, the traffic flow under poor road conditions may lose its original stability, further causing traffic congestion. Therefore, the worse road condition is, the more sensitive to the model will be. This is consistent with previous results of linear stability analysis.

Taking a=1.5, t=1 500 s, the numerical simulation under different road conditions the phase space (space headway -the velocity space) hysteresis loop, as shown in Fig. 3. By the linear stability analysis, it shows, a=1.5 in this case, the road condition and traffic are not stable, thus in the process of evolution, vehicles will change from unimpeded to crowded, and then from crowded back to unimpeded. However, such a process is not the same as closing a reply. There is hysteresis. From Fig. 3, the worse road condition is, the larger the area formed by the phase space in the closed hysteresis loop will be (i.e. the greater the instability region). When the value of a is taken to meet the stability condition, only one dot exists in the phase space not forming a hysteresis loop. This is consistent with the previous results of nonlinear analysis.

Numerical simulation under different road conditions at the same time (t=1 500 s) and different time (t=1 490 s and t=1 500 s) the same road conditions (in a very smooth ice film for example) all vehicles space headway distribution, as shown in Figs. 4, 5 and 6, respectively. Through the non-linear analysis, traffic flow in the unstable region, the critical point (h_c, a_c) density wave is kink-anti-kink solitary wave form. Figs. 4 and 5 illustrate that the amplitude of density waves increases along with the deterioration of road conditions, that is to say, traffic flow easily causes instability in poor road conditions, forming stop-go density waves. Fig. 6 shows that the waveforms in the unstable region at different moments are exactly the same. When density waves propagate upstream, the vehicles in front will decelerate and park, which may lead to queuing and congestion.

6 Conclusions

This paper proposes an improved car-following model considering the coefficient of road friction, based on the influence of different road conditions on propagation characteristics of traffic wave under ice and snowfall conditions. After the linear stability analysis, stability conditions of car following model under different road conditions are identified. In addition, critical stable curves under seven different road conditions are drawn. The results show that when road conditions become harsher, especially the friction coefficient diminishes, the stable area of traffic flow gets much smaller. Traffic congestion is the limiting form of density wave propagation caused by the interaction of the vehicles. It can be concluded from nonlinear analysis that density wave is in kink-anti-kink solitary wave form, considering variation behavior of the slow variables in the unstable region or near the critical point of traffic flow. Through deduction, kink-anti-kink density wave solutions of distance headway and coexistence curve are obtained. Numerical simulation results show that worse road conditions or greater speed fluctuations will lead to greater instability region of traffic flow and less absorption interference ability, even it can cause more congestion and accidents. Ultimately, the conclusions have great significance in solving the problems of traffic congestions and security under ice and snowfall conditions.