A model for a network of conveyor belts with discontinuous speed and capacity

1.   Introduction

Conveyor belts have attained a favored position in transporting bulk materials due to their economy, reliability, safety, versatility and almost unlimited range of variations, conveying a wide variety of materials. However, theoretical study of these systems is non-trivial since the presence of intersections, fixed or smart divert/merge devices, and different production speeds add complexity to the structure and generate interesting dynamics. For a general presentation of the models in the literature, we refer to the monograph [3] and the references therein. As more recent contributions, we mention in particular the works [1,8,10,11,13], where the authors introduce a production (or traffic) model with discontinuous flux and study the properties of solutions. While in [1,10,13], the focus is on the mathematical modeling, and numerical simulation of the discontinuous flux function, the authors in [8,11] prove the existence of solutions using wave-front tracking.

In this work, we consider a production network consisting of multiple linked conveyor belts. We remark that the new network model differs essentially from [10] in the proposed coupling conditions and the concept of solutions. Each arc in the network corresponds to a conveyor belt with a certain constant speed and capacity along the arc. We suppose that the capacity is limited meaning that a maximum value of the density cannot be exceeded. Since we assume no buffers at the intersections, in the case of capacity drop the products get stuck on the belt (but they are still transported with a certain velocity). This is a key difference to traffic flow models [7,9], where the velocity is dependent on the density and vehicles have to stop once a maximal capacity is reached.

The original contribution of the present paper consists in describing a class of entropy solution of special relevance for this model taking into account that uniqueness in general is not guaranteed. Moreover, we show how these solutions can be successfully approximated by appropriate Gordunov schemes.

From an application viewpoint, the network model is particularly appropriate to study production systems for bottling, canning, and packaging. Figure 1 shows a brewery, where beer bottles are transported through a system of conveyor belts. We observe different lanes converging in an accumulator device and some ample space in the middle to stock the units that cannot be absorbed by the system.

The paper is organized as follows: in Section 2 we discuss the basic model and the notion of weak solutions which permits to derive the existence of a solution. We extend the framework to networks in Section 3: for one-to-one junction (Sec. 3.1) as well as diverging and merging intersection (Sect. 3.3 and 3.2) where we build an analytic solution for the model. In Section 4, we introduce a suitable discretization method to tackle the network model and present in Section 5 numerical results for different network settings where we show how the numerical approximations converge to the (previously described) analytic solution of interest.

2.   Basic model

We start recalling the model originally introduced in [1] and generalized to networks in [10]. In both articles, a production system is described by a conservation law with discontinuous flux function and, if different from zero, constant speed $ a>0 $. More precisely, setting the model domain as $ {\mathbb{R}} $ and calling $ \rho: {\mathbb{R}}\times[0, T]\rightarrow [0, \rho^{\text{max}}] $ the product density, the evolution of the system is

where $ H $ is the Heaviside function and the initial data $ \rho^0 $ is a function of bounded variation satisfying $ \rho^0(x)\leq \rho^{\text{max}} $.

We point out that on a single arc (assuming $ \rho^0(x)< \rho^{\text{max}} $ for all $ x\in {\mathbb{R}} $) the equation simply reduces to an advection equation and therefore the solution is simply

This clearly differs from traffic models where the speed typically depends on the local density and shocks can appear on a single arc even with smooth initial data [7,9]. The presence of discontinuities in the flux function, anyway, poses some problems in the case of $ \rho^0(x) = \rho^{\text{max}} $ for some $ x\in {\mathbb{R}} $.

A classic approach to deal with this problem is to use a standard regularization of the flux function using some Friedrichs' mollifiers. As it has been shown in [5], the solutions of the regularized problem converge to a bounded entropy weak solution, in the particular sense of the definitions below. Therefore, we denote by $ \tilde f $ the following multivalued function

Definition 2.1. A function $ \rho\in L^\infty( {\mathbb{R}}\times [0, T]) $ is called weak solution to the Cauchy problem (1) if there exists a function $ v\in L^{\infty}( {\mathbb{R}}\times [0, T]) $ such that $ v(x, t)\in\tilde f(\rho) $ a.e. and

for each $ \phi\in C^1_c( {\mathbb{R}}\times [0, T]) $ (where $ \phi\in C^1_c $ means $ \phi\in C^1 $ with compact support).

Classically, the definition above is completed by the following notion of entropy weak solutions: denote by $ \tilde H $ the following multivalued function

Definition 2.2. A weak solution $ \rho $ of the Cauchy problem (1) is called an entropy weak solution if, for each entropy $ \eta\in C^1( {\mathbb{R}}) $, $ \eta $ convex, there exists a function $ w\in L^{\infty}( {\mathbb{R}}\times [0, T]) $ such that $ w(x, t)\in\tilde H(\rho(x, t)) $ a.e. and

where

Remark. We observe that that the solution (2) is a weak entropy solution in the sense of Definitions 2.1 and 2.2. This can be shown by choosing

