Abstract
For the simulation of material flow problems based on two-dimensional hyperbolic partial differential equations different numerical methods can be applied. Compared to the widely used finite volume schemes we present an alternative approach, namely, the discontinuous Galerkin method, and explain how this method works within this framework. An extended numerical study is carried out comparing the finite volume and the discontinuous Galerkin approach concerning the quality of solutions.
1. Introduction
The modeling and simulation of material flow problems is motivated by engineering plant planning [1–3]. There exist various model approaches ranging from microscopic to macroscopic equations to describe the dynamic behavior of the manufacturing system; see, for example, [4–9], for an overview. The main difference is that microscopic models rely on systems of ordinary differential equations while macroscopic models are based on conservation laws. In particular for a large number of parts macroscopic models have the computational advantage that they are independent of individual parts. Therefore, with regard to the consideration of simulation and optimal control issues, macroscopic models are a beneficial tool for computation.
Macroscopic models, or conservation laws, are used in different engineering areas, for example, traffic flow [10], manufacturing systems [11], and crowd and evacuation dynamics [12, 13]. In the case of material flow problems, active research is recently done on a theoretical level [14] or on numerical solution algorithms [5, 15, 16].
In this work, we are concerned with a two-dimensional nonlocal hyperbolic conservation law that has been originally introduced [5]. This model has been validated against real data measurements and also provides reliable estimates on material flow and throughput rates in manufacturing system. It is even possible to rigorously derive such macroscopic models from microscopic ones via kinetic models; see [8, 9, 16]. From a numerical point of view, solution methods based on finite volume discretizations are an appropriate way to simulate material flow systems; see [5, 15]. However, due to the underlying geometry of the problem discontinuous Galerkin methods represent an alternative approach. This has been successfully shown in [17–23] for similar structured problems in one dimension. Our intention is now to derive a discontinuous Galerkin method in two space dimensions to solve the nonlocal hyperbolic conservation law from [5].
The paper is structured as follows. A short introduction of the macroscopic model taken from [5] is mentioned in Section 2. In Section 3 we discuss two numerical solution approaches for the macroscopic models. In detail, we shortly repeat a finite volume method with dimensional splitting and introduce a discontinuous Galerkin approach. Finally, numerical results are shown in Section 4. In particular, we qualitatively compare the different presented numerical methods.
2. Material Flow Modeling with Conservation Laws
We consider a conveyor belt with the given geometry presented in Figure 1 in two space dimensions, where an obstacle with angle in the middle of the belt is used to redirect parts. A large number of cylindrical-shaped parts are transported from the left to the right. Once the parts interact with the obstacle, queuing effects will occur. We furthermore assume that parts never move out of the , plane.
(a)
(b)
We briefly recall the conservation law for the evolution of the part density derived in [5]. The governing equations in two space dimensions, that is, , are thenwhere denotes the common Heaviside function that is zero for negative arguments and the fixed maximal density. The flux function is often referred to as .
Equation (1a) describes the evolution of the initial part density (1d) depending on velocity field composed of the dynamic velocity field and the time-independent velocity field . The field is induced by the movement of the conveyor belt. For the numerical simulations in Sections 3 and 4 we will use the velocity field as indicated in Figure 1, where the right picture represents a smoothed version of the velocity field. The latter is necessary to avoid problems of well-posedness and stability in the numerical simulations. Note that all angles between 0 and 90 degrees are feasible and that the smoothing is independent of the angle that is considered.
In contrast, the dynamic component in (1c) determines the movement of colliding parts, especially in the case they hit the obstacle. Since colliding parts do not penetrate each other, the density could not be larger than the density of a close-packing of parts . Therefore, we have to prevent situations, where for in a certain time . To activate the density dependent velocity , the Heaviside function in (1b) is introduced. If , that is, , the density dependent velocity is active and inactive otherwise. Hence, the velocity field disperses clouds with and parts are pushed into a direction with lower density.
The nonlocal operator in (1c) is controllable with the constant parameter Here, the negative gradient field is the steepest descent of the convolution , where is a sufficiently smooth mollifier or smoothing function, for example, a Gaussian. So the density dependent force term will only act in a small neighborhood of the boundary of a congested region.
The boundary conditions of (1a), (1b), (1c), and (1d) are imposed by the geometry of the belt (cf. Figure 1). We distinguish between solid boundaries (thick black lines), where free slip conditions are used, and the inflow region at the left boundary, where homogeneous Dirichlet conditions are applied.
3. Numerical Methods
Now we present suitable numerical methods to solve the conservation law (1a), (1b), (1c), and (1d). The first approach is based on a one-dimensional finite volume method which is extended into a two-dimensional problem solver by dimensional splitting [5]. The other approach is a discontinuous Galerkin method which is useful to compute accurate solutions on complex geometries. Since the finite volume method is validated against real data (cf. [5]) we will use the results of this approach as a benchmark solution in Section 4 to test the performance of the discontinuous Galerkin method.
For both simulation approaches, we assume that the discontinuous flux function in (1a) can be rewritten and approximated bywhere is a smoothed version of the Heaviside function; that is,
3.1. Finite Volume Approach with Dimensional Splitting
The following procedure is based on a finite volume method with dimensional splitting; see [24]. We use discretization of the two-dimensional spatial domain in rectangular cells where each cell is identified by the indices , . The center of a cell , is located at . The lengths of the cells are given by the spatial step sizes , and the time is discretized by step size . We use the following space and time grid:with cells . Additionally, grid constants are given by for . The density is defined as a step function Dimensional splitting is applied to split the original two-dimensional problem (1a), (1b), (1c), and (1d) into a sequence of one-dimensional problems. That means, for our purpose, the flux used in the numerical simulation is split in each dimension according to (2a). Applying a Godunov-type scheme, the numerical flux in one dimension, that is, , at points and , is a modified Roe flux combined with the nonlocal term :For at points and time the flux is determined analogously. The term or , respectively, including the convolution is evaluated as follows:whereFurthermore, the static flux in -direction iswhere the discretized static velocity field is given by Again, the flux in -direction is defined analogously.
Summarizing, we have to solve the coupled scheme under the stability condition (also known as CFL condition):for the smoothed Heaviside function (3). For more details and a detailed information on the algorithm we refer to [5].
3.2. Discontinuous Galerkin Methods
We now present an alternative approach to solve (1a), (1b), (1c), and (1d). Discontinuous Galerkin methods (DG methods) play an important role in finding approximations of many physical applications based on hyperbolic partial differential equations. For example, popular applications are found in gas dynamics, compressible and incompressible flows, chemical transports, granular flows, and more. We refer to [25–27] for a short overview. These methods have some interesting benefits; for example, they preserve the flexibility of finite elements in handling complicated geometries and they yield very accurate approximations. As already seen, finite volume methods use constant cell averages. In consideration of upwinding methods, this leads to artificial numerical diffusion which can influence the approximation quality. To avoid this drawback and for further investigations on optimal control issues, we consider other approximation tools such as the discontinuous Galerkin method.
In the following derivation, we assume that the flux function is approximated by polynomials. To ensure numerical stability of the DG method, we need to work again with a continuous flux function as already stated in (2a) and (3). The presentation of the DG method follows the ideas drawn in [18–21].
3.2.1. Space Integration
We consider a finite element discretization of the spatial domain , where is a disjoint union of triangle elements . Also, we assume that the position of each vertices of can only coincide with other vertices of neighboring triangle elements. An example of such finite element discretization or triangulation is given in Figure 2. Note that estimates the “size” of all triangle elements . In this work, denotes the length of the largest triangle edge of all elements .
Let be the solution space and let the approximate space be defined by where is the space of the polynomials of degree . By definition the solutions are discontinuous at the triangle interfaces. For the scheme we characterize all elements by a nodal basis. In this presentation, a nodal basis is a special case of a polynomial basis. Note that a two-dimensional polynomial has degrees of freedom for choosing the coefficients. All polynomials , restricted to a triangle-shaped domain , are constructible by nodal basis functions withwhere are nodal points on the finite element . The polynomials are called Lagrangian basis functions. The nodal points for are distributed on each triangle element as, respectively, shown in Figure 3.
(a)
(b)
(c)
An approximation of the solution (1a), (1b), (1c), and (1d) is given by an element of ; that is,The functions are unknowns and characterize the solution at time . We distinguish that the approximations and the flux fulfill (1a), (1b), (1c), and (1d) in an arbitrary way; that is, where is the residual. Generally, the approximation does not fulfill (1a), (1b), (1c), and (1d) exactly and the residual is not zero in all cases. Furthermore, we must decide in which sense the residual should vanish. Therefore, we choose a test function that is representable as We now require that the residual is orthogonal to all test functions in ; that is, This is true if and only if holds. Thus, we obtainIntegrating (21) by parts yieldswhere represents the local outward pointing normal. The solution at the interfaces between triangle elements is multiply-defined. At this moment, we have a lack of conditions on the local solution and the test functions. Therefore, we need here a correct combination of solutions to reduce the degrees of freedom. We select a numerical flux for the fluxes at the triangle interfaces. An illustrated example is given in Figure 4. Thus, (22) leads to the local statement:In this work, especially, we choose the local Lax-Friedrichs flux for the presented DG method: where , are the interior and exterior solution value. Respectively, is the local maximum of the directional flux The goal is to achieve an ODE system to obtain the quantity . We plug now (16) into (22) and we get the following statement:where denotes the outward pointing normal of the interface of the triangle . The ODE system (26) can be written into a matrix notation; that is,where is a vector of dimension containing the cell unknowns . The local mass matrices and the stiffness matrices , are defined by
Remark 1. The coefficient matrices , , for and depend only on the choice of the basis functions and the corresponding triangulation. Therefore, it is useful to compute these matrices once only for a complete simulation. This can be done by a preprocessing routine.
3.2.2. Discontinuous and Shock Solutions: Filtering
It is well known that nonlinear conservation laws might lead to shocks or discontinuities in solutions. However, the polynomial approximation of solutions of the DG method is not able to prescribe discontinuities so far. If we apply the previous DG method to problems with shock solutions, the following problems will occur:(i)The appearance of artificial and persistent oscillations around the point of discontinuity.(ii)The loss of pointwise convergence at the point of discontinuity.This phenomenon is already known as the Gibbs phenomenon and its behavior is well understood [28].
Note that a high order polynomial basis on the elements gives a high order accuracy of the scheme for smooth solutions. However, the DG method handles discontinuities with persistent oscillations that distort the approximate solution or influence the stability properties. Therefore, we propose the following filter approach in stabilizing the computations and in reducing the oscillations.
The filter approach [29, 30] considers ways to recover some accuracy information hidden in the oscillatory solutions. One possibility is filtering out high frequent redundant oscillations (high order polynomials) in the solutions.
In the following, we consider the canonical basiswhich spans the space of -dimensional polynomials in two variables . Additionally, the spatial variable is restricted to a reference triangle ; that is, . However, it is a complete polynomial basis and it can be orthonormalized through a Gram-Schmidt process. The resulting basis is denoted by . The next step is to transform the basis function back on a triangle element . This is realizable by a linear mapping . Thus, we obtain the basis function on by with the property An approximate solution of an element is given byThe solution above is given in a multidimensional Lagrange polynomial basis . Now we transform into the basis consisting of . Note that the polynomial has the degree . The idea of filtering is to reduce the coefficient of high polynomial basis elements . A popular choice is the exponential filterto obtain the filtered expansion where the filter is characterized by the the maximum damping parameter and the order . It is reasonable to use other filter approaches; see [29, 30]. In this work, we use a filter of the formThe filter (34) is an extension of the exponential filter (32). presents a cutoff; that is, polynomials with degree are left untouched. An example of the filter (34) with different parameters is shown in Figure 5.
(a) ,
(b) ,
Since filtering usage should be used both as minimal as possible and as much as needed, this is necessary to stabilize the method, reduce oscillatory solutions, and reduce artificial viscosity.
3.2.3. Convolution Integration
In particular, the dispersive term of (1a), (1b), (1c), and (1d) depends on the convolution of the density and the gradient of the mollifier ; that is, Hence, it is necessary to include the convolution process into the discontinuous Galerkin scheme. Without loss of generality, we consider the convolution of the approximate solution and in the nodal point of a triangle ; that is, The computation for the convolution of and works analogously.
Remark 2. Note that the weights are time independent. Therefore, can be computed once only before the simulation starts. However, the computation can result in high computational efforts for a large number of triangles and polynomial degree . Under certain circumstances, it is necessary to determine and store a number of weights to evaluate the convolution for all nodal points.
3.2.4. Time Integration
The DG approximation leads to a system of ordinary differential equations (ODEs) over each element . After inverting the local mass matrix , system (27) can be transformed in the following matrix form: where is a vector of dimension containing the cell unknowns . denotes the components of the right hand side of the ODE system (27) multiplied by the inverse mass matrix . The corresponding ODE system can be solved by explicit methods, for example, the forward Euler method.
As a result, the DG computation procedure is illustrated by the following steps:(1)Computation of is given as follows: (2)Reconstruction of the updated solution is given by applying where denotes the filter process that is discussed above.
4. Numerical Results
Finally, we present computational results comparing the finite volume approach and the discontinuous Galerkin method presented in Section 3. In particular, we comment on the quality of solutions and analyze lane and pattern artifacts.
All computations are performed on the same platform, namely, a 3.0 GHz Dual Core computer with 8 GB RAM, and all algorithms are implemented in MATLAB.
4.1. Finite Volume versus Discontinuous Galerkin
First, we compare the quality of the two methods to numerically solve (1a), (1b), (1c), and (1d). The finite volume method and the discontinuous Galerkin method offer their benefits as well as drawbacks that are independently discussed in this section.
4.1.1. Macroscopic Model Settings
The field is given by a smoothed velocity field as indicated in Figure 1(b), where the obstacle angle is set to 60 degrees. We choose a smoothed version of the Heaviside function with according to (3). The mollifier , occurring in the operator , is defined asIn this example, the maximal density is set to . The strength of the term is selected as , where the velocity of the conveyor belt is m/s. Furthermore, the total time horizon is .
4.1.2. Finite Volume Settings
The grid sizes of the finite volume approach with dimensional splitting are chosen as , in the following computation.
4.1.3. Discontinuous Galerkin Settings
The discontinuous Galerkin method uses a triangulation with a maximal triangle edge length . The polynomial degree of each finite element is . ODE system (27) is solved by the explicit Euler method with a time step size . Thereby the filter procedure is called in each computational step of the ODE solver. The filter settings are selected as , , and .
The results are shown in Figure 6. The left column shows the solution computed by the finite volume approach with splitting. The right column shows the results of the discontinuous Galerkin method. Each picture shows the density function as a gray-scaled image plot and each color specifies a density value. Thus, a dark color represents a higher density (black represents the maximal density) and vice versa. In all results, we observe that the parts are transported by the conveyor belt velocity . A formation of congestion can be seen in all results.
(a) Time
(b) Time
(c) Time
(d) Time
(e) Time
(f) Time
(g) Time
(h) Time
In all plots, we recognize a weak dispersing of density (cf. Figures 6(g) and 6(h)). This is caused by the term . The smoothed modification is never zero for . Consequently, the dispersing term is always activated and the quantity drifts apart all the time. This is also true, if the quantity has no connection to the singularizer; a dispersing effect is also recognizable (see Figures 6(a) and 6(b)). Moreover, the term disperses the quantity with addition of artifacts (lane formation). Indeed, lane formations are observable, for example, in Figure 6(g). The solution of the discontinuous Galerkin method seems to be smooth and not accurate in contrast to the results of the finite volume method. This is mainly due to the fact that the DG method uses polynomials on triangle finite elements of degree . However, polynomials are inherently smooth, and it is impossible to approximate accurate shock solutions due to the presented size of the finite elements. Indeed, the quality of the DG method can be improved by refining the triangle mesh grid. Compared to the DG method, the splitting method uses 20 times higher discretization.
The following question rises: what mesh grid sizes and what polynomial degrees are necessary to ensure good approximations due to the discontinuous Galerkin method. In the following, the previous example is computed again by the DG method with different triangulations and polynomial degrees. We test our problem on 3 different mesh grid sizes , , and . The results are shown in Figure 7. For all grid sizes and polynomial degrees, the qualitative behavior of the solution is approximated quite well. A finer grid or a higher polynomial degree generates more precise solutions; that is, quantity shocks and congested formations are drawn in an accurate way.
(a) ,
(b) ,
(c) ,
(d) ,
(e) ,
(f) ,
However, a rough triangulation or a low polynomial degree causes bad approximations (cf. Figure 7(e)). Compared to the other results, the congestion formation in Figure 7(e) looks quite degenerated.
The computation times of the DG method with respect to the mesh-sizes and polynomial degrees are shown in Tables 1 and 2. Furthermore, the computation times are distinguished into preprocessing time (cf. Table 2) and simulation time (cf. Table 1). Preprocessing contains the calculation of the coefficients of the convolution (see Remark 2). The simulation time contains the computation of ODE system (27) by the explicit Euler method.
The computing time required for the calculation of the finite volume approach is about 788.21 s. Consequently, the DG method is quite faster than the finite volume approach for all presented settings. However, the computing times and the memory requirements of the DG preprocessing increase enormously since the computation of the convolution in one nodal point requires at most coefficients. Furthermore, there are nodal points and the convolution is evaluated twice in each dimension. Thus, it is necessary to calculate and store about coefficients. As a consequence, the computer was not able to run the preprocessing routine successfully for small and a large ; for example, and ; see Table 2.
Let us summarize. The discontinuous Galerkin method is able to approximate accurately the extended flow equations on complex geometric domains. However, the presented example consists only of a rectangle-shaped domain and it is not necessary to use methods for complex geometries (cf. regular grids). As already seen, the DG method needs a very time and memory consuming preprocessing due to the convolution. Hence, it is very expensive to apply small step sizes for computation of accurate approximations and evaluating the corresponding convergence behavior.
4.2. Lane and Pattern Formation
The macroscopic model (1a), (1b), (1c), and (1d) is based on an integral-differential equation using a convolution term in the flux function . Similar models are already used for pedestrian flows [12], where certain lane or pattern artifacts are already observed. In this regard, it is not clearly understood why lane or pattern formation occurs. To investigate the phenomena of lane formation in more detail, we are motivated to study these artifacts by applying different numerical schemes.
Pedestrian models as [12] do not limit the influence of the dispersive term to a maximum density in (1a), (1b), (1c), and (1d) and therefore do not need the Heaviside function. If we additionally neglect the static velocity field in (1a), (1b), (1c), and (1d), the conservation law reduces toIn [12], lane formation was observed for the pedestrian model with smooth dispersive term, whereas this effect seems to be much less present in the above presented nonsmooth material flow model. Note that the finite volume method and the discontinuous Galerkin approach as well work with a smoothed Heaviside function. Therefore, we also observe lane artifacts in Figure 6(g).
We solve the simplified equation (41) on the spatial domain . The initial density is set to 1 for ; otherwise . Additionally, we compute the simplified model (41) for three different mollifiers; that is, . The step sizes of the finite volume approach are and .
The results of the finite volume approach are shown in Figure 8. In all plots, we observe that the quantity spreads out in all directions. Figures 8(a)–8(c) correspond to the setting with smoothing function parameter . We recognize a squared-shaped pattern in all time series. In Figures 8(d)–8(f) and 8(g)–8(i), the smoothing function parameter is used. Here, we observe a lane formation with a circular shape. A further increase of the mollifier parameter yields thinner lanes in the solution. However, we recognize the disappearance of the lanes in Figures 8(h) and 8(i). This is caused by the artificial numerical diffusion of the scheme which smears out the thin lanes in the solution.
(a) Time ,
(b) Time ,
(c) Time ,
(d) Time ,
(e) Time ,
(f) Time ,
(g) Time ,
(h) Time ,
(i) Time ,
Figure 9 shows the results of the discontinuous Galerkin method for different triangulations and polynomial degrees; however, the results are plotted for the time and . All plots (exceptional (a) and (d)) lead to the same result and they are similar to the plot of Figure 8(e). Indeed, a low triangulation and a low polynomial degree cause poor results (cf. Figures 9(a) and 9(d)). To get the most solution accuracy, the usage of filters for the DG computations is neglected. Therefore, some high frequent oscillations can appear (cf. Figure 8(c)).
(a) ,
(b) ,
(c) ,
(d) ,
(e) ,
(f) ,
(g) ,
(h) ,
(i) ,
5. Conclusion
We have presented a novel numerical simulation algorithm, the discontinuous Galerkin method, to compute the movement of material flow on conveyor belts. The numerical difficulties arise from the predefined geometry of the setting and the flux function consisting of a nonlocal term including a convolution. We have tested the performance of the discontinuous Galerkin method against a finite volume scheme and observed satisfactory results. In addition to the good qualitative behavior of the numerical results, we also detected and verified solution artifacts as lane formation in both numerical approaches.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was financially supported by the German Research Foundation (DFG), DFG grant OptiFlow (Project-ID GO 1920/3-1). Special thanks go to Veronika Schleper for fruitful and inspiring discussions.