Numerical Simulation

Chapter 4 Numerical Simulation We present the results of the traffic flow simulation for various traffic flow parameters. We develop a computer program for the explicit upwind difference scheme and implement the scheme for artificially generated initial and boundary data and verify the well-known qualitative behaviors of different flow variables. We also estimate the error of the numerical scheme. Finally we try to apply our model for the traffic flow analysis in a 10 (ten) km range of Dhaka-Aricha highway.

4.1 Incorporation of initial & boundary data To obtain the density profile we have used the explicit upwind difference scheme. In order to use the scheme we have made the following assumptions:

 We have considered a highway in a range of 10 (ten) Km.  We consider the number of vehicles of various points at a particular time as initial data and constant boundary ρ( 0, t ) = 21 / 0.1 km . .

 We have estimated maximum density which is parameterized by ρmax in the traffic flow model. To evaluate ρ max , we have used the following equation

ρ max = max ρ ( 0, x ) * c For this we use the initial density for ρ( 0, x ) and take c = 10 as a constant. Using initial and boundary value on explicit upwind difference scheme, we can forecast the traffic flow model.

Numerical Simulation

34

4.2 Algorithm for the numerical solution of traffic flow model To find the numerical solution of the IBVP (2.5) we have to store some variables and values which are presented in the following algorithm.

Input: nx

and

nt the number of spatial and temporal grid points respectively.

tf , the right end point of

xd , the right end point of

[0, T ] . [0, b] .

ρ0 , the initial density, apply as a initial condition.

ρa , the vehicles at starting point, apply as a boundary value. V max , maximum velocity.

Output: ρ( x, t ) Initialization:

the solution matrix. dt =

T −0 , the temporal grid size nt

dx =

b −0 , the spatial grid size. nx

gm = V max*

dt , the courant number. dx

gma =1 − gm

c = 10 , a constant. (c ≥ 3 )

ρ max = max ( ρ 0) * c , maximum density.

Step 1: Calculation for upwind scheme For i =1 to nt For j = 2 to nx + 1 ρ (i + 1, j ) = gma * ρ (i, j ) + gm * ρ (i, j − 1) + (( dt / dx ) * (V max/ ρ max^ 2)) * ( ρ (i, j ).^ 3 − ρ (i, j − 1).^ 3)

End End

Step 2: Calculation for velocity For i =1 to nt Numerical Simulation

35

v( i ) = V max* (1 − ρ (i ).^ 2 / ρ max .^ 2)

End

Step 3: Calculation for flux For i =1 to nt q( i ) = V max* (u (i ) − ρ (i ).^ 3 / ρ max .^ 2)

End

Step 4: Output ρ( x, t ) Step 5: Graph Presentation Step 6: Stop

4.3 Error estimation of the numerical scheme

Numerical Simulation

36

Since the analytic solution of the traffic model is in implicit form, as a result the error is much more difficult to compute. However, we try to estimate the error for the explicit upwind difference scheme for the linear advection equation which is the simplest form of traffic model. We consider a linear advection equation ∂u 1 ∂u  1 + = 0 ; x ∈( 0,4π ), t ∈ 0,  ∂t 2 ∂x  2

− − − −( 4.1)

with initial and boundary conditions u ( 0, x ) = sin x u (t ,0) = t 2

− − − −( 4.2)

Equations (4.1) and (4.2) give an initial boundary value problem (IBVP). The analytic solution of equation (4.1) by using method of characteristics is t  u ( t , x ) = Sin  x −  2 

− − − −( 4.3)

The numerical solution of the equation (4.1) is calculated by using explicit upwind difference scheme. The computer program performed the percentage error as err = 1.466% of the linear advection equation (4.1). Fig: 4.1 represent the analytic solution and numerical solution of the equation (4.1) simultaneously. Analytic Solution

Numerical Solution

Fig 4.1: Error of the numerical solution Fig 4.2 represents the error profile with respect to time for the explicit upwind difference scheme for the linear advection equation.

Numerical Simulation

37

Error in Percentage

1.5

1

0.5

0

0

2

4

6 8 Time [Sec]

10

12

14