and

for any constant $ k\in {\mathbb{R}} $.

At the same time, as it is possible to see adapting the example provided in [4], in case of congested initial solution, a collection of weak solutions are acceptable as long as a condition (derived by the jump on the flux) on the speed of the congested area is verified. This means that in our case, for equation (1) with initial condition equal to

there are at least two distinct entropy solutions equal to

Since we are interested in the more meaningful solution from the point of view of the application, it seems clear that the advection formula provides a good candidate in the case of the presence of a congested region (differently from other applications e.g. traffic models). This point has been discussed in detail in [11], where the authors introduce the additional concept of congested/non congested region in order to select such unique solution. Instead, in the present paper, we focus on building an entropy solution in the more general case of networks and showing how a suitable numerical scheme provides a good approximation of it.

3.   Extension to networks

We now extend the model to a network represented by a directed graph $ \Gamma = (V, E) $ with $ V \neq \emptyset $ the set of vertices and $ E: = \{\Omega_i\} $, $ i\in \{1, ..., M\} $ the set of arcs. For a fixed node $ v \in V $, the sets $ \delta_v^-, \delta_v^+ $ denote the ingoing and outgoing arcs, respectively. We consider the following network problem:

where

is the flux function on arc $ i \in E $ and $ a_{{i}} \in {\mathbb{R}}^+ $ the transport velocity. At intersections, we assume the conservation of flux and some transition conditions which are discussed later in this section. Our main goal is to build an appropriate entropy solution. We underline that, differently from [10], speed and capacity are different for each arc leading to new discontinuities at the intersection points and some non-local phenomena.

The model (3) is difficult to handle in its generality. In the following we discuss separately various cases of junction networks, where only a vertex is considered and the arcs are $ M $ half lines, where each line is isometric to $ [0, +\infty) $.

3.1.   One-to-one junction

The first case we consider is a one-to-one junction. It can be seen as a special case of a one-dimensional problem, as described by (1), with discontinuities in the velocity $ a $ and the capacity $ \rho^{\text{max}} $.

For simplicity in the notation, we consider the problem on $ \Omega = \Omega_1\cup \{0\}\cup \Omega_2 = (-\infty, 0)\cup \{0\}\cup(0, \infty), $ where the intersection is located at $ x = 0 $. The model equations are given by (3), where the junction condition is specified as

with flux function (4) on arc $ i = 1, 2 $. The solution can be easily computed if no congestion occurs during the transportation, i.e.

In this case, we can derive the solution as follows. For $ x \in \Omega_1 $, the characteristics of the problem are simply the straight lines $ y(t) = x-a_{{1}} t $, which lead to the solution

Analogously, we find

which corresponds to the density initially placed on $ \Omega_2 $. The solution for $ (x, t) $ in the case $ 0<x<a_{{2}} t $ is found by the junction condition (5): $ f_1(\rho(x^+, t)) = f_2(\rho(x^-, t)) $ where $ x^\pm $ denotes $ \lim_{h \rightarrow 0} x \pm h $. This leads to

at the interface $ x = 0 $, which gives the solution of the intermediate part with adapted transport velocity. The solution is then described by

The characteristics of this solution are shown in Figure 2 for $ a_{{1}} < a_{{2}} $. Note that the solution might be discontinuous at $ x = 0 $ and $ x = a_{{2}}t $. If the initial data $ \rho^0_1 $ and $ \rho^0_2 $ is continuous, the solution $ \rho(x, t) $ keeps this property in all other points. Along each characteristic $ c \in {\mathbb{R}} $ is constant.

Remark. The interpretation of condition (6) is as follows: We know that the flux is $ a_{{1}} \rho(x, t) $ if the maximal density $ \rho^{\text{max}}_2 $ is not reached. In this case, the capacity of arc 1 before the intersection has no influence. This is due to the positive velocities $ a_{{1}} $ and $ a_{{2}} $. As we have already discovered, congested areas never appear in the interior of an arc, but only in conjunction with intersections.

If condition (6) is not satisfied, congestion arises and the problem becomes more involved. In particular, it is non-trivial to obtain a correct weak entropy solution using only Definitions 2.1 and 2.2.

Let $ t_0 $ denote the first point, where condition (6) is violated, i.e. the first time of congestion:

We track the interface describing the congested area at maximal density $ \rho^{\text{max}}_1 $ that can appear in $ \Omega_1 $ and call the congested region $ \Lambda $, see Figure 3. The interface is a time-dependent function $ g(t) $. The evolution of $ g(t) $, starting at time $ t_0 $, can be derived by integrating the difference between the fluxes entering and exiting the region $ \Lambda $ as well as the current density. The entering flux at time $ t $ is given by $ a_{{1}} \rho^0_1(-a_{{1}}t) $, the exiting flux by $ a_{{2}} \rho^{\text{max}}_2 $ since (6) is violated and the maximal density on the outgoing arc is reached. The resulting density in the congested region is $ \rho^0_1(y) $ for $ g(t)-a_{{1}}t \leq y \leq -a_{{1}}t $. Summarizing, this leads to

