# Discontinuous Galerkin Method for Material Flow Problems

Authors: Simone Göttlich; Patrick Schindler
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 Figure1 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 x1, x2 plane.Figure 1
(a) Geometry of the conveyor belt. (b) Static velocity field used for numerical simulations.
(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, x∈R2, are then(1a)∂tρ+∇·ρvdynρ+vstatx=0,(1b)vdynρ=Hρ-ρmax⁡·Iρ,(1c)Iρ=-ϵ∇η∗ρ1+∇η∗ρ22,(1d)ρx,0=ρ0x,x∈R2,where ρ=ρ(x,t),H denotes the common Heaviside function that is zero for negative arguments and ρmax⁡ the fixed maximal density. The flux function is often referred to as F(ρ)=(ρ(vdyn(ρ)+vstat(x))).Equation (1a) describes the evolution of the initial part density (1d) depending on velocity field composed of the dynamic velocity field vdyn(ρ) and the time-independent velocity field vstat(x). The field vstat(x) is induced by the movement of the conveyor belt. For the numerical simulations in Sections 3 and 4 we will use the velocity field vstat(x) 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 componentvdyn(ρ) 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 ρmax⁡. Therefore, we have to prevent situations, where ρ>ρmax⁡ for ρ0(x1,x2)<ρmax⁡ in a certain time t>0. To activate the density dependent velocity vdyn(ρ), the Heaviside function in (1b) is introduced. If ρ>ρmax⁡, that is, H(ρ-ρmax⁡)=1, the density dependent velocity is active and inactive otherwise. Hence, the velocity field vdyn(ρ) disperses clouds with ρ>ρmax⁡ and parts are pushed into a direction with lower density.The nonlocal operatorI(ρ) in (1c) is controllable with the constant parameter ϵ>0. 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 I(ρ) 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 by(2a)Fρ=F1ρ,x,F2ρ,xT=ρv1stat+v1dyn,ρv2stat+v2dynT,(2b)vddynρ=H~ρ-ρmax⁡Idρ,d=1,2,where H~ is a smoothed version of the Heaviside function; that is,(3)H~u=1πarctan⁡γu+12,γ>0.


### 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 i, j. The center of a cell i, j is located at xi,j=(x1,i,x2,j)T. The lengths of the cells are given by the spatial step sizes Δx1, Δx2 and the time t is discretized by step size Δt. We use the following space and time grid:(4)x1,i=iΔx1,i=1,…,Nx1,x2,j=jΔx2,j=1,…,Nx2,tk=kΔt,k=1,…,Ntwith cells Qi,j=[x1,i-1/2,x1,i+1/2]×[x2,j-1/2,x2,j+1/2]. Additionally, grid constants are given by λd=Δt/Δxd for d=1,2. The density ρ is defined as a step function (5)ρx,tk=ρi,jk∈Rforx∈Qi,j.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 ρ(vdyn(ρ)+vstat(x)) 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, d=1, at points xi+1/2,j and tk, is a modified Roe flux combined with the nonlocal term I(ρ)=(I1(ρ),I2(ρ))T:(6)F1ρ,ρi,jk,ρi+1,jk,xi+1/2,j=ρi,jkH~ρi,jk-ρmax⁡I1ρxi+1/2,j,I1ρxi+1/2,j≥0ρi+1,jkH~ρi+1,jk-ρmax⁡I1ρxi+1/2,j,I1ρxi+1/2,j≤0.For d=2 at points xi,j+1/2 and time tk the flux F2(ρ,ρijk,ρi,j+1k,xi,j+1/2) is determined analogously. The term I1(ρ) or I2(ρ), respectively, including the convolution is evaluated as follows:(7)I1ρxi+1/2,j=-ϵDx1ρi,j1+Dx1ρi,j2+Dx2ρi,j2,I2ρxi,j+1/2=-ϵDx2ρi,j1+Dx1ρi,j2+Dx2ρi,j2,where(8)Dx1ρi,j≔∑p,qρp,qk·∫Qp+1/2,q∂x1ητdτ,Dx2ρi,j≔∑p,qρp,qk·∫Qp,q+1/2∂x2ητdτ.Furthermore, the static flux in x1-direction is(9)G1ρi,jk,ρi+1,jk,vi+1/2,jstat=ρi,jkv1,i+1/2,jstat,v1,i+1/2,jstat≥0,ρi+1,jkv1,i+1/2,jstat,v1,i+1/2,jstat≤0,where the discretized static velocity field is given by (10)vi+1/2,jstat≔v1,i+1/2,jstat,v2,i+1/2,jstatT≔vstatxi+1/2,j.Again, the flux in x2-direction G2(ρi,jk,ρi,j+1k,vi,j+1/2stat) is defined analogously.Summarizing, we have to solve the coupled scheme(11)ρ~i,jk=ρi,jk-λ1F1+-F1-,F1+=F1ρ,ρi,jk,ρi+1,jk,xi+1/2,j+G1ρi,jk,ρi+1,jk,vi+1/2,jstat,F1-=F1ρ,ρi-1,jk,ρi,jk,xi-1/2,j+G1ρi-1,jk,ρi,jk,vi-1/2,jstat,ρi,jk+1=ρ~i,jk-λ2F2+-F2-,F2+=F2ρ,ρ~i,jk,ρ~i,j+1k,xi,j+1/2+G2ρ~i,jk,ρ~i,j+1k,vi,j+1/2stat,F2-=F2ρ,ρ~i,j-1k,ρ~i,jk,xi,j-1/2+G2ρ~i,j-1k,ρ~i,jk,vi,j-1/2statunder the stability condition (also known as CFL condition):(12)ΔtΔxdmaxρ⁡∂∂ρρvdynρ+vstatx∞≤1,ford=1,2for 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 functionF(ρ) 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Ω≃Ωh=⋃˙k=1KDk, where Ωh is a disjoint union of triangle elements Dk. Also, we assume that the position of each vertices of Dk 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 h estimates the “size” of all triangle elements Dk. In this work, h denotes the length of the largest triangle edge of all elements Dk.Figure 2
A finite element discretization (triangulation) of a domainΩ (ellipse).LetV=L2(Ω,R+) be the solution space and let the approximate space Vh⊂V be defined by (13)Vh≔v∈V:vDk∈PN,k=1,…,K,where PN is the space of the polynomials of degree N. By definition the solutions v are discontinuous at the triangle interfaces. For the scheme we characterize all elements v∈Vh 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 (14)Np≔N+1N+22degrees of freedom for choosing the coefficients. All polynomials vDk, restricted to a triangle-shaped domain Dk, are constructible by nodal basis functions lik(x)∈PN with(15)likxjk=1,i=j,0,i≠j,∀i,j=1,…,Np,where xjk∈Dk are nodal points on the finite element k. The polynomials lik(x) are called Lagrangian basis functions. The nodal points xik for i=1,…,Np are distributed on each triangle element Dk as, respectively, shown in Figure 3.Figure 3
Nodal points of the basis for linear, quadratic, and cubic triangle elementsDk.
(a)
N
=
1
(b)
N
=
2
(c)
N
=
3An approximation of the solution (1a), (1b), (1c), and (1d) is given by an element of Vh; that is,(16)ρhkx,t≔∑i=1Npρiktlikx,Fhkx,t≔∑i=1NpFρiktlikx,∀x∈Dk.The functions ρik(t) are unknowns and characterize the solution ρhk at time t. We distinguish that the approximations ρhk and the flux Fhk fulfill (1a), (1b), (1c), and (1d) in an arbitrary way; that is, (17)∂tρhkx,t+∇·Fhkx,t=Rhkx,t,∀x∈Dk,where Rhk(x,t) is the residual. Generally, the approximation ρhk 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 ϕ(x)∈Vh that is representable as (18)ϕhkx≔∑i=1Npϕiklikx,∀x∈Dk.We now require that the residual is orthogonal to all test functions in Vh; that is, (19)∫DkRhkx,tϕhkxdx=0.This is true if and only if (20)∫DkRhkx,tljkxdx=0,∀j=1,…,Npholds. Thus, we obtain(21)∫Dk∂tρhkx,t+∇·Fhkx,tljkxdx=0.Integrating (21) by parts yields(22)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·Fhkx,tljkxdx∀j=1,…,Np,where n 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 F∗ for the fluxes at the triangle interfaces. An illustrated example is given in Figure 4. Thus, (22) leads to the local statement:(23)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·F∗ljkxdx∀j=1,…,Np.In this work, especially, we choose the local Lax-Friedrichs flux for the presented DG method: (24)F∗ρ+,ρ-=Fρ++Fρ-2+C2nρ+-ρ-,where ρ+, ρ- are the interior and exterior solution value. Respectively, C is the local maximum of the directional flux (25)C=maxρ∈ρ+,ρ-⁡nx∂F1∂ρ+ny∂F2∂ρ.The goal is to achieve an ODE system to obtain the quantity ρik(t). We plug now (16) into (22) and we get the following statement:(26)∑i=1Np∂ρikt∂t∫Dklikxljkxdx-Fρikt·∫Dklikx∇ljkxdx=-∫∂Dkn·∑i=1NpF∗likxljkxdx=-∑e=13∫interfaceene·∑i=1NpF∗likxljkxdx,where ne denotes the outward pointing normal of the interface e of the triangle Dk. The ODE system (26) can be written into a matrix notation; that is,(27)Mk∂ρkt∂t+S1kF1ρkt+S2kF2ρkt=-∑e=13Mk,ene·F∗,where ρk is a vector of dimension Np containing the cell unknowns ρik. The local mass matrices Mk and the stiffness matrices S1k, S2k are defined by (28)Mi,jk=∫Dklikxljkxdx,Sd,i,jk=∫Dklikx∂xdljkxdx,d=1,2,∀i,j=1,…,Np,k=1,…,K,Mi,jk,e=∫interfaceelikxljkxdx,e=1,2,3.Figure 4
Interface of two neighboring trianglesDk and Dl. The position of the nodal points xik (red) and xjl (blue) coincides; that is, xik=xjl. The interior and exterior densities ρ+ and ρ- define the numerical flux F∗ at the transition point xik=xjl.Remark 1.
The coefficient matricesMi,jk, Sd,i,jk, Mi,jk,e for d=1,2 and e=1,2,3 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 basis(29)ψmr=r1ir2j,i,j≥0;i+j≤N,m≔j+N+1i+1-i2i-1,i,j≥0;i+j≤N,which spans the space of N-dimensional polynomials in two variables r=(r1,r2). Additionally, the spatial variable r is restricted to a reference triangle I; that is, r∈I≔{(r1,r2):r1,r2≥-1,r1+r2≤0}. However, it is a complete polynomial basis and it can be orthonormalized through a Gram-Schmidt process. The resulting basis is denoted by ψ~m(r). The next step is to transform the basis function ψ~m(r) back on a triangle element Dk. This is realizable by a linear mapping Ψ:I→Dk. Thus, we obtain the basis function on Dk by ψ~mk(x)≔ψ~(Ψ-1(x)) with the property (30)∫Dkψ~mxψ~nxdx=δm,n.An approximate solution of an element Dk is given by(31)ρhkx=∑i=1Npρiklikx=∑m=1Npρ~mkψ~mk.The solution above is given in a multidimensional Lagrange polynomial basis lik. Now we transform ρhk(x) into the basis consisting of ψ~mk. Note that the polynomial ψ~mk has the degree deg(ψ~mk)=i+j. The idea of filtering is to reduce the coefficient ρ~mk of high polynomial basis elements ψ~mk. A popular choice is the exponential filter(32)ςω=exp⁡-βω2sto obtain the filtered expansion (33)ρhk,Fx=∑i,j≥0i+j≤Nςi+jNρ~mkψ~mk,where the filter is characterized by the the maximum damping parameter β>0 and the order s>0. It is reasonable to use other filter approaches; see [29, 30]. In this work, we use a filter of the form(34)ςω=1,0≤ω≤ωc=NcN,exp⁡-βω-ωc1-ωc2s,ωc≤ω≤1.The filter (34) is an extension of the exponential filter (32). Nc presents a cutoff; that is, polynomials ρ~mk with degree deg(ρ~mk)≤Nc are left untouched. An example of the filter (34) with different parameters is shown in Figure 5.Figure 5
Examples of how the filter function (34) varies from the three parameters: the order s, the cutoff Nc=Nωc, and the maximum damping parameter β.
(a)
β
=
36, ωc=0.5
(b)
s
=
2, ωc=0.25Since 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 termI(ρ) of (1a), (1b), (1c), and (1d) depends on the convolution of the density ρ and the gradient of the mollifier η; that is, (35)∇η∗ρ=∂x1η∗ρ,∂x2η∗ρT.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 ρh∈Vh and ∂1η in the nodal point xik of a triangle k; that is, (36)∂x1η∗ρhxik=∫Ωηxik-τρhτdτ=∑l=1K∫Dlηxik-τρhlτdτ=∑l=1K∫Dlηxik-τ∑j=1Npρjlljlτdτ=∑l=1K∑j=1Npρjl∫Dlηxik-τljlτdτ︸≔ci,jk,l=∑l=1K∑j=1Npρjlci,jk,l.The computation for the convolution of ρh∈Vh and ∂2η works analogously.Remark 2.
Note that the weightsci,jk,l are time independent. Therefore, ci,jk,l can be computed once only before the simulation starts. However, the computation can result in high computational efforts for a large number of triangles K and polynomial degree N. Under certain circumstances, it is necessary to determine and store a number of O((NpK)2) weights to evaluate the convolution (∂x1η∗ρh) for all nodal points.


#### 3.2.4. Time Integration

The DG approximation leads to a system ofNp ordinary differential equations (ODEs) over each element Dk. After inverting the local mass matrix Mk, system (27) can be transformed in the following matrix form: (37)dρktdt=Aρk,where ρk(t) is a vector of dimension Np containing the cell unknowns ρik. A(ρk) denotes the components of the right hand side of the ODE system (27) multiplied by the inverse mass matrix Mi,jk. 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ρ~k is given as follows: (38)ρ~k=ρktn+ΔtAρktn,∀k=1,…,K.(2)
Reconstruction of the updated solutionρ~k is given by applying (39)ρktn+1=Uρ~k,∀k=1,…,K,whereU denotes the filter process that is discussed above.


## 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 i, j. The center of a cell i, j is located at xi,j=(x1,i,x2,j)T. The lengths of the cells are given by the spatial step sizes Δx1, Δx2 and the time t is discretized by step size Δt. We use the following space and time grid:(4)x1,i=iΔx1,i=1,…,Nx1,x2,j=jΔx2,j=1,…,Nx2,tk=kΔt,k=1,…,Ntwith cells Qi,j=[x1,i-1/2,x1,i+1/2]×[x2,j-1/2,x2,j+1/2]. Additionally, grid constants are given by λd=Δt/Δxd for d=1,2. The density ρ is defined as a step function (5)ρx,tk=ρi,jk∈Rforx∈Qi,j.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 ρ(vdyn(ρ)+vstat(x)) 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, d=1, at points xi+1/2,j and tk, is a modified Roe flux combined with the nonlocal term I(ρ)=(I1(ρ),I2(ρ))T:(6)F1ρ,ρi,jk,ρi+1,jk,xi+1/2,j=ρi,jkH~ρi,jk-ρmax⁡I1ρxi+1/2,j,I1ρxi+1/2,j≥0ρi+1,jkH~ρi+1,jk-ρmax⁡I1ρxi+1/2,j,I1ρxi+1/2,j≤0.For d=2 at points xi,j+1/2 and time tk the flux F2(ρ,ρijk,ρi,j+1k,xi,j+1/2) is determined analogously. The term I1(ρ) or I2(ρ), respectively, including the convolution is evaluated as follows:(7)I1ρxi+1/2,j=-ϵDx1ρi,j1+Dx1ρi,j2+Dx2ρi,j2,I2ρxi,j+1/2=-ϵDx2ρi,j1+Dx1ρi,j2+Dx2ρi,j2,where(8)Dx1ρi,j≔∑p,qρp,qk·∫Qp+1/2,q∂x1ητdτ,Dx2ρi,j≔∑p,qρp,qk·∫Qp,q+1/2∂x2ητdτ.Furthermore, the static flux in x1-direction is(9)G1ρi,jk,ρi+1,jk,vi+1/2,jstat=ρi,jkv1,i+1/2,jstat,v1,i+1/2,jstat≥0,ρi+1,jkv1,i+1/2,jstat,v1,i+1/2,jstat≤0,where the discretized static velocity field is given by (10)vi+1/2,jstat≔v1,i+1/2,jstat,v2,i+1/2,jstatT≔vstatxi+1/2,j.Again, the flux in x2-direction G2(ρi,jk,ρi,j+1k,vi,j+1/2stat) is defined analogously.Summarizing, we have to solve the coupled scheme(11)ρ~i,jk=ρi,jk-λ1F1+-F1-,F1+=F1ρ,ρi,jk,ρi+1,jk,xi+1/2,j+G1ρi,jk,ρi+1,jk,vi+1/2,jstat,F1-=F1ρ,ρi-1,jk,ρi,jk,xi-1/2,j+G1ρi-1,jk,ρi,jk,vi-1/2,jstat,ρi,jk+1=ρ~i,jk-λ2F2+-F2-,F2+=F2ρ,ρ~i,jk,ρ~i,j+1k,xi,j+1/2+G2ρ~i,jk,ρ~i,j+1k,vi,j+1/2stat,F2-=F2ρ,ρ~i,j-1k,ρ~i,jk,xi,j-1/2+G2ρ~i,j-1k,ρ~i,jk,vi,j-1/2statunder the stability condition (also known as CFL condition):(12)ΔtΔxdmaxρ⁡∂∂ρρvdynρ+vstatx∞≤1,ford=1,2for 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 functionF(ρ) 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Ω≃Ωh=⋃˙k=1KDk, where Ωh is a disjoint union of triangle elements Dk. Also, we assume that the position of each vertices of Dk 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 h estimates the “size” of all triangle elements Dk. In this work, h denotes the length of the largest triangle edge of all elements Dk.Figure 2
A finite element discretization (triangulation) of a domainΩ (ellipse).LetV=L2(Ω,R+) be the solution space and let the approximate space Vh⊂V be defined by (13)Vh≔v∈V:vDk∈PN,k=1,…,K,where PN is the space of the polynomials of degree N. By definition the solutions v are discontinuous at the triangle interfaces. For the scheme we characterize all elements v∈Vh 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 (14)Np≔N+1N+22degrees of freedom for choosing the coefficients. All polynomials vDk, restricted to a triangle-shaped domain Dk, are constructible by nodal basis functions lik(x)∈PN with(15)likxjk=1,i=j,0,i≠j,∀i,j=1,…,Np,where xjk∈Dk are nodal points on the finite element k. The polynomials lik(x) are called Lagrangian basis functions. The nodal points xik for i=1,…,Np are distributed on each triangle element Dk as, respectively, shown in Figure 3.Figure 3
Nodal points of the basis for linear, quadratic, and cubic triangle elementsDk.
(a)
N
=
1
(b)
N
=
2
(c)
N
=
3An approximation of the solution (1a), (1b), (1c), and (1d) is given by an element of Vh; that is,(16)ρhkx,t≔∑i=1Npρiktlikx,Fhkx,t≔∑i=1NpFρiktlikx,∀x∈Dk.The functions ρik(t) are unknowns and characterize the solution ρhk at time t. We distinguish that the approximations ρhk and the flux Fhk fulfill (1a), (1b), (1c), and (1d) in an arbitrary way; that is, (17)∂tρhkx,t+∇·Fhkx,t=Rhkx,t,∀x∈Dk,where Rhk(x,t) is the residual. Generally, the approximation ρhk 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 ϕ(x)∈Vh that is representable as (18)ϕhkx≔∑i=1Npϕiklikx,∀x∈Dk.We now require that the residual is orthogonal to all test functions in Vh; that is, (19)∫DkRhkx,tϕhkxdx=0.This is true if and only if (20)∫DkRhkx,tljkxdx=0,∀j=1,…,Npholds. Thus, we obtain(21)∫Dk∂tρhkx,t+∇·Fhkx,tljkxdx=0.Integrating (21) by parts yields(22)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·Fhkx,tljkxdx∀j=1,…,Np,where n 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 F∗ for the fluxes at the triangle interfaces. An illustrated example is given in Figure 4. Thus, (22) leads to the local statement:(23)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·F∗ljkxdx∀j=1,…,Np.In this work, especially, we choose the local Lax-Friedrichs flux for the presented DG method: (24)F∗ρ+,ρ-=Fρ++Fρ-2+C2nρ+-ρ-,where ρ+, ρ- are the interior and exterior solution value. Respectively, C is the local maximum of the directional flux (25)C=maxρ∈ρ+,ρ-⁡nx∂F1∂ρ+ny∂F2∂ρ.The goal is to achieve an ODE system to obtain the quantity ρik(t). We plug now (16) into (22) and we get the following statement:(26)∑i=1Np∂ρikt∂t∫Dklikxljkxdx-Fρikt·∫Dklikx∇ljkxdx=-∫∂Dkn·∑i=1NpF∗likxljkxdx=-∑e=13∫interfaceene·∑i=1NpF∗likxljkxdx,where ne denotes the outward pointing normal of the interface e of the triangle Dk. The ODE system (26) can be written into a matrix notation; that is,(27)Mk∂ρkt∂t+S1kF1ρkt+S2kF2ρkt=-∑e=13Mk,ene·F∗,where ρk is a vector of dimension Np containing the cell unknowns ρik. The local mass matrices Mk and the stiffness matrices S1k, S2k are defined by (28)Mi,jk=∫Dklikxljkxdx,Sd,i,jk=∫Dklikx∂xdljkxdx,d=1,2,∀i,j=1,…,Np,k=1,…,K,Mi,jk,e=∫interfaceelikxljkxdx,e=1,2,3.Figure 4
Interface of two neighboring trianglesDk and Dl. The position of the nodal points xik (red) and xjl (blue) coincides; that is, xik=xjl. The interior and exterior densities ρ+ and ρ- define the numerical flux F∗ at the transition point xik=xjl.Remark 1.
The coefficient matricesMi,jk, Sd,i,jk, Mi,jk,e for d=1,2 and e=1,2,3 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 basis(29)ψmr=r1ir2j,i,j≥0;i+j≤N,m≔j+N+1i+1-i2i-1,i,j≥0;i+j≤N,which spans the space of N-dimensional polynomials in two variables r=(r1,r2). Additionally, the spatial variable r is restricted to a reference triangle I; that is, r∈I≔{(r1,r2):r1,r2≥-1,r1+r2≤0}. However, it is a complete polynomial basis and it can be orthonormalized through a Gram-Schmidt process. The resulting basis is denoted by ψ~m(r). The next step is to transform the basis function ψ~m(r) back on a triangle element Dk. This is realizable by a linear mapping Ψ:I→Dk. Thus, we obtain the basis function on Dk by ψ~mk(x)≔ψ~(Ψ-1(x)) with the property (30)∫Dkψ~mxψ~nxdx=δm,n.An approximate solution of an element Dk is given by(31)ρhkx=∑i=1Npρiklikx=∑m=1Npρ~mkψ~mk.The solution above is given in a multidimensional Lagrange polynomial basis lik. Now we transform ρhk(x) into the basis consisting of ψ~mk. Note that the polynomial ψ~mk has the degree deg(ψ~mk)=i+j. The idea of filtering is to reduce the coefficient ρ~mk of high polynomial basis elements ψ~mk. A popular choice is the exponential filter(32)ςω=exp⁡-βω2sto obtain the filtered expansion (33)ρhk,Fx=∑i,j≥0i+j≤Nςi+jNρ~mkψ~mk,where the filter is characterized by the the maximum damping parameter β>0 and the order s>0. It is reasonable to use other filter approaches; see [29, 30]. In this work, we use a filter of the form(34)ςω=1,0≤ω≤ωc=NcN,exp⁡-βω-ωc1-ωc2s,ωc≤ω≤1.The filter (34) is an extension of the exponential filter (32). Nc presents a cutoff; that is, polynomials ρ~mk with degree deg(ρ~mk)≤Nc are left untouched. An example of the filter (34) with different parameters is shown in Figure 5.Figure 5
Examples of how the filter function (34) varies from the three parameters: the order s, the cutoff Nc=Nωc, and the maximum damping parameter β.
(a)
β
=
36, ωc=0.5
(b)
s
=
2, ωc=0.25Since 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 termI(ρ) of (1a), (1b), (1c), and (1d) depends on the convolution of the density ρ and the gradient of the mollifier η; that is, (35)∇η∗ρ=∂x1η∗ρ,∂x2η∗ρT.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 ρh∈Vh and ∂1η in the nodal point xik of a triangle k; that is, (36)∂x1η∗ρhxik=∫Ωηxik-τρhτdτ=∑l=1K∫Dlηxik-τρhlτdτ=∑l=1K∫Dlηxik-τ∑j=1Npρjlljlτdτ=∑l=1K∑j=1Npρjl∫Dlηxik-τljlτdτ︸≔ci,jk,l=∑l=1K∑j=1Npρjlci,jk,l.The computation for the convolution of ρh∈Vh and ∂2η works analogously.Remark 2.
Note that the weightsci,jk,l are time independent. Therefore, ci,jk,l can be computed once only before the simulation starts. However, the computation can result in high computational efforts for a large number of triangles K and polynomial degree N. Under certain circumstances, it is necessary to determine and store a number of O((NpK)2) weights to evaluate the convolution (∂x1η∗ρh) for all nodal points.


### 3.2.4. Time Integration

The DG approximation leads to a system ofNp ordinary differential equations (ODEs) over each element Dk. After inverting the local mass matrix Mk, system (27) can be transformed in the following matrix form: (37)dρktdt=Aρk,where ρk(t) is a vector of dimension Np containing the cell unknowns ρik. A(ρk) denotes the components of the right hand side of the ODE system (27) multiplied by the inverse mass matrix Mi,jk. 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ρ~k is given as follows: (38)ρ~k=ρktn+ΔtAρktn,∀k=1,…,K.(2)
Reconstruction of the updated solutionρ~k is given by applying (39)ρktn+1=Uρ~k,∀k=1,…,K,whereU denotes the filter process that is discussed above.


## 3.2.1. Space Integration

We consider a finite element discretization of the spatial domainΩ≃Ωh=⋃˙k=1KDk, where Ωh is a disjoint union of triangle elements Dk. Also, we assume that the position of each vertices of Dk 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 h estimates the “size” of all triangle elements Dk. In this work, h denotes the length of the largest triangle edge of all elements Dk.Figure 2
A finite element discretization (triangulation) of a domainΩ (ellipse).LetV=L2(Ω,R+) be the solution space and let the approximate space Vh⊂V be defined by (13)Vh≔v∈V:vDk∈PN,k=1,…,K,where PN is the space of the polynomials of degree N. By definition the solutions v are discontinuous at the triangle interfaces. For the scheme we characterize all elements v∈Vh 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 (14)Np≔N+1N+22degrees of freedom for choosing the coefficients. All polynomials vDk, restricted to a triangle-shaped domain Dk, are constructible by nodal basis functions lik(x)∈PN with(15)likxjk=1,i=j,0,i≠j,∀i,j=1,…,Np,where xjk∈Dk are nodal points on the finite element k. The polynomials lik(x) are called Lagrangian basis functions. The nodal points xik for i=1,…,Np are distributed on each triangle element Dk as, respectively, shown in Figure 3.Figure 3
Nodal points of the basis for linear, quadratic, and cubic triangle elementsDk.
(a)
N
=
1
(b)
N
=
2
(c)
N
=
3An approximation of the solution (1a), (1b), (1c), and (1d) is given by an element of Vh; that is,(16)ρhkx,t≔∑i=1Npρiktlikx,Fhkx,t≔∑i=1NpFρiktlikx,∀x∈Dk.The functions ρik(t) are unknowns and characterize the solution ρhk at time t. We distinguish that the approximations ρhk and the flux Fhk fulfill (1a), (1b), (1c), and (1d) in an arbitrary way; that is, (17)∂tρhkx,t+∇·Fhkx,t=Rhkx,t,∀x∈Dk,where Rhk(x,t) is the residual. Generally, the approximation ρhk 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 ϕ(x)∈Vh that is representable as (18)ϕhkx≔∑i=1Npϕiklikx,∀x∈Dk.We now require that the residual is orthogonal to all test functions in Vh; that is, (19)∫DkRhkx,tϕhkxdx=0.This is true if and only if (20)∫DkRhkx,tljkxdx=0,∀j=1,…,Npholds. Thus, we obtain(21)∫Dk∂tρhkx,t+∇·Fhkx,tljkxdx=0.Integrating (21) by parts yields(22)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·Fhkx,tljkxdx∀j=1,…,Np,where n 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 F∗ for the fluxes at the triangle interfaces. An illustrated example is given in Figure 4. Thus, (22) leads to the local statement:(23)∫Dk∂ρhkx,t∂tljkx-Fhkx,t·∇ljkxdx=-∫∂Dkn·F∗ljkxdx∀j=1,…,Np.In this work, especially, we choose the local Lax-Friedrichs flux for the presented DG method: (24)F∗ρ+,ρ-=Fρ++Fρ-2+C2nρ+-ρ-,where ρ+, ρ- are the interior and exterior solution value. Respectively, C is the local maximum of the directional flux (25)C=maxρ∈ρ+,ρ-⁡nx∂F1∂ρ+ny∂F2∂ρ.The goal is to achieve an ODE system to obtain the quantity ρik(t). We plug now (16) into (22) and we get the following statement:(26)∑i=1Np∂ρikt∂t∫Dklikxljkxdx-Fρikt·∫Dklikx∇ljkxdx=-∫∂Dkn·∑i=1NpF∗likxljkxdx=-∑e=13∫interfaceene·∑i=1NpF∗likxljkxdx,where ne denotes the outward pointing normal of the interface e of the triangle Dk. The ODE system (26) can be written into a matrix notation; that is,(27)Mk∂ρkt∂t+S1kF1ρkt+S2kF2ρkt=-∑e=13Mk,ene·F∗,where ρk is a vector of dimension Np containing the cell unknowns ρik. The local mass matrices Mk and the stiffness matrices S1k, S2k are defined by (28)Mi,jk=∫Dklikxljkxdx,Sd,i,jk=∫Dklikx∂xdljkxdx,d=1,2,∀i,j=1,…,Np,k=1,…,K,Mi,jk,e=∫interfaceelikxljkxdx,e=1,2,3.Figure 4
Interface of two neighboring trianglesDk and Dl. The position of the nodal points xik (red) and xjl (blue) coincides; that is, xik=xjl. The interior and exterior densities ρ+ and ρ- define the numerical flux F∗ at the transition point xik=xjl.Remark 1.
The coefficient matricesMi,jk, Sd,i,jk, Mi,jk,e for d=1,2 and e=1,2,3 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 basis(29)ψmr=r1ir2j,i,j≥0;i+j≤N,m≔j+N+1i+1-i2i-1,i,j≥0;i+j≤N,which spans the space of N-dimensional polynomials in two variables r=(r1,r2). Additionally, the spatial variable r is restricted to a reference triangle I; that is, r∈I≔{(r1,r2):r1,r2≥-1,r1+r2≤0}. However, it is a complete polynomial basis and it can be orthonormalized through a Gram-Schmidt process. The resulting basis is denoted by ψ~m(r). The next step is to transform the basis function ψ~m(r) back on a triangle element Dk. This is realizable by a linear mapping Ψ:I→Dk. Thus, we obtain the basis function on Dk by ψ~mk(x)≔ψ~(Ψ-1(x)) with the property (30)∫Dkψ~mxψ~nxdx=δm,n.An approximate solution of an element Dk is given by(31)ρhkx=∑i=1Npρiklikx=∑m=1Npρ~mkψ~mk.The solution above is given in a multidimensional Lagrange polynomial basis lik. Now we transform ρhk(x) into the basis consisting of ψ~mk. Note that the polynomial ψ~mk has the degree deg(ψ~mk)=i+j. The idea of filtering is to reduce the coefficient ρ~mk of high polynomial basis elements ψ~mk. A popular choice is the exponential filter(32)ςω=exp⁡-βω2sto obtain the filtered expansion (33)ρhk,Fx=∑i,j≥0i+j≤Nςi+jNρ~mkψ~mk,where the filter is characterized by the the maximum damping parameter β>0 and the order s>0. It is reasonable to use other filter approaches; see [29, 30]. In this work, we use a filter of the form(34)ςω=1,0≤ω≤ωc=NcN,exp⁡-βω-ωc1-ωc2s,ωc≤ω≤1.The filter (34) is an extension of the exponential filter (32). Nc presents a cutoff; that is, polynomials ρ~mk with degree deg(ρ~mk)≤Nc are left untouched. An example of the filter (34) with different parameters is shown in Figure 5.Figure 5
Examples of how the filter function (34) varies from the three parameters: the order s, the cutoff Nc=Nωc, and the maximum damping parameter β.
(a)
β
=
36, ωc=0.5
(b)
s
=
2, ωc=0.25Since 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 termI(ρ) of (1a), (1b), (1c), and (1d) depends on the convolution of the density ρ and the gradient of the mollifier η; that is, (35)∇η∗ρ=∂x1η∗ρ,∂x2η∗ρT.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 ρh∈Vh and ∂1η in the nodal point xik of a triangle k; that is, (36)∂x1η∗ρhxik=∫Ωηxik-τρhτdτ=∑l=1K∫Dlηxik-τρhlτdτ=∑l=1K∫Dlηxik-τ∑j=1Npρjlljlτdτ=∑l=1K∑j=1Npρjl∫Dlηxik-τljlτdτ︸≔ci,jk,l=∑l=1K∑j=1Npρjlci,jk,l.The computation for the convolution of ρh∈Vh and ∂2η works analogously.Remark 2.
Note that the weightsci,jk,l are time independent. Therefore, ci,jk,l can be computed once only before the simulation starts. However, the computation can result in high computational efforts for a large number of triangles K and polynomial degree N. Under certain circumstances, it is necessary to determine and store a number of O((NpK)2) weights to evaluate the convolution (∂x1η∗ρh) for all nodal points.


## 3.2.4. Time Integration

The DG approximation leads to a system ofNp ordinary differential equations (ODEs) over each element Dk. After inverting the local mass matrix Mk, system (27) can be transformed in the following matrix form: (37)dρktdt=Aρk,where ρk(t) is a vector of dimension Np containing the cell unknowns ρik. A(ρk) denotes the components of the right hand side of the ODE system (27) multiplied by the inverse mass matrix Mi,jk. 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ρ~k is given as follows: (38)ρ~k=ρktn+ΔtAρktn,∀k=1,…,K.(2)
Reconstruction of the updated solutionρ~k is given by applying (39)ρktn+1=Uρ~k,∀k=1,…,K,whereU 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 Section3. 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 fieldvstat(x) 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 γ=25 according to (3). The mollifier η, occurring in the operator I(ρ), is defined as(40)ηx=σ2πexp⁡-12σx22,σ=2500.In this example, the maximal density is set to ρmax⁡=1. The strength of the term I(ρ) is selected as ϵ=2vT, where the velocity of the conveyor belt is vT=0.395m/s. Furthermore, the total time horizon is T=7.


#### 4.1.2. Finite Volume Settings

The grid sizes of the finite volume approach with dimensional splitting are chosen asΔx1=Δx2=5·10-3, Δt=1.25·10-3 in the following computation.


#### 4.1.3. Discontinuous Galerkin Settings

The discontinuous Galerkin method uses a triangulationΩh with a maximal triangle edge length h=0.1. The polynomial degree of each finite element is N=10. ODE system (27) is solved by the explicit Euler method with a time step size Δt=10-3. Thereby the filter procedure is called in each computational step of the ODE solver. The filter settings are selected as β=36, s=6, and Nc=1.The results are shown in Figure6. 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 vT. A formation of congestion can be seen in all results.Figure 6
Results of the finite volume method with splitting:Δx1=Δx2=5⋅10-3 (left) and results of the discontinuous Galerkin method: h≈0.1, N=10 (right).
(a)
Timet=0.50
(b)
Timet=0.50
(c)
Timet=1.00
(d)
Timet=1.00
(e)
Timet=1.50
(f)
Timet=1.50
(g)
Timet=2.00
(h)
Timet=2.00In all plots, we recognize a weak dispersing of density (cf. Figures6(g) and 6(h)). This is caused by the term vdyn=H~(ρ-ρmax⁡)I(ρ). The smoothed modification H~(ρ-ρmax⁡) is never zero for ρ<ρmax⁡. Consequently, the dispersing term I(ρ) 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 I(ρ) 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 N=10. 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 sizesh=0.1, h=0.06, and h=0.04. 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.Figure 7
Results of the discontinuous Galerkin method with different triangulations (h=0.1,0.06,0.04) and polynomials degrees N. All plots show the solution at time t=1.
(a)
h
≈
0.1, N=5
(b)
h
≈
0.1, N=11
(c)
h
≈
0.06, N=3
(d)
h
≈
0.06, N=7
(e)
h
≈
0.04, N=1
(f)
h
≈
0.04, N=3However, a rough triangulation or a low polynomial degree causes bad approximations (cf. Figure7(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 Tables1 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.Table 1
Computation times of the discontinuous Galerkin method (simulation process) with different grid sizesh and polynomial degrees N. The time is measured in seconds.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
7.14 s
12.30 s
13.42 s
3
9.30 s
18.63 s
51.70 s
5
17.31 s
46.29 s
—
7
30.10 s
111.94 s
—
9
49.58 s
—
—
11
88.70 s
—
—Table 2
Computation times in seconds for the convolution preprocessing due to the grid sizeh and polynomial degree N.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
0.06 s
1.12 s
1.73 s
3
0.48 s
3.88 s
79.99 s
5
2.01 s
60.82 s
—
7
5.44 s
420.08 s
—
9
47.11 s
—
—
11
183.86 s
—
—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 mostNp·K coefficients. Furthermore, there are Np·K nodal points and the convolution is evaluated twice in each dimension. Thus, it is necessary to calculate and store about 2·(Np·K)2 coefficients. As a consequence, the computer was not able to run the preprocessing routine successfully for small h and a large N; for example, N=11 and h=0.06; 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 sizesh 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 F(ρ). 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 I(ρ) 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 vstat in (1a), (1b), (1c), and (1d), the conservation law reduces to(41)∂tρ+∇·ρIρ=0,Iρ=-∇η∗ρ1+∇η∗ρ22.In [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 Ω=[-1,1]2. The initial density ρ(x,0) is set to 1 for x∈[-1/2,1/2]2; otherwise ρ(x,0)=0. Additionally, we compute the simplified model (41) for three different mollifiers; that is, σ=25,100,400. The step sizes of the finite volume approach are Δx1=Δx2=0.01 and Δt=0.005.The results of the finite volume approach are shown in Figure8. 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 σ=25. 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 σ=100,400 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.Figure 8
Numerical solution of the simplified model (41) computed by the finite volume approach with dimensional splitting. Visualized for times t=0.1,0.2,0.3 and smoothing function parameter σ=25,100,400.
(a)
Timet=0.1, σ=25
(b)
Timet=0.2, σ=25
(c)
Timet=0.3, σ=25
(d)
Timet=0.1, σ=100
(e)
Timet=0.2, σ=100
(f)
Timet=0.3, σ=100
(g)
Timet=0.1, σ=400
(h)
Timet=0.2, σ=400
(i)
Timet=0.3, σ=400Figure9 shows the results of the discontinuous Galerkin method for different triangulations and polynomial degrees; however, the results are plotted for the time t=2 and σ=100. 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)).Figure 9
Comparison of the results of the discontinuous Galerkin method with different mesh-sizesh and polynomial degrees N. The plots show the solution at time t=0.2 for a smoothness parameter σ=100.
(a)
h
≈
0.4, N=1
(b)
h
≈
0.4, N=5
(c)
h
≈
0.4, N=11
(d)
h
≈
0.2, N=1
(e)
h
≈
0.2, N=3
(f)
h
≈
0.2, N=7
(g)
h
≈
0.1, N=1
(h)
h
≈
0.1, N=2
(i)
h
≈
0.1, N=4


## 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 fieldvstat(x) 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 γ=25 according to (3). The mollifier η, occurring in the operator I(ρ), is defined as(40)ηx=σ2πexp⁡-12σx22,σ=2500.In this example, the maximal density is set to ρmax⁡=1. The strength of the term I(ρ) is selected as ϵ=2vT, where the velocity of the conveyor belt is vT=0.395m/s. Furthermore, the total time horizon is T=7.


### 4.1.2. Finite Volume Settings

The grid sizes of the finite volume approach with dimensional splitting are chosen asΔx1=Δx2=5·10-3, Δt=1.25·10-3 in the following computation.


### 4.1.3. Discontinuous Galerkin Settings

The discontinuous Galerkin method uses a triangulationΩh with a maximal triangle edge length h=0.1. The polynomial degree of each finite element is N=10. ODE system (27) is solved by the explicit Euler method with a time step size Δt=10-3. Thereby the filter procedure is called in each computational step of the ODE solver. The filter settings are selected as β=36, s=6, and Nc=1.The results are shown in Figure6. 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 vT. A formation of congestion can be seen in all results.Figure 6
Results of the finite volume method with splitting:Δx1=Δx2=5⋅10-3 (left) and results of the discontinuous Galerkin method: h≈0.1, N=10 (right).
(a)
Timet=0.50
(b)
Timet=0.50
(c)
Timet=1.00
(d)
Timet=1.00
(e)
Timet=1.50
(f)
Timet=1.50
(g)
Timet=2.00
(h)
Timet=2.00In all plots, we recognize a weak dispersing of density (cf. Figures6(g) and 6(h)). This is caused by the term vdyn=H~(ρ-ρmax⁡)I(ρ). The smoothed modification H~(ρ-ρmax⁡) is never zero for ρ<ρmax⁡. Consequently, the dispersing term I(ρ) 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 I(ρ) 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 N=10. 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 sizesh=0.1, h=0.06, and h=0.04. 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.Figure 7
Results of the discontinuous Galerkin method with different triangulations (h=0.1,0.06,0.04) and polynomials degrees N. All plots show the solution at time t=1.
(a)
h
≈
0.1, N=5
(b)
h
≈
0.1, N=11
(c)
h
≈
0.06, N=3
(d)
h
≈
0.06, N=7
(e)
h
≈
0.04, N=1
(f)
h
≈
0.04, N=3However, a rough triangulation or a low polynomial degree causes bad approximations (cf. Figure7(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 Tables1 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.Table 1
Computation times of the discontinuous Galerkin method (simulation process) with different grid sizesh and polynomial degrees N. The time is measured in seconds.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
7.14 s
12.30 s
13.42 s
3
9.30 s
18.63 s
51.70 s
5
17.31 s
46.29 s
—
7
30.10 s
111.94 s
—
9
49.58 s
—
—
11
88.70 s
—
—Table 2
Computation times in seconds for the convolution preprocessing due to the grid sizeh and polynomial degree N.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
0.06 s
1.12 s
1.73 s
3
0.48 s
3.88 s
79.99 s
5
2.01 s
60.82 s
—
7
5.44 s
420.08 s
—
9
47.11 s
—
—
11
183.86 s
—
—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 mostNp·K coefficients. Furthermore, there are Np·K nodal points and the convolution is evaluated twice in each dimension. Thus, it is necessary to calculate and store about 2·(Np·K)2 coefficients. As a consequence, the computer was not able to run the preprocessing routine successfully for small h and a large N; for example, N=11 and h=0.06; 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 sizesh for computation of accurate approximations and evaluating the corresponding convergence behavior.


## 4.1.1. Macroscopic Model Settings

The fieldvstat(x) 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 γ=25 according to (3). The mollifier η, occurring in the operator I(ρ), is defined as(40)ηx=σ2πexp⁡-12σx22,σ=2500.In this example, the maximal density is set to ρmax⁡=1. The strength of the term I(ρ) is selected as ϵ=2vT, where the velocity of the conveyor belt is vT=0.395m/s. Furthermore, the total time horizon is T=7.


## 4.1.2. Finite Volume Settings

The grid sizes of the finite volume approach with dimensional splitting are chosen asΔx1=Δx2=5·10-3, Δt=1.25·10-3 in the following computation.


## 4.1.3. Discontinuous Galerkin Settings

The discontinuous Galerkin method uses a triangulationΩh with a maximal triangle edge length h=0.1. The polynomial degree of each finite element is N=10. ODE system (27) is solved by the explicit Euler method with a time step size Δt=10-3. Thereby the filter procedure is called in each computational step of the ODE solver. The filter settings are selected as β=36, s=6, and Nc=1.The results are shown in Figure6. 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 vT. A formation of congestion can be seen in all results.Figure 6
Results of the finite volume method with splitting:Δx1=Δx2=5⋅10-3 (left) and results of the discontinuous Galerkin method: h≈0.1, N=10 (right).
(a)
Timet=0.50
(b)
Timet=0.50
(c)
Timet=1.00
(d)
Timet=1.00
(e)
Timet=1.50
(f)
Timet=1.50
(g)
Timet=2.00
(h)
Timet=2.00In all plots, we recognize a weak dispersing of density (cf. Figures6(g) and 6(h)). This is caused by the term vdyn=H~(ρ-ρmax⁡)I(ρ). The smoothed modification H~(ρ-ρmax⁡) is never zero for ρ<ρmax⁡. Consequently, the dispersing term I(ρ) 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 I(ρ) 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 N=10. 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 sizesh=0.1, h=0.06, and h=0.04. 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.Figure 7
Results of the discontinuous Galerkin method with different triangulations (h=0.1,0.06,0.04) and polynomials degrees N. All plots show the solution at time t=1.
(a)
h
≈
0.1, N=5
(b)
h
≈
0.1, N=11
(c)
h
≈
0.06, N=3
(d)
h
≈
0.06, N=7
(e)
h
≈
0.04, N=1
(f)
h
≈
0.04, N=3However, a rough triangulation or a low polynomial degree causes bad approximations (cf. Figure7(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 Tables1 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.Table 1
Computation times of the discontinuous Galerkin method (simulation process) with different grid sizesh and polynomial degrees N. The time is measured in seconds.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
7.14 s
12.30 s
13.42 s
3
9.30 s
18.63 s
51.70 s
5
17.31 s
46.29 s
—
7
30.10 s
111.94 s
—
9
49.58 s
—
—
11
88.70 s
—
—Table 2
Computation times in seconds for the convolution preprocessing due to the grid sizeh and polynomial degree N.
N
h
=
0.1
h
=
0.06
h
=
0.04
1
0.06 s
1.12 s
1.73 s
3
0.48 s
3.88 s
79.99 s
5
2.01 s
60.82 s
—
7
5.44 s
420.08 s
—
9
47.11 s
—
—
11
183.86 s
—
—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 mostNp·K coefficients. Furthermore, there are Np·K nodal points and the convolution is evaluated twice in each dimension. Thus, it is necessary to calculate and store about 2·(Np·K)2 coefficients. As a consequence, the computer was not able to run the preprocessing routine successfully for small h and a large N; for example, N=11 and h=0.06; 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 sizesh 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 F(ρ). 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 I(ρ) 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 vstat in (1a), (1b), (1c), and (1d), the conservation law reduces to(41)∂tρ+∇·ρIρ=0,Iρ=-∇η∗ρ1+∇η∗ρ22.In [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 Ω=[-1,1]2. The initial density ρ(x,0) is set to 1 for x∈[-1/2,1/2]2; otherwise ρ(x,0)=0. Additionally, we compute the simplified model (41) for three different mollifiers; that is, σ=25,100,400. The step sizes of the finite volume approach are Δx1=Δx2=0.01 and Δt=0.005.The results of the finite volume approach are shown in Figure8. 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 σ=25. 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 σ=100,400 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.Figure 8
Numerical solution of the simplified model (41) computed by the finite volume approach with dimensional splitting. Visualized for times t=0.1,0.2,0.3 and smoothing function parameter σ=25,100,400.
(a)
Timet=0.1, σ=25
(b)
Timet=0.2, σ=25
(c)
Timet=0.3, σ=25
(d)
Timet=0.1, σ=100
(e)
Timet=0.2, σ=100
(f)
Timet=0.3, σ=100
(g)
Timet=0.1, σ=400
(h)
Timet=0.2, σ=400
(i)
Timet=0.3, σ=400Figure9 shows the results of the discontinuous Galerkin method for different triangulations and polynomial degrees; however, the results are plotted for the time t=2 and σ=100. 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)).Figure 9
Comparison of the results of the discontinuous Galerkin method with different mesh-sizesh and polynomial degrees N. The plots show the solution at time t=0.2 for a smoothness parameter σ=100.
(a)
h
≈
0.4, N=1
(b)
h
≈
0.4, N=5
(c)
h
≈
0.4, N=11
(d)
h
≈
0.2, N=1
(e)
h
≈
0.2, N=3
(f)
h
≈
0.2, N=7
(g)
h
≈
0.1, N=1
(h)
h
≈
0.1, N=2
(i)
h
≈
0.1, N=4


## 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.



---

*Source: 341893-2015-11-08.xml*