Fig 4.2: Error profile This error can be reduced for further smaller discretization parameters ∆t , ∆x and this justifies the correctness of the computer program. Therefore we now apply this program for explicit upwind difference scheme of our model considered for traffic simulation.

4.4 Numerical results for traffic flow Numerical Simulation

38

According to focusing on some parameters we present some specific cases of traffic flow which are as follows:

Case A: In this case we choose maximum velocity, Vmax = 60 km / hour . It is notified that for satisfying the CFL condition we pick the unit of velocity as km / sec . We consider ρmax = 250 / km and perform numerical experiment for 6 minutes in ∆t = 90 time steps (temporal grid size) for a highway of 10 km in 51 spatial grid points with step sizes ∆x = 100 meters. We consider the initial density ρ( 0, x ) as presented in Fig: 4.3 and the constant boundary value ρ( t ,0 ) = 21 / 0.1 km .

24 2 3 .5

Density of Car /0.1km

23 2 2 .5 22 2 1 .5 21 2 0 .5 20

0

1

2

3

4 5 6 H ig h w a y [k m ]

7

8

9

10

Fig: 4.3 Initial traffic density in 0.1 km

We run the program and attain the density profiles as presented in Fig: 4.3, the computation time is calculated as 6.1892 seconds.

Numerical Simulation

39

Fig: 4.4 show the initial position of car as well as the position of car after 6 minutes with respect to the certain points of 10 km highway.

In it ia l 6 m in u t e s

24 2 3 .5

Density of Car /0.1km

23 2 2 .5 22 2 1 .5 21 2 0 .5 20 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig 4.4: Comparison between initial position & after 6 minutes

In fig: 4.5, the curve marked by solid line represents the density of car at 2 minutes, and the curve manifested as “-o-“ represents the density profile at 4 minutes also the curve manifested as “-*-“ represents the density profile at 6 minutes correspondingly.

Numerical Simulation

40

2 2 .6

2 m in 4 m in 6 m in

2 2 .4 2 2 .2

Density of Car /0.1km

22 2 1 .8 2 1 .6 2 1 .4 2 1 .2 21 2 0 .8 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig: 4.5 Density of Car for 6 minutes

Fig: 4.6 represent the density profile of car of 10 km highway for 9 minutes. The curve marked by solid line represents the density of car at 3 minutes, and the curve marked as ‘-o-‘ represents the density of car at 6 minutes, also the curve marked as ‘-*-‘ represents the density of car at 9 minutes respectively.

Numerical Simulation

41

2 2 .6

3 m in 6 m in 9 m in

2 2 .4 2 2 .2

Density of Car /0.1km

22 2 1 .8 2 1 .6 2 1 .4 2 1 .2 21 2 0 .8 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig: 4.6 Density of Car for 9 minutes

Fig: 4.7 represent the density profile of car of 10 km highway for 12 minutes. The curve marked by solid line represents the density of car at 4 minutes, and the curve marked as ‘-o-‘ represents the density of car at 8 minutes, also the curve marked as ‘-*-‘

Numerical Simulation

42

represents the density of car at 12 minutes respectively. We observe that there is no obstacle in the next location of the highway.

2 2 .6 4 m in 8 m in 1 2 m in

2 2 .4 2 2 .2

Density of Car /0.1km

22 2 1 .8 2 1 .6 2 1 .4 2 1 .2 21 2 0 .8 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig: 4.7 Car Density for 12 minutes

Fig: 4.8 represent the respective computed velocity profile according to the certain points of the highway. The velocity is computed by the following relation

 ρ2 V ( ρ ) = V max 1 − 2 ρ max 

 .  

Numerical Simulation

43

0.16 2 m in 4 m in 6 m in

0 .1 5 9 8

V elocity of Car [0.1km /Sec]

0 .1 5 9 6 0 .1 5 9 4 0 .1 5 9 2 0.159 0 .1 5 8 8 0 .1 5 8 6 0 .1 5 8 4 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig 4.8: Velocity Profile

Now, we are known about the density and velocity for certain points. So we can calculate the flux with the aid of the relation, q = ρ * V . Fig 4.9 represents the computed flux with respect to the distance.