Rearranging the terms, we can describe the congested region $ \Lambda $ as

with interface $ g $ defined as $ g:[0, \infty) \rightarrow {\mathbb{R}}_{\leq 0} $, where $ t $ is mapped to the solution of

if this is negative, and to zero otherwise. The final time of congestion is $ t^E = \min\left\{t\geq t_0 \hbox{ such that }g(t) = 0\right\}. $ Notice that $ g'(t) $ is not necessarly monotone, thus the congestion region can decrease and than grow again without disappearing.

Figure 3 shows the trajectories of the problem in the case of congestion for $ a_{{1}} > a_{{2}} $. Outside the region $ \Lambda $, they correspond to the characteristics. If the initial data $ \rho^0_1 $ and $ \rho^0_2 $ are smooth, the solution may still become discontinuous in the interface $ x = g(t) $, in $ x = 0 $ and in $ x = a_{{2}}t $. The congested region $ \Lambda $ is highlighted in gray.

It is straightforward to see that the whole congested region can be constituted by multiple non-connected sets, if the congestion disappears and condition (6) is violated again. In that case, we set $ t_0^E : = t^E $ and there exists a $ t_k $ such that

The procedure to build this second connected subset of the whole congested region is the same as for the first one and so on. Therefore, it is sufficient to only consider one connected set $ \Lambda $.

Inside the region $ \Lambda $, the transport velocity $ \bar{a} $ is such that the coupling condition (5) holds true, i.e. the inflow $ \bar{a} \rho^{\text{max}}_1 $ equals the outflow $ a_{{2}} \rho^{\text{max}}_2 $ at $ x = 0 $. Therefore, the velocity inside the region $ \Lambda $ is

Remark. Since we assumed that condition (6) is violated, we find a value $ \bar x\in\Omega_1 $ such that $ a_{{1}} \rho^0_1(\bar x)>a_{{2}} \rho^{\text{max}}_2 $ and $ \rho^0_1(\bar{x})\leq \rho^{\text{max}}_1, $ where the second term is due to the choice $ \rho^0(x)\leq\rho^{\max} $. This implies

Dividing the first and the last term of this inequality by $ \rho^{\text{max}}_1 $, we obtain $ \bar a<a_{{1}} $, which means that the velocity always decreases as soon as the mass enters the maximal density area $ \Lambda $. This confirms the intuitive assumption that the transport velocity is reduced in the congested region $ \Lambda $.

Next, we derive a general solution for this case (shown in Figure 3) for the modifications in the characteristics.

We follow an approach considering the associated Riemann problem. We calculate the solution for a general initial condition, approximated by piecewise constant functions. The same approach is discussed in [11] for more general problems.

We briefly discuss the Riemann problem on $ [0, T)\times {\mathbb{R}} $, consisting of equation (3) and the initial data

If we allow waves with negative velocity, we can solve this Riemann problem by distinguishing the following cases:

a) $ f_1 \leq f^{\text{max}}_2 $: no congestion arises, (6) is verified and the solution is (7) with piecewise constant initial data $ \rho^0_1(x) = \rho_l $ and $ \rho^0_2(x) = \rho_r $.

b) $ f_1 > f^{\text{max}}_2 $: congestion arises since (6) is not verified. In this case, the solution consists of three shock waves (see Figure 4) starting at $ x = 0 $:

ⅰ) At $ x = 0 $, a shock with velocity $ s = 0 $ arises, where the solution $ \rho $ jumps from $ \rho^{\text{max}}_1 $ to $ \rho^{\text{max}}_2 $. For the special case $ \rho^{\text{max}}_1 = \rho^{\text{max}}_2 $, there is no jump in the solution at $ x = 0 $.

ⅱ) We obtain a left-going shock wave, where the density jumps from $ \rho_l $ to $ \rho^{\text{max}}_1 $ with negative velocity $ s_l < 0 $ (computed according to the Rankine-Hugoniot condition):

This shock wave describes the left boundary of the congested area $ \Lambda $, meaning that it follows the function $ g(t) $. Therefore, the transport right of this shock is with velocity $ \bar{a} $ as defined in (12). In the case $ \rho_l \rightarrow \rho^{\text{max}}_1 $, the shock velocity tends to $ - \infty $.

ⅲ) We obtain a right-going shock wave, where the density jumps from $ \rho^{\text{max}}_2 $ to $ \rho_r $ with positive velocity $ s_r = a_{{2}} > 0 $ (computed according to the Rankine-Hugoniot condition):

In the case $ t_0 > 0 $, we are in the non-congested case for $ t < t_0 $. At $ t = t_0 $, congestion starts and the shock waves appear.

In the general case with non-constant initial data, we obtain the following solution:

If $ g(t) = 0 $ for all $ t \in [0, T] $, we have $ \Lambda = \emptyset $ and recover the previously described solution (7).

Contrary to traffic flow models, rarefaction waves do not appear in the conveyor belt problem. This is due to the linear flux function up to the discontinuity and drops to zero.

3.2.   One-to-two junction

Now, we consider the case of a splitting intersection, where one arc is separated into two. We denote by $ i = 1 $ the incoming and by $ i = 2, 3 $ the outgoing arcs (see Figure 5).

The choice of a distribution rule depends on the application: if the intersection is "passive", we impose a fixed rate between the outgoing fluxes which is kept constant during the transportation process via a device $ D $ (see the left picture in Figure 5). In this situation, congestion arises even if the outgoing arcs are not both congested. A different approach is an "active" junction (see right picture in Figure 5). In this case, the ratio between the outgoing fluxes can change at the intersection. The aim is the maximization of the total outgoing flux and the reduction of congestion. We underline that this behavior partially differs from the standard coupling conditions introduced for vehicular traffic fluxes (cf. [7]), since the choice of the vehicles (left/right at a junction) cannot be determined by the local status of the traffic.

We consider the problem on $ \Omega = \Omega_1\, e_1 \cup \{0\} \cup \Omega_2\, e_1\cup\Omega_3\, e_2 = (-\infty, 0)\, e_1\cup\{0\} \cup(0, \infty)\, e_1\cup(0, \infty)\, e_2, $ where $ (e_1, e_2) $ is the standard base of $ {\mathbb{R}}^2 $, and we identify the element $ x \in \Omega_1 $ with the vector $ (x, 0)^T $, and analogously for arcs $ i = 2, 3 $. The equations that we consider are (3), where the junction condition is specified as

equipped with flux function (4) on arc $ i = 1, 2, 3 $.

Case 1: "passive" junction. At first, we consider the case of same flux rates between the two exiting arcs. An intersection device $ D $ keeps the ratio of the two outgoing fluxes constant, i.e., for a fixed distribution parameter $ \mu\in[0, 1] $, it holds

even if only one outgoing arc is congested. The case $ \mu \in \{0, 1\} $ reduces the problem to the one-to-one junction, so we consider $ \mu \in (0, 1) $ now. The relation (16) states that a fixed rate

is kept between the two outgoing fluxes during the evolution of the system, independent of the incoming flux. If no congestion arises, i.e.,

holds true, the solution is obtained in a similar way as in Section 3.1. We skip this point and draw our attention directly to a general formula for the solution (with or without congested areas).

Due to the constant rate (17) between the two outgoing fluxes $ f_2 $ and $ f_3 $, it might happen that only one outgoing arc reaches the maximal density before congestion on the incoming arc $ 1 $ arises. Therefore, we define the actual density on arc $ 2 $ and $ 3 $ as

This is the minimum of the density the arc is supposed to take due to the distribution parameter, and the maximal density possible on this arc. The first time of congestion $ t_0 $ can be determined as

We can introduce the interface $ g(t) $, in analogy to (11), with exiting flux $ \bar{\rho}_2 a_{{2}} + \bar{\rho}_3 a_{{3}} $. The interface is defined as $ g:[0, \infty) \rightarrow {\mathbb{R}}_{\leq 0} $, where $ t $ is mapped to the solution of equation

if this is negative, and to zero otherwise. The region of congestion $ \Lambda $ on $ \overline{\Omega}_1 $ is given by

analogously to (10). The general solution $ \rho $ on $ \Omega $ is given by

where $ \alpha_2 = \mu \hbox{ and } \alpha_3 = 1-\mu $ for a compact notation.

Remark. Congestion occurs if the maximal density of one of the two exiting arcs $ i = 2, 3 $ is reached. The other arc, even if the maximal density is not reached, shows a similar congested behavior, i.e., a value less than $ \rho^{\text{max}}_i $ is reached and kept. Since congestion arises without using the full capacity of both outgoing arcs, the duration of the congested phase is prolonged.

Case 2: "active" junction. We consider the possibility of a diverter interpreted as a device that keeps a constant ratio among the two outgoing fluxes as long as no arc is congested. If congestion arises, the device adapts the fluxes to ensure the maximal total outgoing flux, i.e., $ f_2 = f^{\text{max}}_2 $ and $ f_3 = f^{\text{max}}_3 $.

As before, we fix a parameter $ \mu \in [0, 1] $ in order to set a constant ratio in the non-congested case. As for the passive junction case (16), we have $ \mu f_1 = f_2 \hbox{ and } (1-\mu)f_1 = f_3. $ In order to get a unique solution also in the congested case, we define parameters $ \beta_i $ corresponding to $ \alpha_i $ in (21) with the following properties:

a) The flux conservation $ f_1 = f_2+f_3 $ at the coupling is satisfied.

b) The parameters $ \beta_i $ are equal to $ \alpha_i $, $ i = 2, 3 $ as in the previous case if no congestion occurs, i.e., condition (18) holds true.