Numerical Simulation

44

3.75 2 m in 4 m in 6 m in

3 .7

Flux

3.65

3 .6

3.55

3 .5

3.45 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig: 4.9 Flux of traffic

In the fig: 4.10, the computed velocity is plotted with respect to the density profile for the given traffic data which gives a nonlinear curve as desired.

Numerical Simulation

45

V e lo c i t y a s a fu n c t io n o f d e n s it y ( N o n - li n e a r ) 0 .1 8 0 .1 6 0 .1 4

Velocity [0.1 km /sec]

0 .1 2 0 .1 0 .0 8 0 .0 6 0 .0 4 0 .0 2 0

0

5

10 15 D e n s it y / 0 . 1 k m

20

25

Fig: 4.10 Traffic velocity as a function of density

In the fig: 4.11, the computed flux is also plotted with respect to the density profile for the given traffic data and the fundamental diagram is verified showing traffic flow as a concave shape in the range 0 ≤ ρ ≤ ρmax .

Numerical Simulation

46

F u n d a m e n t a l D ia g r a m 0.1 8 0.1 6 0.1 4 0.1 2

Flow

0 .1 0.0 8 0.0 6 0.0 4 0.0 2 0

0

5

10 15 D e n s it y / 0 . 1 k m

Fig: 4.11 Flux as a function of density

Case B:

20

25

0 ≤ ρ ≤ ρ max

In this case, we reduce the parameter of maximum velocity

V max = 30 km / hour but treating the same ρmax = 250 / km with the same initial density

as in the case A. As maximum velocity is reduced by a factor 2(two) from the previous case, the density is also increased that is the speed is decreased. The computed

Numerical Simulation

47

density profile as shown in fig: 4.12 and desired traffic waves are moving much slower than that in the previous case.

23 2 m in 4 m in 6 m in

Density of Car /0.1km

22.5

22

21.5

21

20.5 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

Fig: 4.12 Car Density waves moves slower for smaller

Case C:

10

V max

If we take a source term after 5 km of our considered 10 km highway by

which some vehicles are entered but suppose the same maximum velocity and density also the initial density as in case A. The nature of the computed density profile is presented in fig: 4.13. We see that after 5 km, density is increased and after 6 minutes it is backed to the previous situation likely. Numerical Simulation

48

24

2 m in 4 m in 6 m in

23 .5

Density of Car /0.1km

23

22 .5

22

21 .5

21

20 .5 0

1

2

3

4 5 6 H ig h w a y [ k m ]

7

8

9

10

Fig 4.13: Density profile for a source term

4.5 Implementation of the traffic flow model To implement the traffic flow model for a real situation, we need some traffic data of a particular highway. For this purpose, we collected some data from the Roads and

Numerical Simulation

49

Highways Department (RHD), Manikgonj division in a range of 10 (ten) km of DhakaAricha highway. Those data are in tabular form as follows:

Table-1: Number of vehicles Time 6.00-

N.Nagar

Nayarhat

Dhamrai

Kalampur

7.00 AM

218

205

201

197

Dulivita

Arepara

180

171

Table-2: Time

Number of vehicles at Nabinagar.

6.00 AM

210

7.00 AM

217

8.00 AM

221

9.00 AM

236

Source: Roads and Highways Department (RHD), Manikgonj division. The data of table-1 and table-2 gives the flux not densities at that point. So we have to generate the density and velocity. We choose the density by approximation based on the given flux as follows: Time 6.00 AM

Density / km

142

144

153

160

167

We compute the velocity using the relation, V = Time 6.00 AM

172

q

ρ

180

189

205

209

0.878

0.818

as given below.

Velocity km/Sec 1.422

1.424

1.392

1.325

1.305

1.25

1.117

1.042

We want to test the data for nonlinear density-velocity relationship. That is why we plot the velocity according to the density.

Numerical Simulation

50

Fig 4.14: Density-velocity relation of Dhaka-Aricha highway From the above figure, Fig:4.14 we say that the curve can be approximated by a nonlinear density –velocity function because those two curves are very close to each other and therefore this data can be used for our considered model with concave shape relation.

Numerical Simulation

51