c) If only one outgoing arc $ i $ is congested, i.e., $ \rho^0_1(x-a_{{1}} t)> \frac{a_{{i}}}{\alpha_i a_{{1}}} \rho^{\text{max}}_i $ at time $ t $ the parameter $ \beta_i $ changes to

which is the smallest value to avoid congestion.

d) A further change is necessary if the value

is reached. This is the optimal ratio to maximize the flux through the junction. At this point the congestion starts.

Having stated the conditions on $ \beta_i, \; i = 2, 3 $, we can now derive the solution. No congestion arises as long as the following inequality holds true:

If condition (22) is not fulfilled, the time $ t_0 $, meaning the first time of congestion, is independent of the distribution parameter $ \mu $. It is then given by

In this case, the interface $ g $ is defined as $ g:[0, \infty) \rightarrow {\mathbb{R}}_{\leq 0} $, where $ t $ is mapped to the solution of equation

if this solution is negative, and to zero otherwise. The region of congestion $ \Lambda $ on $ \Omega_1 $ is given by (20). Contrary to the passive junction case, the function $ g(t) $ considers the maximal flux on arcs $ 2 $ and $ 3 $ but no longer the distribution parameter $ \mu $.

We define for $ i = 2, 3 $

with $ \alpha_2 = \mu $ and $ \alpha_3 = 1-\mu $. The general solution on $ \Omega $ is then

with $ i = 2, 3 $.

Remark. Compared to the "passive" junction (21), congestion only occurs if the maximal capacity of both exiting arcs is reached. This implies that the choice of $ \mu \in \{0, 1\} $ does no longer reduce to a one-to-one junction.

3.3.   Two-to-one junction

In this part, we focus on the case of a merging junction. We know from traffic flow that in the free flow regime no additional information is needed. Conversely, in the congested case, we need a priority rule between the two incoming arcs, i.e., how to use released capacities of the outgoing arc. We denote by $ i = 1, 2 $ the incoming and by $ i = 3 $ the outgoing arcs.

We consider the problem on $ \Omega = \Omega_1\, e_1\cup \Omega_2\, e_2\cup\{0\}\cup\Omega_3\, e_1 = (-\infty, 0)\, e_1\cup (-\infty, 0)\, e_2\cup\{0\}\cup(0, \infty)\, e_1 $ with the same interpretation as before, where the system is given by (3), the coupling condition reads as

and the flux function is again (4) on arcs $ i = 1, 2, 3 $. The solution can be directly computed if it holds

i.e. no congestion arises. Then, the solution is given by

and the incoming mass can be totally absorbed by the outgoing arc. This is independent of the capacity of the incoming arcs, and no further priority rule is needed to obtain a unique solution.

By similar considerations as in the previous subsection, we obtain a solution also in the congested case if condition (27) is not verified. To obtain a unique solution, we modify the priority rule described in [7] for a traffic flow. Here, the main goal is to use the whole capacity of the outgoing arc $ i = 3 $, which implies

We set the merging parameter $ q \in [0, 1] $ such that

This leads to

describing a ratio of the actual fluxes on the corresponding arcs. The admissible region for the fluxes is

and shaded gray in Figure 6. If condition (29) is not fulfilled by only considering the ratio (31), the parameter $ q $ is adapted to obtain a unique solution. This is shown in Figure 6:

a) The intersection point $ P $ between the maximal outgoing flux (the line $ f_1+f_2 = f^{\text{max}}_3 $) and the priority ratio (31) is inside the admissible set $ \Theta $. In this case, we keep $ {q \in (0, 1)} $ fix and we have (30).

b) The intersection point $ P $ is outside $ \Theta $. We choose the closest point $ Q $ inside $ \Theta $ on the line $ f_1 + f_2 = f^{\text{max}}_3 $, which guarantees maximal throughput, i.e., $ f_1 = f^{\text{max}}_1 $. The merging parameter changes to $ q = f^{\text{max}}_1/ f^{\text{max}}_3. $ Then, the resulting fluxes are $ f_1 = f^{\text{max}}_1 $ and $ f_2 = f^{\text{max}}_3 - f^{\text{max}}_1. $

We call $ \Lambda_i $ the congested region (defined as in (20)) on arc $ i = 1, 2 $, $ g_i(t) $ its interface and $ q_{{i}} \in \{q, 1-q\} $ the corresponding merging parameters. The solution is then

If only one arc is congested, the solution holds true with $ \Lambda_i = \emptyset $ for the non-congested arc $ i $. The shape of the congested region $ \Lambda_i $ depends on the merging parameter $ q_{{i}} $. It is described by the interface $ g_i:[0, \infty) \rightarrow {\mathbb{R}}_{\leq 0} $, where $ t $ is mapped to the solution of equation

if this is negative, and to zero otherwise, analogously to the previous cases.

Remark. The time $ t_0 $ is unique since congestion starts (independent on $ q $) if the outgoing belt $ 3 $ is not able to absorb all the incoming flux $ f_1 + f_2 $. This is due to (29), which ensures that congestion arises if the maximal capacity is reached. At the same time the evolution of the function $ g_i(t) $ depends on the parameter $ q $. Therefore, it is possible that $ g_1(\hat t) = 0 $ for some $ \hat t>t_0 $ if $ g_2(\hat t)>0 $ (or vice-versa). This implies that for one arc $ i \in \{1, 2\} $, the function $ g_i $ may start from zero while the other one starts from a negative value, i.e., at least one arc is congested, and the description of the $ \Lambda_i $ is not completely separated.

We now draw our attention the numerical treatment of the formerly stated problems.

4.   Numerical approximation

In this section, we present numerical experiments to illustrate and confirm our theoretical results for different network configurations. The numerical scheme we propose is an adaptation of the scheme in [12] to networks. This scheme has also been successfully applied to pedestrian networks in [2].

We discretize each arc $ \Omega_i = ( \bar{\mathfrak{a}}_i, \bar{\mathfrak{b}}_i) \subset {\mathbb{R}} $ by $ X_i = (x_{i, 0} = \bar{\mathfrak{a}}_i, x_{i, 1}, \dots, x_{i, m_i} = \bar{\mathfrak{b}}_i) $ with constant discretization step $ \Delta x = {|x_{i, j} - x_{i, j-1}|}, \; j = 1, \ldots, { m_i} $. The spatial grid cells are defined as $ C_{i, j} = (x_{i, j-1/2}, x_{i, j+1/2} ) \subset \mathbb{R}, $ where the index $ i $ refers to the corresponding arc. We discretize the time set $ [0, T] $ with $ \{t^n = n\Delta t, \; n = 0, 1, \ldots, T/\Delta t\}, $ where $ \Delta t \in [0, T] $. We define the piecewise constant approximation of the solution $ \rho $ as $ \tilde{\rho}^n_{i, j}\approx\rho(x_{i, j}, t^n) $ with $ \tilde{\rho}^n $ constant in each grid cell $ C_{i, j} $. The scheme to update the approximation in each time step is for all $ i \in E, j \in {2, \ldots, m_i-1} $ given by

For $ j = 1 $ and all $ i \in E $, the update rule reads as

with $ h^{n, \text{in}}_i $ describing an inflow condition at time $ t^n $. Similarly, for $ j = m_i $ and outflow condition $ h^{n, \text{out}}_i $ at time $ t^n $, the update rule is for all $ i \in E $

For the nodes $ v \in V $ with $ \delta_v^- = \emptyset $, we need an inflow condition which is set to zero, i.e., $ h^{n, \text{in}}_{ \hat{i}} = 0 \text{ for } \hat{i} \in \delta_v^+. $ For the nodes $ v \in V $ with $ \delta_v^+ = \emptyset $, we set the outflow as $ h^{n, \text{out}}_{ i} = f_{ i}( \tilde{\rho}^n_{ i, m_{i}}) \text{ for } i \in \delta_v^-. $ For all interior nodes $ v \in V $ with $ \delta_v^- \neq \emptyset $ and $ \delta_v^+ \neq \emptyset $, we need to impose a junction rule at $ v $ depending on the type of junction. The inflow and outflow conditions are defined as follows.

a) For an one-to-one junction, i.e. $ |\delta_v^-| = |\delta_v^+| = 1 $, we set

b) For an one-to-two junction, i.e. $ 1 = |\delta_v^-|\neq|\delta_v^+| = 2 $, we choose in the non-congested case (18)

so that the flux is distributed according to the parameter $ \alpha_{ \hat{i}} $. In the congested case, i.e., if (18) is violated, we consider the junction properties described in Section 3.2. If we are in the case of a "passive" junction, we set the outgoing flux to

where the incoming fluxes on the outgoing arcs $ \hat{i} \in \delta_v^+ = \{2, 3\} $ are defined by

If we consider an "active" junction instead, we set

c) For a two-to-one junction, i.e. $ 2 = |\delta_v^-|\neq|\delta_v^+| = 1 $, we choose in the non-congested case (27)

In the congested case, we apply the merging parameter $ q_{{i}} $ described in Section 3.3 and set

To define $ h( \tilde{\rho}^n_{i, j}, \tilde{\rho}^n_{i, j+1}) $ in (34), we look for a function $ h $ satisfying

with continuous function $ m: {\mathbb{R}} \to {\mathbb{R}} $ and $ m_- = \min(m, 0) $, $ m_+ = \max (m, 0) $. Conditions (35) and (36) enable to use results from [2,12]. Since $ f_i $ is discontinuous in $ \rho $, we need a suitable regularization. We define a Friedrichs mollifier $ \varphi \in C_0^{\infty}( {\mathbb{R}}) $ with compact support in $ [-1, 1] $ such that

In our case, we use the mollifier $ {\varphi(y): = \max(0, 1-|y|)} $ and define $ \varphi_ \xi(y): = \frac{2}{ \xi}\varphi(\tfrac{2y}{ \xi}) $ for a small parameter $ \xi>0 $, which implies that $ \varphi_ \xi(y) $ has compact support in $ [-\frac{ \xi}{2}, \frac{ \xi}{2}] $. For each arc $ i \in E $, we introduce the following smooth regularization of the flux function (4)

see Figure 7. The function coincides with the original flux function (4) in $ x\in[0, \rho^{\text{max}}_i] $ but we observe that there is a continuously differentiable connection to the value $ {f_{ \xi, i}( \rho^{\text{max}}_i+ \xi) = 0} $, i.e., the regularized flux function $ f_{ \xi, i} $ itself is continuously differentiable. We note that $ (f_{ \xi, i}( \rho^{\text{max}}_i+ \xi))' = 0 $ and $ (f_{ \xi, i}(\bar \rho))' = \bar{a}, $ i.e. the transport velocity inside the congested area $ \Lambda $. Moreover, the derivative is bounded by $ |(f_{ \xi, i})'(\rho)|\leq \frac{2}{ \xi} $ for small $ \xi $. In the limit $ \xi\rightarrow 0^+ $, we recover the original discontinuous flux function (4).

Note that the second equality of condition (35) translates to

in the regularized case.

We choose the numerical flux function $ h $ as Godunov flux

and condition (36) is then satisfied with $ m(\rho) = \big(f_{ \xi, i}\big)'(\rho). $

The scheme (34) is stable, if the following CFL condition

is fulfilled, cf. [12]. From inequality (40), we can establish a relation between the regularization parameter $ \xi $ and the discretization steps $ \Delta t $ and $ \Delta x $. In particular, if $ \max|\delta_v^-|\leq 2 $ for a fixed $ v \in V $, the CFL condition reduces to

Knowing that $ \xi $ is supposed to be small, this is a quite restrictive condition. However, in the next section we will see that the choice of very small parameters $ \xi $ does not improve the solution significantly.

5.   Tests

In this section, we present numerical results for different network settings. Throughout this section, we compute the numerical solution based on scheme (34) with Godunov flux (39). If not stated otherwise, the space step size is $ \Delta x = 5\cdot 10^{-3} $, the time step size $ \Delta t = 10^{-5} $ and the smoothing parameter is fixed to $ \xi = 10^{-2}. $ For simplicity, we choose $ \rho^{\text{max}} = 1 $ in all experiments.

5.1.   One-to-one junction

First, we study the situation described in Section 3.1. The linear network is given by $ \Omega_1 = (-\infty, 0) $ and $ \Omega_2 = (0, +\infty) $, i.e. the coupling is at $ x = 0 $. We fix the initial solution $ \rho^0 $ on $ \Omega_1 \cup \Omega_2 $ as

Clearly, for computational purposes, we need to approximate the arcs $ \Omega_1 $, $ \Omega_2 $ with some bounded intervals. We notice that choosing $ \Omega_1\approx [-5, 0] $, $ \Omega_2\approx [0, 5] $, the solution has compact support in the approximated domain at every instant of its evolution. The same principle is applied elsewhere in the rest of the section.

Test 1: free-flow case. This is the non-congested case, i.e. condition (6) is holds true and the analytical solution is simply (7). The results are displayed in Figure 8. It shows the evolution of the initial density at different time steps for the different velocities $ a_{{1}} = 1 $ and $ a_{{2}} = 2 $. We see that the analytical and the numerical solution match very well. A discontinuity appears in the solution at $ x = 0 $. There, the doubling of the velocity has the effect of "spreading" the initial solution. Due to the mass conservation, the local density behind the junction point is halved. This can be seen by comparing the maximal arising density $ \rho = 1 $ in front of the junction and the maximal arising density behind the junction which is $ \rho = 0.5 $.

Test 2: congested case. For the second test, we only vary the velocity of the arcs and we set $ a_{{1}} = 2 $ and $ a_{{2}} = 1 $. Then, the condition (6) is not true anymore and we are in the congested case, where the analytical solution is given by (14). Figure 9 shows the comparison between the analytical and numerical solution. After the time $ t = 1.7 $, the evolution continues on arc $ 2 $ as linear transport. Obviously, the discontinuity backward wave, modeled by the function $ g(t) $ in (14), is correctly tracked by the numerical solution. It should be noticed that the density located in the congested region $ \Lambda $ in $ \Omega_1 $ moves with the velocity $ \bar{a} = a_{{2}} $, cf. (12).

Figure 10 shows the space-time diagram for the numerical and analytical solution. Note that the highlighted congested region $ \Lambda $ in the middle is correctly tracked by the numerical scheme. However, we observe the diffusive effect of the Godunov scheme. The latter effect could be reduced by the use of other flux approximations (as e.g. proposed in [6]) but not avoided completely.

For the same setting, the influence of the discretization step sizes is shown in Table 1 (left). Since the analytical solution is known, we can evaluate the $ L^2 $-error to the numerical solution. According to the space step size, the time step size is adapted to satisfy the CFL condition (40). As expected, the error tends to zero with decreasing step sizes.

To study the influence of the smoothing parameter $ \xi $, the fixed discretization is chosen such that the CFL condition (41) is fulfilled for the smallest value of $ \xi $. This leads to a time step $ {\Delta t = 2\cdot 10^{-6}} $ with a space step size $ \Delta x = 5\cdot 10^{-3} $. The result is shown in Table 1 (right). The error turns out to be only slightly smaller for the smallest value of $ \xi $.

5.2.   One-to-two junction

We pass to the one-to-two junction described in Section 3.2. We consider the domain $ \Omega_1 = (-5, 0)\times \{0\}, \; \{0\}\times\Omega_2 = (0, -5) , \; \Omega_3 = (0, 5)\times\{0\}, \; t\in [0, 2] $ with intersection point $ x = (0, 0) $. We choose the velocities $ a_{{1}} = 4, a_{{2}} = 1, a_{{3}} = 2 $ and distribution parameter $ \mu = 0.5 $. The initial solution $ \rho^0 $ is again (42). As already mentioned, the flux conservation (15) is not sufficient to ensure uniqueness of the solution at the intersection and therefore "passive" and "active" junctions are considered. Due to our choice of parameters, there exists a time $ t\in[0, 2] $ such that the conditions (18) and (22) are violated and thus congestion arises.

Test 3: "passive" junction. In the passive junction case, the solution is given by (21). The comparison of the numerical and analytical solution are displayed in Figure 11.

In the first column, the density distribution on arc $ 1 $ at different time steps is shown. In the second and the third column, the density distribution on arcs $ i = 2, 3 $ are presented. Congestion starts at about $ t_0 = 0.3 $, so at time $ t = 0.5 $, we are already in the congested phase. On arc $ i = 3 $, a density value of $ 0.5 $ is kept. This is due to the constant ratio of the outgoing fluxes, even if the maximal capacity is not used. As we see here, this leads to shocks in the solution also if the corresponding arc is not congested. At time $ t = 1.0 $, all mass passed the junction and is transported by the outgoing arcs. The final time of congestion is about $ t^E = 0.9 $ in this scenario.

Test 4: "active" junction. Here, the solution is given by (25). The results in Figure 12 are again for the congested case. Note that now congestion starts at about $ t_0 = 0.4 $. At time $ t = 0.5 $, the maximal capacity of both outgoing arcs $ i = 2, 3 $ is reached and we are in the congested phase. Compared to the passive junction case, congestion is reduced on the incoming arc $ i = 1 $. At time $ t = 1.0 $, all mass passed the junction and is transported by the outgoing arcs. Now, the final time of congestion is about $ t^E = 0.7 $ which is less than in the previous case.

5.3.   Two-to-one junction

The last test is the merging junction discussed in Section 3.3. We consider the domain $ \Omega_1 = (-5, 0)\times \{0\}, \; \Omega_2 = \{0\}\times(-5, 0), \; \Omega_3 = (0, 5)\times \{0\}, \; t\in [0, 5]. $ The initial data $ \rho^0_i $ on each incoming arc $ i = 1, 2 $ is (42). On the outgoing arc $ i = 3 $, we set $ \rho^0_3 = 0. $. All velocities are fixed to $ a_i = 1 $ for all arcs. This setting also allows to recover the results developed in [10], where the capacity and the speed are assumed to be equal for all arcs. Figure 13 shows the result of the evolution of the density at various time steps with merging parameter $ q = 0.3 $. The latter leads to a prioritization of arc $ 2 $ and a non-symmetric transportation on the two incoming arcs. We observe how the density initially placed on $ \Omega_1 $ and $ \Omega_2 $ is transported till $ x = (0, 0) $ is reached and congestion forms. At time $ t = 2 $, both arcs are congested. Due to the prioritization of arc $ 2 $, congestion is less than on arc $ 1 $. At time $ t = 4 $, all mass is absorbed by the outgoing arc. It is worth noticing that, once all density is absorbed by the outgoing arc, the configuration on the outgoing arc is the same, independent of the merging parameter $ q $. This is due to the condition (29), which implies that the outgoing arc always absorbs as much mass as possible.

To conclude, one can say that the numerical method presented in this section is a powerful tool to approximate the most meaningful solution of the material flow problem in the network case given by (3) - (4). This is particularly relevant in cases where we can no longer directly compute the analytical solution of the problem. To validate the results, in our numerical simulation study we considered special cases, for which we already derived the analytical solution. For those, the numerical results obtained match very well the analytical solution and catches the behavior of the evolution of the congested area before the junction point. Moreover, an error analysis for the case of a one-to-one junction has shown that the numerical solution converges against the analytical one with respect to the discretized version of the time-averaged $ L^2 $-norm.