Loading [MathJax]/extensions/TeX/boldsymbol.js
Research article

A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows

  • Received: 29 August 2022 Revised: 05 November 2022 Accepted: 10 November 2022 Published: 28 November 2022
  • A triangular system of conservation laws with discontinuous flux that models the one-dimensional flow of two disperse phases through a continuous one is formulated. The triangularity arises from the distinction between a primary and a secondary disperse phase, where the movement of the primary disperse phase does not depend on the local volume fraction of the secondary one. A particular application is the movement of aggregate bubbles and solid particles in flotation columns under feed and discharge operations. This model is formulated under the assumption of a variable cross-sectional area. A monotone numerical scheme to approximate solutions to this model is presented. The scheme is supported by three partial theoretical arguments. Firstly, it is proved that it satisfies an invariant-region property, i.e., the approximate volume fractions of the three phases, and their sum, stay between zero and one. Secondly, under the assumption of flow in a column with constant cross-sectional area it is shown that the scheme for the primary disperse phase converges to a suitably defined entropy solution. Thirdly, under the additional assumption of absence of flux discontinuities it is further demonstrated, by invoking arguments of compensated compactness, that the scheme for the secondary disperse phase converges to a weak solution of the corresponding conservation law. Numerical examples along with estimations of numerical error and convergence rates are presented for counter-current and co-current flows of the two disperse phases.

    Citation: Raimund Bürger, Stefan Diehl, M. Carmen Martí, Yolanda Vásquez. A difference scheme for a triangular system of conservation laws with discontinuous flux modeling three-phase flows[J]. Networks and Heterogeneous Media, 2023, 18(1): 140-190. doi: 10.3934/nhm.2023006

    Related Papers:

    [1] Raimund Bürger, Kenneth H. Karlsen, John D. Towers . On some difference schemes and entropy conditions for a class of multi-species kinematic flow models with discontinuous flux. Networks and Heterogeneous Media, 2010, 5(3): 461-485. doi: 10.3934/nhm.2010.5.461
    [2] Raimund Bürger, Christophe Chalons, Rafael Ordoñez, Luis Miguel Villada . A multiclass Lighthill-Whitham-Richards traffic model with a discontinuous velocity function. Networks and Heterogeneous Media, 2021, 16(2): 187-219. doi: 10.3934/nhm.2021004
    [3] Raimund Bürger, Stefan Diehl, María Carmen Martí . A conservation law with multiply discontinuous flux modelling a flotation column. Networks and Heterogeneous Media, 2018, 13(2): 339-371. doi: 10.3934/nhm.2018015
    [4] Mauro Garavello, Roberto Natalini, Benedetto Piccoli, Andrea Terracina . Conservation laws with discontinuous flux. Networks and Heterogeneous Media, 2007, 2(1): 159-179. doi: 10.3934/nhm.2007.2.159
    [5] Giuseppe Maria Coclite, Lorenzo di Ruvo, Jan Ernest, Siddhartha Mishra . Convergence of vanishing capillarity approximations for scalar conservation laws with discontinuous fluxes. Networks and Heterogeneous Media, 2013, 8(4): 969-984. doi: 10.3934/nhm.2013.8.969
    [6] Maya Briani, Emiliano Cristiani . An easy-to-use algorithm for simulating traffic flow on networks: Theoretical study. Networks and Heterogeneous Media, 2014, 9(3): 519-552. doi: 10.3934/nhm.2014.9.519
    [7] Christophe Chalons, Paola Goatin, Nicolas Seguin . General constrained conservation laws. Application to pedestrian flow modeling. Networks and Heterogeneous Media, 2013, 8(2): 433-463. doi: 10.3934/nhm.2013.8.433
    [8] Clément Cancès . On the effects of discontinuous capillarities for immiscible two-phase flows in porous media made of several rock-types. Networks and Heterogeneous Media, 2010, 5(3): 635-647. doi: 10.3934/nhm.2010.5.635
    [9] Wen Shen . Traveling wave profiles for a Follow-the-Leader model for traffic flow with rough road condition. Networks and Heterogeneous Media, 2018, 13(3): 449-478. doi: 10.3934/nhm.2018020
    [10] Adriano Festa, Simone Göttlich, Marion Pfirsching . A model for a network of conveyor belts with discontinuous speed and capacity. Networks and Heterogeneous Media, 2019, 14(2): 389-410. doi: 10.3934/nhm.2019016
  • A triangular system of conservation laws with discontinuous flux that models the one-dimensional flow of two disperse phases through a continuous one is formulated. The triangularity arises from the distinction between a primary and a secondary disperse phase, where the movement of the primary disperse phase does not depend on the local volume fraction of the secondary one. A particular application is the movement of aggregate bubbles and solid particles in flotation columns under feed and discharge operations. This model is formulated under the assumption of a variable cross-sectional area. A monotone numerical scheme to approximate solutions to this model is presented. The scheme is supported by three partial theoretical arguments. Firstly, it is proved that it satisfies an invariant-region property, i.e., the approximate volume fractions of the three phases, and their sum, stay between zero and one. Secondly, under the assumption of flow in a column with constant cross-sectional area it is shown that the scheme for the primary disperse phase converges to a suitably defined entropy solution. Thirdly, under the additional assumption of absence of flux discontinuities it is further demonstrated, by invoking arguments of compensated compactness, that the scheme for the secondary disperse phase converges to a weak solution of the corresponding conservation law. Numerical examples along with estimations of numerical error and convergence rates are presented for counter-current and co-current flows of the two disperse phases.



    It is the purpose of this work to introduce, and in part analyze, a numerical scheme for a system of conservation laws with source terms of the type

    t(A(z)(ϕψ))+z(A(z)(J(ϕ,z,t)˜F(ϕ,ψ,z,t)))=Kk=1QF,k(t)(ϕF,k(t)ψF,k(t))δ(zzF,k), (1.1)

    where t is time, z is spatial position, and ϕ and ψ are the volume fractions of the primary and secondary disperse phases, respectively. Both disperse phases move within the continuous phase of the one-dimensional flow. We let A(z) denote a variable cross-sectional area. The flux functions J and ˜F are discontinuous across the positions z=zU<zF,1<<zF,K<zE, and due to constitutive assumptions of the model, are nonlinear functions of ϕ and ψ. The right-hand side of Eq (1.1) describes singular sources located at z=zF,k, k=1,,K, and is composed of given functions. It is assumed that QF,k(t) is the volumetric bulk flow of the mixture (of the continuous and two disperse phases) injected at z=zF,k, and that ϕF,k(t) and ψF,k(t) are the volume fractions of the primary and secondary disperse phases in the feed flow, respectively. The system Eq (1.1) is posed on ΠT:=R×(0,T) together with initial conditions

    ϕ(z,0)=ϕ0(z)for all zR, (1.2a)
    ψ(z,0)=ψ0(z)for all zR, (1.2b)

    where we assume that

    0ϕ0(z)1,0ψ0(z)1ϕ0(z)for allzR (1.3)

    along with

    TV(ϕ0)<,TV(ψ0)<. (1.4)

    Likewise, we assume that ϕF,k and ψF,k are piecewise continuous functions of bounded variation with a finite number of discontinuities and that they satisfy the bounds

    0ϕF,k(t)1,0ψF,k(t)1ϕF,k(t)for allk=1,,Kandt[0,T]. (1.5)

    (In later parts of the analysis we will assume that these functions and the bulk flows QF,k are constants.) If θ denotes the volume fraction of the continuous phase, then we assume that

    0ϕ,ψ,θ1;ϕ+ψ+θ=1, (1.6)

    which motivates assumptions Eq (1.3) and Eq (1.5). (Of course, satisfaction of Eq (1.6) by exact or numerical solutions of Eq (1.1), Eq (1.2) on ΠT needs to be proved.)

    A specific application that gives rise to the system Eq (1.1) is a model of a flotation column [8,9], where ϕ denotes the volume fraction of bubbles and ψ that of solid particles (Figure 1). The bottom of the column has the coordinate zU (the underflow) and the top zE (the effluent). The primary disperse phase of bubbles – specifically, aggregate bubbles, to which hydrophobic valuable particles (minerals) are attached – is assumed to flow through the suspension of solid particles and liquid independently of the volume fraction of solids. The secondary disperse phase consists of solid hydrophilic particles (ore) that move in the remaining space outside the bubbles. If the solid particles of the secondary disperse phase have a density larger than that of the fluid, the two disperse phases undergo counter-current, and otherwise, co-current flow. The distinction between primary and secondary disperse phase also becomes evident in the flux functions: the flux J of the primary disperse phase depends on ϕ only (besides z and t), while that of the secondary disperse phase, ˜F, depends both on ϕ and ψ. Thus, the system Eq (1.1) is triangular; however, it is generally non-strictly hyperbolic; see [9], where a counter-current model of the form Eq (1.1) is studied.

    Figure 1.  Schematic of a one-dimensional column with K=3 inlets and K+1=4 zones, where QU is the downwards volumetric outflow, QF,j is the volumetric flow at the inlet zF,j, for each j=1,,K, and QE is the upwards volumetric outflow. Note that the distances between the inlets/outlets are arbitrary and the cross-sectional area A=A(z) may vary piecewise continuously (although the figure shows a piecewise constant example).

    The main contribution of this work is an easily implemented explicit monotone numerical scheme for Eq (1.1). To properly address the theoretical support we are able to provide for this scheme, we refer to the complete model Eq (1.1), with all assumptions stated so far in effect, as "Model 1", and the corresponding scheme (that handles Model 1) as "Scheme 1". Additional properties of Scheme 1 can be established for two successively simplified versions of Model 1, named "Model 2" and "Model 3", for which the corresponding versions of Scheme 1 are addressed as "Scheme 2" and "Scheme 3", respectively.

    The scheme is supported by three partial theoretical arguments. Firstly, it is proved that Scheme 1 satisfies an invariant-region property, i.e., the approximate volume fractions satisfy a discrete analogue of Eq (1.6) at every point. Secondly, the assumption of a constant cross-sectional area, i.e.,

    Aconstant,A>0, (1.7)

    and time-independent feed and volume rates defines Model 2, and it is shown that the corresponding scheme for the primary disperse phase (the "ϕ-scheme" of Scheme 2; in short, "ϕ-Scheme 2") converges to a suitably defined entropy solution. Thirdly, we additionally assume that there are no flux discontinuities, so that Model 2 is reduced to the triangular system of conservation laws ("Model 3")

    tϕ+zJ(ϕ)=0, (1.8a)
    tψ+z˜F(ϕ,ψ)=0,(z,t)ΠT, (1.8b)

    where J and ˜F are z- and t-independent versions of the fluxes arising in Eq (1.1) and Eq (1.8) is equipped with the initial conditions Eq (1.2), where assumptions Eq (1.3) remain in effect. The corresponding reduced version of Scheme 2 that handles Model 3 is called "Scheme 3." Under these additional assumptions, we may invoke arguments of compensated compactness to prove that the scheme for the secondary disperse phase (the "ψ-Scheme 3") converges to a weak solution of the corresponding conservation law Eq (1.8b). Summarizing all arguments, we prove that Scheme 3 converges to a weak solution of the system Eq (1.8) in the sense of the following definition.

    Definition 1.1. The pair (ϕ,ψ) is called a weak solution of the initial value problem Eq (1.8), Eq (1.2) if

    (ⅰ) The functions ϕ and ψ belong to L(ΠT).

    (ⅱ) The functions ϕ and ψ satisfy Eq (1.8), Eq (1.2) in the sense of distributions on ΠT, that is, for each smooth test function ζ with compact support in ΠT, the following identities hold:

    ΠT(ϕtζ+J(ϕ)zζ)dzdt+Rϕ0(z)dz=0, (1.9)
    ΠT(ψtζ+˜F(ϕ,ψ)zζ)dzdt+Rψ0(z)dz=0. (1.10)

    (ⅲ) The function ϕ is an entropy solution of the single conservation law Eq (1.8a), that is, for each smooth and nonnegative test function ζ with compact support in ΠT, the following inequality holds for all cR:

    ΠT(|ϕc|tζ+sgn(ϕc)(J(ϕ)J(k))zζ)dzdt+R|ϕ0(z)c|dz0. (1.11)

    Numerical experiments illustrate that Scheme 1 for the full model Eq (1.1) (Model 1), approximates expected solution behaviour for counter-current and co-current flows and that approximate numerical errors tend to zero as the mesh is refined.

    The system Eq (1.1) models the evolution of the primary unknown ϕ independently of the secondary unknown ψ. One application of such triangular systems is the process of column flotation, which is a solid-liquid separation process used in mineral processing, environmental and chemical engineering [10,11,27,28,42,45]. The model Eq (1.1) restricted to three-phase counter-current flow in a flotation column was originally proposed in [9]. Its nonlinear constitutive assumptions come from the drift-flux theory (used to analyze the bubbly and froth regions [41,50,51]) and the solids-flux theory (for particles in a liquid [24,25,38]). In [9], the construction of steady-state solutions is detailed, where conservation laws with discontinuous flux are a key ingredient with a specific entropy condition [2,22,29]. The most interesting desired steady states are classified in [9] and visualized in graphical so-called operating charts that show how the control variables QU, QF:=QF,1 (feed mixture of gas, solids and water) and QW:=QF2 (feed washwater) should be chosen.

    The mathematical and numerical difficulties associated with Eq (1.1) are twofold; namely, one has to deal with discontinuities of the fluxes with respect to z, as well as with the definition of the governing model by a (triangular) system of conservation laws (in contrast to otherwise similar, known scalar two-phase models arising in flotation or sedimentation [8,13,22,23]). The well-known difficulty of conservation laws with discontinuous flux lies in the appropriate formulation of admissibility conditions of jumps of the solution across discontinuities of the flux such that the resulting concept of weak (discontinuous) solutions supplied with an entropy condition would admit a uniqueness result. There exist many criteria for selecting unique solutions (e.g., [1,22]), each of which corresponds to a particular physical reality and relies on specific assumptions on the fluxes adjacent to a discontinuity. A unified treatment of this problem is advanced in [2]. While scalar conservation laws with discontinuous flux have been studied widely, only a few analyses of systems with discontinuous flux are available (e.g., [16,47]). That said, its triangular nature makes Eq (1.1) potentially easier to treat than a full 2×2 system of conservation laws (where the flux of each component would depend on both unknowns).

    The triangular system with discontinuous flux studied in [9] was solved numerically with a staggered-grid scheme that utilizes the triangular structure of Eq (1.1). Such a semi-Godunov scheme for general triangular hyperbolic systems is one of the two suggested schemes by Karlsen et al. [32,33], who proved convergence of the numerical solutions under certain assumptions on the flux functions. We here propose a simpler numerical scheme (on a single grid) that is easier to implement and analyze. The analysis (of the scheme proposed under simplifying assumptions) relies on the aligned version of the scheme introduced in [33] and in particular on the convergence analysis of an Engquist-Osher scheme for multi-dimensional triangular system of conservation laws by Coclite et al. [18]. (The proof of convergence done in [18] is motivated by the more easily proven convergence of a vanishing viscosity approximation for the same model, see [17].) These analyses, and the present treatment for the reduced model Eq (1.8), rely on compactness techniques that use discrete entropy inequalities and the compensated compactness framework.

    Further applications and results on the analysis of triangular systems include two-component chromatography [3]. Furthermore, polymer flooding in oil recovery is modelled by a 2×2 system [31], which can be converted to a triangular system in Lagrange coordinates [49]. In [20,40,48], the authors study the delta shock wave formation in solutions of triangular system of conservation laws from the so-called generalized pressureless gas dynamics model. Bressan et al. [5] established the existence and uniqueness of vanishing viscosity solutions for scalar conservation laws for a Cauchy problem and their results can be applied to a triangular system under suitable assumptions. The results of Karlsen et al. [32,33] for general triangular systems can be applied to models of three-phase flows in porous media, for example, in oil recovery.

    The remainder of this paper is organized as follows. In Section 2, the model of [9] of gas-solid-liquid three-phase flow in a flotation column from is written in a slightly more general form. Starting from the balance equations of the three phases we outline the derivation of the algebraic forms of the fluxes J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) arising in the governing PDE system Eq (1.1). In Section 3 the numerical method (Scheme 1) proposed for the approximation of solutions to the initial value problem Eq (1.1), Eq (1.2) (Model 1) is detailed, where computational effort is essentially reduced to the interior of the vessel (cf. Figure 1). After outlining the discretization in Section 3.1, we specify the numerical fluxes and update formulas for the primary and secondary disperse phases in Sections 3.2 and 3.3, respectively. Both formulas are adapted to the particular algebraic form of the fluxes J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) and involve upwind discretizations, a particular monotone discretization for "concentration times velocity" fluxes from [6], and the Engquist-Osher numerical flux [26]. We then prove in Section 3.4 that Scheme 1 is monotone and that the numerical solutions satisfy a so-called invariant-region property (Theorems 3.1 and 3.2), that is, a discrete analogue of Eq (1.6), provided, of course, that the initial data satisfy Eq (1.3) and the time step and spatial meshwidth satisfy a CFL condition. Section 4 provides further partial results of the convergence analysis of the numerical scheme based on additional simplifying assumptions, namely those of a constant cross-sectional area A and constant bulk and feed flows QU, QF,k, ϕF,k and ψF,k (k=1,,K). Thus, the scheme under discussion is Scheme 2. We can then prove convergence of the ϕ-Scheme 2 (the one that discretizes the ϕ-component of the governing PDE; Section 4.1) and L1 Lipschitz continuity of the ψ-Scheme 2 (Section 4.2). If in addition all z-dependent flux discontinuities are removed, so that the governing PDE system is Eq (1.8) (Model 3) and the scheme reduces to Scheme 3, we may apply compensated compactness techniques to prove convergence of the ψ-scheme (Section 5). For the simplified problems, the initial conditions Eq (1.2) and assumptions Eq (1.3) and Eq (1.5) are imposed, so Theorems 3.1 and 3.2 remain in effect. While in that case the convergence of (the monotone) ϕ-Scheme 3 to an entropy solution of Eq (1.8a) follows by standard arguments (for monotone schemes), the principal result of Section 5 is convergence of ψ-Scheme 3 to a weak solution of Eq (1.8b) (Lemma 5.5 and Theorem 5.1). Estimations of errors and convergence order of the numerical method can be found in Section 6.1. Some numerical examples are presented in Section 6, starting with preliminaries (Section 6.1). First, in Section 6.3, we use a smooth solution to estimate the order of convergence. Later on, we present two numerical examples that illustrate the model predictions for counter-current (Section 6.4) and co-current flows (Sections 6.5 and 6.6). Finally, some conclusions are drawn in Section 7, and Appendix A contains the proofs of Lemma 5.5 and Theorem 5.1.

    The density of each phase is assumed constant, so the conservation of mass can be expressed by the balance equations (vϕ, vψ, and vθ are the phase velocities)

    t(A(z)ϕ)+z(A(z)ϕvϕ)=Kk=1QF,k(t)ϕF,k(t)δ(zzF,k), (2.1)
    t(A(z)ψ)+z(A(z)ψvψ)=Kk=1QF,k(t)ψF,k(t)δ(zzF,k), (2.2)
    t(A(z)θ)+z(A(z)θvθ)=Kk=1QF,k(t)(1ϕF,k(t)ψF,k(t))δ(zzF,k), (2.3)

    where the right-hand sides contain Dirac symbols, the feed volume fractions ϕF,k and ψF,k of the disperse phases, and the corresponding volume fraction 1ϕF,k(t)ψF,k(t) of the continuous phase, at the inlet located at z=zF,k, k=1,,K.

    We define the volume-average velocity, or bulk velocity, of the mixture by

    q:=ϕvϕ+ψvψ+θvθ,

    and replace Eq (2.3) by the sum of Eqs (2.1)–(2.3), which is

    z(A(z)q)=Kk=1QF,k(t)δ(zzF,k), (2.4)

    hence q varies with z due to the K inlets and the variable cross-sectional area. We define Q(z,t):=A(z)q(z,t) and integrate Eq (2.4) from any point z0<zU to obtain

    Q(z,t)=Q(z0,t)+Kk=1QF,k(t)H(zzF,k),

    where H() is the Heaviside function. If the volumetric underflow QU(t) is given, then Q(z,t)=QU(t) for z<zU, and

    Q(z,t)=QU(t)+Kk=1QF,k(t)H(zzF,k)=QU(t)+k:zF,kzQF,k(t).

    This continuity equation of the mixture replaces Eq (2.3). Next, Eq (2.1) and Eq (2.2) are rewritten in terms of q and two constitutive functions. We refer to the continuous phase and the secondary disperse phase as "secondary mixture", and define the volume fraction of the secondary disperse phase within the secondary mixture as

    φ:=ψψ+θ=ψ1ϕ(when ϕ<1),

    where 0φ1 by Eq (1.6). The volume-average velocity of the secondary mixture is

    qs:=ψvψ+θvθψ+θ=φvψ+1ϕψ1ϕvθ=φvψ+(1φ)vθ.

    It is then assumed that within [zU,zE), the relative velocity vϕs:=vϕqs of the primary disperse phase with respect to the secondary mixture is a given constitutive function ˜vϕs(ϕ), while outside that interval, both phases move at the same velocity, so their velocity difference is zero. Thus, in terms of the characteristic function

    γ(z):=χ[zU,zE)(z):={1for z[zU,zE),0for z[zU,zE),

    this assumption can be expressed as vϕs=γ(z)˜vϕs(ϕ). Within [zU,zE), the relative velocity of the secondary disperse phase with respect to the continuous phase vψθ:=vψvθ is supposed to be a given function ˜vψθ of φ, that is, vψθ=γ(z)˜vψθ(φ).

    The definitions of all velocities imply the identities

    ϕvϕ=ϕq+γ(z)ϕ(1ϕ)˜vϕs(ϕ),ψvψ=ψq+γ(z)ψ((1φ)˜vψθ(φ)ϕ˜vϕs(ϕ)) (2.5)

    for the (unweighted) fluxes ϕvϕ and ψvψ arising in Eq (2.1) and Eq (2.2), respectively. It is then useful to introduce the velocity and flux functions

    W(ϕ):=(1ϕ)˜vϕs(ϕ),V(φ):=σ(1φ)˜vψθ(φ),j(ϕ):=ϕW(ϕ),f(φ):=φV(φ), (2.6)

    where σ=±1 is chosen depending on the application such that V(φ),f(φ)0 (for standard convenience, e.g., when plotting their graphs); σ=1 for co-current flows (upwards) and σ=1 for counter-current flows. The velocity and flux of the secondary disperse phase with respect to z are therefore σV(φ) and σf(φ), respectively. We assume that W,V0 and V(1)=W(1)=0, as well as that

    fhas one local maximumωand one inflection point˜ω,0<ω<˜ω<1. (2.7)

    Combining Eq (2.5) and Eq (2.6) we obtain the expressions

    ϕvϕ=ϕq+γ(z)ϕW(ϕ)=:J(ϕ,z,t),ψvψ=(1ϕ)φq+γ(z)((1ϕ)φσV(φ)φϕW(ϕ))=:F(ϕ,φ,z,t) (2.8)

    for the total fluxes of Eq (2.1) and Eq (2.2). For ϕ<1, we define the final flux function

    ˜F(ϕ,ψ,z,t):=F(ϕ,ψ1ϕ,z,t)=ψq+γ(z)(ψσV(ψ1ϕ)ψϕW(ϕ)1ϕ), (2.9)

    whereas for ϕ=1, we set ˜F(1,ψ,z,t):=0 (since F(1,φ,z,t)=0 for all φ[0,1]). Substituting Eq (2.8) and Eq (2.9) into Eq (2.1) and Eq (2.2), respectively, we obtain the final governing PDE system Eq (1.1).

    Illustrations and numerical examples are based on the expressions

    W(ϕ)=vterm,p(1ϕ)npfor 0ϕ1,np>1, (2.10)
    V(φ)=vterm,s(1φ)nsfor 0φ1,ns>1 (2.11)

    (see [46]), where vterm,p and vterm,s are the terminal velocities of a single particle of the primary and secondary disperse phases, respectively, in an unbounded fluid. We set np=3.2, vterm,p=2.7cm/s, ns=2.5, and vterm,s=0.5cm/s along with σ=1. These values are used in Applications 1 and 2 in Section 6. The resulting nonlinearities of J(ϕ,z,t) and ˜F(ϕ,ψ,z,t) in the different zones of the column are illustrated in [52,Figure 3.2].

    The discretization of the model is based on the triangularity of the system of conservation laws Eq (1.1). The numerical fluxes for ϕ are based on the particular treatment of conservation laws with fluxes having an explicit "concentration times velocity" structure [6]. In each time step, an approximate solution ϕ of the first PDE of Eq (1.1) is obtained and used as a given piecewise constant function in space and time in the second PDE of Eq (1.1), which is updated accordingly.

    We define a computational domain [0,zend) (to be used for the error calculation; see Section 6.1) consisting of N cells by covering the vessel with N2 cells and placing one cell each below and above; see Figure 2 (left). This setup, with a finite spatial domain, is introduced for practical reasons and is the minimal spatial domain that captures the interior of the tank and the concentrations in the underflow and effluent zones. The formulation of the scheme (i.e., Scheme 1) and subsequent proof of invariant region property are referred to this computational domain, but for the convergence analysis the model is specified as the initial value problem Eq (1.1), Eq (1.2) with the initial data posed on the real line. This distinction is merely a formal one since on (,0) and (zend,) the model reduces to linear advection equations describing that matter is transported away from the unit at constant velocity (if no changes in A in these zones arise).

    Figure 2.  (Left) Discretization of ϕ and ψ in the application to flotation, where the height of the vessel is H=zEzU, there are K inlets, and the cross-sectional area A(z) has two values separated by a discontinuity at z=zF,2; cf. the examples in Sections 6.4 and 6.5, (right) enlarged view illustrating cell division for error computations when Δrz is the discretization of the reference solution (see Section 6.1).

    Given the column height H, we define Δz:=H/(N2), the cell boundaries zi:=iΔz, i=0,1,,N, and the cells (intervals) Ii1/2:=[zi1,zi) and Ii:=[zi1/2,zi+1/2). We place the column between zU:=Δz=z1 and zE:=zU+H=(N1)Δz=zN1. Then the length of the interval of error calculation is zend:=H+2Δz=NΔz. Each injection point zF,k is assumed to belong to one cell Ii1/2 and we define the dimensionless quantity

    δk,i1/2:=Ii1/2δzF,k(z)dz:={1if zF,kIi1/2,0otherwise. (3.1)

    The cross-sectional area A=A(z) is allowed to have a finite number of discontinuities and it is discretized by

    Ai:=1ΔzIiA(z)dz,Ai+1/2:=1ΔzIi+1/2A(z)dz.

    We simulate NT time steps up to the final time T:=NTΔt, with the fixed time step Δt satisfying the Courant-Friedrichs-Lewy (CFL) condition

    Δt(2Q,TAmin+M(max{V(0),V}+W+W))Δz, (3.2)

    where

    M:=maxi=1,2,,N{Ai1Ai1/2,AiAi1/2},Amin:=mink=0,1/2,1,3/2,,NAk,Q,T:=max0tTKk=1QF,k(t),W:=max0ϕ1|W(ϕ)|.

    Finally, we set tn:=nΔt for n=0,1,,NT.

    The time-dependent feed functions are discretized as

    QnF,k:=1Δttn+1tnQF,k(t)dt,ϕnF,k:=1Δttn+1tnϕF,k(t)dt,

    for k=1,,K, and the same is made for ψF,k.

    The first equation of Eq (1.1) is discretized by combining upwind discretizations of qϕ with the particular scheme proposed in [6] for models with a "concentration times velocity" flux, as is the case for the term ϕW(ϕ).

    The initial data are discretized by

    ϕ0i1/2:=1Ai1/2ΔzIi1/2ϕ(z,0)A(z)dz.

    To advance from tn to tn+1 from given values ϕni1/2, i=1,,N, we define the numerical flux at z=zi by

    Jni:={ϕn1/2qn0for i=0,ϕni1/2qn+i+ϕni+1/2qni+γiϕni1/2W(ϕni+1/2)for i=1,,N1,ϕnN1/2qn+Nfor i=N, (3.3)

    where the notation

    a+:=max{a,0},a:=min{a,0},γi:=γ(zi),andqn+i:=(q(zi,tn))+

    is used. Since the bulk fluxes above and below the tank are directed away from it,

    ϕn1/2qn+0=0andϕnN+1/2qnN=0for any values of ϕn1/2andϕnN+1/2.

    To simplify the presentation, we use the middle line of Eq (3.3) as the definition of Jni together with ϕn1/2:=0 and ϕnN+1/2:=0. With the notation λ:=Δt/Δz and Qn+i:=Aiqn+i etc., the conservation law on Ii1/2 implies the update formula

    ϕn+1i1/2=ϕni1/2+λAi1/2(Ai1Jni1AiJni+Kk=1QnF,kϕnF,kδk,i1/2)=:Hi1/2(ϕni3/2,ϕni1/2,ϕni+1/2),i=1,,N. (3.4)

    Then we define the piecewise constant approximate solution ϕΔz on R×[0,T) by

    ϕΔz(z,t):=i,nχIi1/2(z)χ[tn,tn+1)(t)ϕni1/2, (3.5)

    where χΩ denotes the characteristic function of the set Ω.

    We discretize the initial data by

    ψ0i1/2:=1Ai1/2ΔzIi1/2ψ(z,0)A(z)dz.

    The well-known Engquist-Osher numerical flux [26] for a given continuous, piecewise differentiable flux function g and real values a and b on the left/right is given by

    G(g;a,b):=g(0)+a0max{0,g(s)}ds+b0min{0,g(s)}ds. (3.6)

    Then a consistent numerical flux corresponding to Eq (2.9) is

    Fni:=ψni1/2qn+i+ψni+1/2qni+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/21ϕni+1/2W(ϕni+1/2)),i=0,,N,

    where we set ψn1/2:=0 and ψnN+1/2:=0 with the same motivation as for ϕ above (these values are irrelevant). Here

    Gni(ψni1/2,ψni+1/2):=G(σfni;ψni1/2,ψni+1/2) (3.7)

    is the Engquist-Osher numerical flux associated with the function

    σfni(ψ):=σψ˜V(ψψnmax,i),˜V(u):={V(u)for u<1,0for u1, (3.8)

    where (ab:=min{a,b}, ab:=max{a,b})

    ψnmax,i:=(1ϕni1/2)(1ϕni+1/2)=1(ϕni1/2ϕni+1/2).

    The function ψσfni(ψ) is unimodal. Let ˆψni denote the maximum point of fni. For a given function ˜V the values ˆψni and ψnmax,i are related by the following lemma.

    Lemma 3.1. Assume that 0<ω<˜ω<1 are the unique local maximum and inflection point, respectively, of f(φ)=φV(φ) (cf. Eq (2.7)). Then ˆψni=ωψnmax,i for all i and n and all possible values 0ψnmax,i1. Moreover, the unique inflection point ψninfl,i(ˆψni,ψnmax,i) of fni satisfies ψninfl,i=˜ωψnmax,i for all i and n and all possible values 0ψnmax,i1. (See Figure 3.)

    Figure 3.  Illustration of Lemma 3.1.

    Proof. Assume that 0<ψnmax,i1. Since ˆψni is the unique solution ˆψni<ψnmax,i of

    ddψ(ψ˜V(ψψnmax,i))=0˜V(ψψnmax,i)+ψψnmax,i˜V(ψψnmax,i)=0,

    it follows that ω is the unique solution in (0,1) of ˜V(ω)+ω˜V(ω)=0 (cf. Eq (2.7)). By a similar argument, ˜ω is the unique solution of 2˜V(˜ω)+˜ω˜V(˜ω)=0.

    The Engquist-Osher numerical flux Eq (3.7), which in this form appears in Scheme 1 as well as in its reduced versions, Schemes 2 and 3, can now be computed as follows, where we recall that fni(0)=0. For σ=1 we get

    ψni1/20max{0,(fni)(s)}ds={fni(ψni1/2)if ψni1/2ˆψni,fni(ˆψni)if ψni1/2>ˆψni,ψni+1/20min{0,(fni)(s)}ds={0if ψni+1/2ˆψni,fni(ψni+1/2)fni(ˆψni)if ψni+1/2>ˆψni, (3.9)

    hence

    G(fni;ψni1/2,ψni+1/2)={fni(ψni1/2)if ψni1/2ˆψniandψni+1/2ˆψni,fni(ψni1/2)+fni(ψni+1/2)fni(ˆψni)if ψni1/2ˆψniandψni+1/2>ˆψni,fni(ˆψni)if ψni1/2>ˆψniandψni+1/2ˆψni,fni(ψni+1/2)if ψni1/2>ˆψniandψni+1/2>ˆψni. (3.10)

    By analogous reasoning we obtain for σ=1

    G(fni;ψni1/2,ψni+1/2)={fni(ψni+1/2)if ψni1/2ˆψniandψni+1/2ˆψni,fni(ˆψni)if ψni1/2ˆψniandψni+1/2>ˆψni,fni(ˆψni)fni(ψni1/2)fni(ψni+1/2)if ψni1/2>ˆψniandψni+1/2ˆψni,fni(ψni1/2)if ψni1/2>ˆψniandψni+1/2>ˆψni. (3.11)

    We define the difference operators Δai:=aiai1 and Δ+ai:=ai+1ai. Then the marching formula for ψ-Scheme 1 is

    ψn+1i1/2=ψni1/2+λAi1/2(Ai1Fni1AiFni+Kk=1QnF,kψnF,kδk,i1/2)=ψni1/2λAi1/2(Δ(ψni1/2Qn+i+ψni+1/2Qni+(Aγ)i(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/21ϕni+1/2W(ϕni+1/2)))Kk=1QnF,kψnF,kδk,i1/2),i=1,,N. (3.12)

    Then we define the piecewise constant approximate solution ψΔz on R×[0,T) by

    ψΔz(z,t):=i,nχIi1/2(z)χ[tn,tn+1)(t)ψni1/2. (3.13)

    We prove that Scheme 1, defined by the update formulas Eq (3.4) and Eq (3.12), is monotone, a property which then is used to prove the invariant-region property that the approximate solutions are positive and bounded.

    Theorem 3.1. If the CFL condition Eq (3.2) is satisfied, then the update formula for ϕ, Eq (3.4) (that is, ϕ-Scheme 1) is monotone and 0ϕni1/21 for i=1,,N and n=1,,NT.

    Proof. We recall the assumption Eq (1.3). We first prove monotonicity of the three-point scheme for ϕ Eq (3.4), i.e, that ϕn+1i1/2/ϕnk1/20 for all i=1,,N and k=i1,i,i+1. We have

    ϕn+1i1/2ϕni3/2=λAi1/2(Qn+i1+(Aγ)i1W(ϕni1/2))0,ϕn+1i1/2ϕni+1/2=λAi1/2(Qni(Aγ)iϕni1/2W(ϕni+1/2))0,ϕn+1i1/2ϕni1/2=1+λAi1/2(Qni1+(Aγ)i1ϕni3/2W(ϕni1/2)Qn+i(Aγ)iW(ϕni+1/2))1λ(2Q,TAmin+M(W+W))0,

    where we have used the CFL condition Eq (3.2).

    We now prove that if 0ϕni1/21 for all i, then 0ϕn+1i1/21 for all i. Clearly, Eq (1.3) implies that 0ϕ0i1/21 for all i. Since the scheme Eq (3.4) is monotone, Hi1/2 is non-decreasing in each argument. Since by assumption W(1)=0, we get the following estimation (where we use a++a=a):

    0λAi1/2Kk=1QnF,kϕnF,kδk,i1/2=Hi1/2(0,0,0)ϕn+1i1/2=Hi1/2(ϕni3/2,ϕni1/2,ϕni+1/2)Hi1/2(1,1,1)=1+λAi1/2((Qni1Qni)+Kk=1QnF,kϕnF,kδk,i1/2)1+λAi1/2(Kk=1(QnF,k)δk,i1/2+Kk=1QnF,kδk,i1/2)=1.

    Lemma 3.2. The function fni (cf. Eq (3.8)) satisfies (fni)max{V(0),V}.

    Proof. By Eq (2.7), the function f(φ)=φV(φ) has a single inflection point ˜ω(0,1) and by Lemma 3.1, fni has the inflection point ˜ωψnmax,i(0,ψnmax,i). We have (fni)(0)=V(0), (fni)(φ)=0 for ψnmax,iφ1 and the lowest (and negative) value of (fni) is obtained at its only critical point ˜ωψnmax,i, for which

    (fni)(˜ωψnmax,i)=˜V(˜ω)+˜ω˜V(˜ω)V.

    This concludes the proof.

    Lemma 3.3. There holds Gni(1ϕni1/2,1ϕni+1/2)=0 for all i and n.

    Proof. Assume that 0<ψnmax,i=(1ϕni1/2)(1ϕni+1/2)1. By Lemma 3.1, ˆψni<ψnmax,i, hence Eq (3.10), Eq (3.11), and

    ˜V((1ϕni1/2)/ψnmax,i)=˜V((1ϕni+1/2)/ψnmax,i)=0

    imply that

    Gni(1ϕni1/2,1ϕni+1/2)={fni(1ϕni+1/2)=0if σ=1,fni(1ϕni1/2)=0if σ=1.

    Theorem 3.2. Under the assumptions of Theorem 3.1, the update formula for ψ, Eq (3.12) (i.e., ψ-Scheme 1) is monotone and along with Eq (3.4) produces approximate solutions that satisfy 0ψni1/21ϕni1/2 for all i and n.

    Proof. Assumptions Eq (1.3) and Eq (1.5) imply that 0ψ0i1/21ϕ0i1/2 for all i and

    ψnF,k1ϕnF,kfor alln. (3.14)

    To prove that the scheme Eq (3.12) is monotone, we write it as

    ψn+1i1/2=Kni1/2(ψni3/2,ψni1/2,ψni+1/2) (3.15)

    and show that this expression is non-decreasing in each of its arguments.

    Since 0ϕni1/21 for a given n and all i, and appealing to Eq (3.6), we have

    ψn+1i1/2ψni3/2=λAi1/2(Qn+i1+(Aγ)i1Gni1ψni3/2)0,ψn+1i1/2ψni+1/2=λAi1/2(Qni(Aγ)iGniψni+1/2+(Aγ)iϕni1/21ϕni+1/2W(ϕni+1/2))0,ψn+1i1/2ψni1/2=1+λAi1/2(Qni1Qn+i+(Aγ)i1(Gni1ψni1/2ϕni3/2W(ϕni1/2)1ϕni1/2)(Aγ)iGniψni1/2)1λ(2Q,TAmin+M(Gniψni1/2Gni1ψni1/2+W(ϕni1/2)1ϕni1/2)).

    By Eq (3.6) and Lemma 3.2 we also obtain

    Gniψni1/2Gni1ψni1/2=(fni)(ψni1/2)+(fni)(ψni1/2)=|(fni)(ψni1/2)|(fni)max{V(0),V},

    and for the remaining term, we use that W(1)=0 to get

    W(ϕni1/2)1ϕni1/2=W(ϕni1/2)W(1)1ϕni1/2=W(ξ)Wfor some ξ(ϕni1/2,1).

    Hence, the CFL condition Eq (3.2) implies

    ψn+1i1/2ψni1/21λ(2Q,TAmin+M(max{V(0),V}+W))0.

    The inequalities proved imply that Kni1/2 is non-decreasing in each of its arguments. Now we use that 0ψni1/21ϕni1/2 for all i and Lemma 3.3 to obtain

    0λAi1/2Kk=1QnF,kψnF,kδk,i1/2=Hi1/2(0,0,0)ψn+1i1/2=Hi1/2(ψni3/2,ψni1/2,ψni+1/2)Hi1/2(1ϕni3/2,1ϕni1/2,1ϕni+1/2)=1ϕni1/2+λAi1/2(Ai1Fni1(1ϕni3/2,1ϕni1/2)AiFni(1ϕni1/2,1ϕni+1/2)+Kk=1QnF,kψnF,kδk,i1/2)=1ϕni1/2+λAi1/2((1ϕni3/2)Qn+i1+(1ϕni1/2)Qni1(Aγ)i1ϕni3/2W(ϕni1/2)(1ϕni1/2)Qn+i(1ϕni+1/2)Qni+(Aγ)iϕni1/2W(ϕni+1/2)+Kk=1QnF,kψnF,kδk,i1/2).

    Appealing to Eq (3.14) and the update formula for ϕ Eq (3.4), we get

    ψn+1i1/21ϕn+1i1/2+λAi1/2(Qn+i1+Qni1Qn+iQni+Kk=1QnF,kδk,i1/2)=1ϕn+1i1/2+λAi1/2{Qni1Qni+Kk=1QnF,kδk,i1/2}=1ϕn+1i1/2.

    The last equality holds since {}=0 irrespective of whether there is a source in the cell; Qni1Qni+QnF,k=0, or not; Qni1Qni=0.

    For ease of the argument, let us focus on the case of a constant interior cross-sectional area A, i.e., assume that Eq (1.7) is in effect. In addition, we assume that QnF,k, ϕnF,k, and ψnF,k (k=1,,K) are constant and therefore do not depend on n. The same is assumed for the underflow volumetric flow QU. (That is, we now study Scheme 2 suitable for Model 2.) Then Eq (3.4) and Eq (3.12) take the forms

    ϕn+1i1/2=ϕni1/2λΔ(ϕni1/2q+i+ϕni+1/2qi+γiϕni1/2W(ϕni+1/2))+λKk=1qF,kϕF,kδk,i1/2, (4.1)
    ψn+1i1/2=ψni1/2λΔ(ψni1/2q+i+ψni+1/2qi+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2ψni+1/2W(ϕni+1/2)1ϕni+1/2))+λKk=1qF,kψF,kδk,i1/2, (4.2)

    where qF,k:=QF,k/A. To embed the treatment into available analyses of schemes for conservation laws with discontinuous flux, we absorb the feed terms into the numerical flux. That is, we define ik:=i if δk,i1/2=1 (see Eq (3.1)). Then

    qi={qUif ii11,qU+qF,1++qF,lif iliil+11,l=1,,K1,qU+qF,1++qF,Kfor iiK. (4.3)

    Furthermore, we define the feed flux

    hF,i:={0if ii11,qF,1ϕF,1++qF,lϕF,lif iliil+11,l=1,,K1,qF,1ϕF,1++qF,KϕF,Kfor iiK, (4.4)

    such that

    hF,ihF,i1=Kk=1qF,kϕF,kδk,i1/2.

    Consequently, we may write the scheme (i.e., ϕ-Scheme 2) as

    ϕn+1i1/2=ϕni1/2λΔ(ϕni+1/2qi+ϕni1/2q+i+γiϕni1/2W(ϕni+1/2)+hF,i). (4.5)

    For later use we define the piecewise constant functions

    q(z):=qkandhF(z):=hF,kfor zF,k<z<zF,k+1,k=0,,K,

    where zF,0:=, zF,K+1:=, and we define the function

    h(z,v,u):=q(z)v+q+(z)u+γ(z)uW(v)+hF(z) (4.6)

    that allows us to write Eq (4.5) as

    ϕn+1i1/2=ϕni1/2λΔh(zi,ϕni+1/2,ϕni1/2). (4.7)

    The PDE for ϕ within Model 2, that is when the simplification Eq (1.7) is applied to Model 1, is the conservation law

    tϕ+zJ(ϕ,z)=0,(z,t)ΠT (4.8)

    with discontinuous flux

    J(ϕ,z)={q(z)ϕKk=1qF,kϕF,kfor z>zE,q(z)ϕKk=1qF,kϕF,k+j(ϕ)for zF,K<z<zE,q(z)ϕlk=1qF,kϕF,k+j(ϕ)for zF,l<z<zF,l+1,l=1,,K1,qUϕ+j(ϕ)for zU<z<zF,1,qUϕfor z<zU. (4.9)

    posed along with the initial condition Eq (1.2a).

    The choice of the appropriate solution concept for weak solutions, and the ways we may relate the model to the available theory of conservation laws with discontinuous flux, requires verifying whether J(ϕ,z) as given by Eq (4.9) satisfies the so-called "crossing condition" across each discontinuity

    zZ:={zU,zF,1,,zF,K,zE}. (4.10)

    Certain early well-posedness (existence, stability, and uniqueness) results for conservation laws with discontinuous flux (and related equations) rely on satisfaction of this condition [35], although later developments advance solution concepts that do not rely on satisfaction of the crossing condition [4,36,39]. In the present context this condition is satisfied for a particular discontinuity at z if the adjacent fluxes to the right and the left, J(ϕ,z+) and J(ϕ,z), satisfy

    ϕ1,ϕ2[0,1]:J(ϕ1,z+)J(ϕ1,z)<0<J(ϕ2,z+)J(ϕ2,z)ϕ1<ϕ2, (4.11)

    which means either the graphs of J(,z) and J(,z+) do not intersect, or if they do, there is at most one flux crossing ϕχ and the graph of J(,z) lies above that of J(,z+) to the left of ϕχ. For J(ϕ,z) as given by Eq (4.9) this condition is clearly satisfied for z{zE,zU} (considering that j(ϕ)>0 for 0<ϕ<1 implies that J(,z) and J(,z+) do not intersect in this case), while

    J(ϕ,z+F,l)J(ϕ,zF,l)=qF,l(ϕϕF,l)for l=1,,K.

    Thus, the crossing condition is satisfied also for z=zF,l, l=1,,K, since either ϕF,l=0 and the adjacent fluxes do not intersect in (0,1), or the intersection takes place at ϕχ=ϕF,l and Eq (4.11) holds since qF,l>0 for all l. The preceding consideration is analogous to the one for the simpler clarifier-thickener model (equivalent to K=1 in the present notation) studied e.g. in [13,14]. With the present analysis it is clear that the crossing condition is satisfied at each flux discontinuity zZ.

    Some of the available analyses refer to initial-value problems of the type

    tu+xF(u,x)=0for (x,t)ΠT,u(x,0)=u0(x)for xR,whereF(u,x):=H(x)g(u)+H(x)f(u) (4.12)

    where f and g are Lipschitz continuous functions of u denoting the "right" and "left" flux adjacent to a flux discontinuity across x=0 and H denotes the Heavyside function. The model problem Eq (4.12) features, of course, only one flux discontinuity (sitting at x=0), while Eq (4.9), Eq (1.2a) includes several of them at separate spatial locations. The study of Eq (4.12) is, however, sufficient for the analysis of each single flux discontinuity.

    Here we start from the concept of entropy solutions of type V introduced by Karlsen and Towers [36]. This concept does not appeal to the existence of traces of the unknown with respect to the interfaces zZ across which J(ϕ,z) is discontinuous. To state its adaptation to the situation at hand, we define the sets

    Π(K+3/2)T:=(zE,)×(0,T),Π(K+1/2)T:=(zF,K,zE)×(0,T),Π(k1/2)T:=(zF,k1,zF,k)×(0,T),k=2,,K,Π(1/2)T:=(zU,zF,1)×(0,T),Π(1/2)T:=(,zU)×(0,T).

    Definition 4.1. A measurable function ϕ=ϕ(z,t)L1(ΠT) is an entropy solution of type V of the initial-value problem Eq (4.8), Eq (1.2a) if it satisfies the following conditions:

    (ⅰ) The function ϕ belongs to L(ΠT); for a.e. (z,t)ΠT there holds ϕ(z,t)[0,1].

    (ⅱ) The function ϕ is a weak solution of Eq (4.8), i.e., for all smooth test functions ζ with compact support in ΠT,

    ΠT(ϕtζ+J(ϕ,z)zζ)dzdt=0. (4.13)

    (ⅲ) For all l=0,,K+2, for any nonnegative smooth test function ζ(l) with compact support in Π(l)T and all c[0,1] there holds

    ΠT(|ϕc|tζ(l)+sgn(ϕc)(J(ϕ,z)J(c,z))zζ(l))dzdt+R|ϕ0c|ζ(l)(z,0)dt0. (4.14)

    (ⅳ) The following Kružkov-type [37] entropy inequality holds for all nonnnegative smooth test functions ζ with compact support in ΠT and all constants cR:

    ΠT(|ϕc|tζ+sgn(ϕc)(J(ϕ,z)J(c,z))zζ)dzdt+T0zZ|J(c,z+)J(c,z)|ζ(z,t)dt0. (4.15)

    Note that the entropy inequality Eq (4.15) does not imply the weak formulation Eq (4.13). The standard derivation of the weak formulation from the Kružkov entropy inequality (e.g., [30,Section 2.1]) does not apply here since some of the flux differences |J(c,z+)J(c,z)| are not compactly supported with respect to c, cf. [13,Rem. 1.1].

    Lemma 4.1. Consider ϕ-Scheme 2 applied to Eq (4.8), Eq (1.2a). There exists a constant C1, depending on TV(ϕ0), such that

    ΔziZ|ϕn+1i1/2ϕni1/2|ΔziZ|ϕ1i1/2ϕ0i1/2|C1Δt.

    Proof. Subtracting from Eq (4.1) its version from the previous time step, we get

    ϕn+1i1/2ϕni1/2=(ϕni3/2ϕn1i3/2)λBni1/2+(ϕni1/2ϕn1i1/2){1λBni+1/2+λCni1/2}+(ϕni+1/2ϕn1i+1/2){λCni+1/2},

    where we define

    Bni1/2:=q+i1+γi1W(ϕni1/2),Cni+1/2:={qi+γiϕn1i1/2W(ϕni+1/2)W(ϕn1i+1/2)ϕni+1/2ϕn1i+1/2if ϕni+1/2ϕn1i+1/2,0otherwise.

    Clearly Bni1/20, Cni+1/20, and due to the CFL condition,

    1λBni+1/2+λCni1/20,

    hence taking absolute values and summing over iZ we get, by appealing to standard arguments, that

    ΔziZ|ϕn+1i1/2ϕni1/2|ΔziZ|ϕni1/2ϕn1i1/2|ΔziZ|ϕ1i1/2ϕ0i1/2|.

    Furthermore, following the lines e.g. of the proof of [13,Lemma 3.2], we get that there exists a constant C2 that is independent of (Δt,Δz) such that

    iZ|ϕ1i1/2ϕ0i1/2|C2(TV(ϕ0)+TV(q)+TV(γ)),

    which completes the proof.

    A straightforward calculation yields that we can write the scheme in the form

    ϕn+1i1/2=ϕni1/2+CniΔ+ϕni1/2Dni1Δϕni1/2θni,

    where we define

    Cni:={λqiλγi1ϕni3/2ΔW(ϕni+1/2)Δ+ϕni1/2if Δ+ϕni1/20,λqiotherwise, Dni1:=λq+i+λγi1W(ϕni+1/2),θni:=λ(ϕni1/2Δqi+ϕni3/2Δq+i+ϕni1/2W(ϕni+1/2)ΔγiΔhF,i).

    The incremental coefficients satisfy Cni0 and Dni0; furthermore, the CFL condition ensures that Cni+Dni1 (in all cases for all i and n). Notice that θni=0 with the possible exception for those indices i at which Δqi0, Δq+i0, or Δγi0. According to the definition of γi and that of qi, see Eq (4.3), this may occur at most at a finite number of indices. Precisely, we may assert that (see Eq (4.10))

    θni=0if zi1,ziZ,

    hence for all indices i with the exception of finitely many indices i such that |zjζ|Δz for some ζZ, the scheme is given by the incremental form

    ϕn+1i1/2=ϕni1/2+CniΔ+ϕni1/2Dni1Δϕni1/2

    with incremental coefficients Cni0, Dni0, and Cni+Dni1. This property, in conjunction with Lemma 4.1, shows that we may apply [15,Lemma 5.3] (which is essentially Lemma 4.2 of [6], where a proof can be found) to the situation at hand. From [15,Lemma 5.3] we deduce the following lemma, where Vba(g) denotes the total variation of a function zg(z) over the interval (a,b).

    Lemma 4.2. Consider ϕ-Scheme 2 applied to Eq (4.8), Eq (1.2a). For any interval [a,b] such that [a,b]Z= and any t[0,T] there exists a total variation bound

    Vba(ϕΔz(,t))C(a,b),

    where C(a,b) is independent of (Δx,Δt) and t for t[0,T].

    Finally, we have shown in Theorem 3.1 that ϕ-Scheme 1 Eq (3.4) is monotone. This applies, in particular, to the reduced ϕ-Scheme 2 Eq (4.1) or equivalently, Eq (4.5) or Eq (4.7). Thus, ϕ-Scheme 2 satisfies a discrete entropy inequality. The proof of the following lemma is identical to that of [36,Lemma 5.2], and is therefore omitted.

    Lemma 4.3. The scheme Eq (4.7) (ϕ-Scheme 2) satisfies the following entropy inequality for any ci3/2,ci1/2,ci+1/2[0,1]:

    |ϕn+1i1/2ci1/2||ϕni1/2ci1/2|λΔHniλsgn(ϕn+1i1/2ci1/2)Δh(zi,ϕni+1/2,ϕni1/2),

    where h is defined in Eq (4.6) and the numerical entropy flux Hni is defined by

    Hni:=h(zi,ϕni+1/2ci+1/2,ϕni1/2ci1/2)h(zi,ϕni+1/2ci+1/2,ϕni1/2ci1/2).

    We now may appeal to the results of [36] and argue as follows. Theorem 3.1 and Lemmas 4.1–4.3 ensure convergence of the functions ϕΔz to a weak solution of Eq (4.8), Eq (1.2a) that satisfies items (ⅰ), (ⅱ) and (ⅲ) of Definition 4.1. It also satisfies the entropy inquality Eq (4.15) arising in part (ⅳ) of that definition by utilizing the discrete entropy inequality stated in Lemma 4.3. Thus, we have proved the following theorem.

    Theorem 4.1. Suppose that assumptions Eq (1.3) to Eq (1.5) are in effect and that ϕΔz is defined by Eq (3.5), where the values ϕni1/2 are defined by the scheme Eq (4.5) (that is, ϕ-Scheme 2). Let Δt,Δz0 with λ=Δt/Δz=const. such that the CFL condition Eq (3.2) is satisfied. Then ϕΔz converges in L1loc(ΠT) and a.e. in ΠT to an entropy solution of type V of the initial-value problem Eq (4.8), Eq (1.2a).

    Next, we deal with the marching formula Eq (4.2). To this end, we define a feed flux ˜hF,i exactly as in Eq (4.4) but with ϕF,i replaced by ψF,i for i=1,,K. Furthermore, we recall that ˜vϕs(ϕ)=W(ϕ)/(1ϕ). Thus, the scheme can be written as

    ψn+1i1/2=ψni1/2λΔ(˜hF,i+ψni1/2q+i+ψni+1/2qi+γi(Gni(ψni1/2,ψni+1/2)ϕni1/2˜vϕs(ϕni+1/2)ψni+1/2)). (4.16)

    Lemma 4.4 (Crandall and Tartar [19]). Assume that (Ω,μ) is some measure space and that D is a subset of L1(Ω) with the property that if u,vD, then (uv)=max{u,v}D. Assume that T is a map T:DuT(u)D such that

    ΩT(u)dμ=Ωudμfor all uD.

    Then the following statements, valid for all u,vD, are equivalent:

    (ⅰ) If uv, then T(u)T(v).

    (ⅱ) Ω((T(u)T(v))0)dμΩ((uv)0)dμ.

    (ⅲ) Ω|T(u)T(v)|dμΩ|uv|dμ.

    Following, for instance, [18], we utilize Lemma 4.4 for the following mapping. Assume that DL1(R) is the set of piecewise constant functions and that are constant on the intervals Ii1/2 for iZ, and that with the marching formula Eq (3.15) we associate an operator Kn:DD such that if ψΔz(,tn) is the piecewise constant function defined by Eq (3.13) for t=tn, we may write ψ-Scheme 2 as

    ψΔz(,tn+1)=Kn(ψΔz(,tn)).

    Clearly, the monotonicity of the scheme implies that if u,vD, then

    uvKn(u)Kn(v).

    Thus, Lemma 4.4 (ⅰ) holds. For u=ψΔz(,tn) and v=ψΔz(,tn1), Lemma 4.4 (ⅲ) implies that

    ΔziZ|ψn+1i1/2ψni1/2|=R|ψΔz(z,tn+1)ψΔz(z,tn)|dzR|ψΔz(z,tn)ψΔz(z,tn1)|dz=ΔziZ|ψni1/2ψn1i1/2|

    and therefore

    ΔziZ|ψn+1i1/2ψni1/2|ΔziZ|ψ1i1/2ψ0i1/2|.

    However, we may assert that there exists a constant C3, which is independent of (Δt,Δx), such that

    iZ|ψ1i1/2ψ0i1/2|=iZ|Δ(ψ0i1/2q+i+ψ0i+1/2qi+γi(G0i(ψ0i1/2,ψ0i+1/2)ϕni1/2ψ0i+1/2W(ϕ0i+1/2)1ϕ0i+1/2))λKk=1QF,kAψ0F,kδk,i1/2|C3.

    Since Eq (1.4) is a sufficient condition for this bound on the initial discrete divergence to hold, we get

    ΔziZ|ψn+1i1/2ψni1/2|ΔzC3=ΔtλC3.

    Consequently, we have proved the following lemma.

    Lemma 4.5. There exists a constant C4 that is independent of (Δt,Δz) such that the numerical approximations produced by ψ-Scheme 2 satisfy

    ΔziZ|ψn+1i1/2ψni1/2|C4Δt.

    To write down Scheme 1 in the simplest setting possible, we consider the model and numerical scheme under the assumptions before, and additionally assume a constant bulk velocity q, that the feed terms (giving rise to the singular source) are zero, and set the parameter γ=1. Thus, the model reduces to the triangular system of conservation laws Eq (1.8) with the initial conditions Eq (1.2) (Model 3), where we recall that assumptions Eq (1.3) are in effect, and Scheme 2, which in turn is a reduced form of Scheme 1, now is further reduced to Scheme 3.

    Assume now that η=η(ψ) is a smooth convex entropy function and Q=Q(ϕ,ψ) is the corresponding compatible entropy flux compatible with Eq (1.8b), i.e.,

    ψQ(ϕ,ψ)=η(ψ)ψ˜F(ϕ,ψ). (5.1)

    In what follows, we refer to (η,Q) as an entropy pair for Eq (1.8b). In particular we denote by (η0,Q0) the Kružkov entropy pair [37], that is

    η0(ψ)=|ψc|,Q0(ϕ,ψ)=sgn(ψc)(˜F(ϕ,ψ)˜F(ϕ,c)), (5.2)

    where cR is a constant.

    The convergence proof is based on the following lemma, slightly adapted from [18,Lemma 2.2], which in turn is an adaptation of [43,Theorem 5] (see also [44]).

    Lemma 5.1. Let ϕ be the unique entropy solution of the initial-value problem Eq (1.8a), Eq (1.2a), and assume that {ψν}ν>0 is a family of functions defined on ΠT. If {ψν} is bounded in L(ΠT) and {tη0(ψν)+zQ0(ϕ,ψν)}ν>0 lies in a compact set of H1loc(ΠT) for all constants c, then there exists a sequence {νn}nN such that νn0 as n and a function ψL(ΠT) such that

    ψνnψa.e. and in Lploc(ΠT),1p<.

    Consistently with Eq (4.6), Eq (4.7) we assume that the scheme employed to approximate entropy solutions of Eq (1.8a) is ϕ-Scheme 3, that is,

    ϕn+1i1/2=ϕni1/2λΔh(ϕni+1/2,ϕni1/2),h(v,u):=qv+q+u+uW(v).

    Clearly, under a suitable CFL condition, the ϕ-Scheme 3 converges to the unique entropy solution of Eq (1.8a), Eq (1.2a). Our goal is to establish convergence of the corresponding scheme for ψ (ψ-Scheme 3). We here write the scheme as

    ψn+1i1/2=ψni1/2λΔF(ϕni1/2,ϕni+1/2,ψni1/2,ψni+1/2)ψni1/2λΔF(ϕni,ψni), (5.3)

    where we define the four-argument numerical flux

    F(a,b,u,v):=q+u+qv+(G(a,b,u,v)a˜vϕs(b)v), (5.4)

    denote pairs of neighboring ϕ- and ψ-values by

    ϕni:=(ϕni1/2,ϕni+1/2)andψni:=(ψni1/2,ψni+1/2),

    and replace the arguments "ϕni1/2,ϕni+1/2" by ϕni (analogously for ψ). In Eq (5.4) a and b play the roles of ϕni1/2 and ϕni+1/2, and u and v those of ψni1/2 and ψni+1/2, respectively, and we define G(a,b,u,v) as follows (cf. Eq (3.7), Eq (3.8)). Let

    f(a,b,ψ):=ψ˜V(ψ1(ab)),

    then G(a,b,,) is the Engquist-Osher numerical flux [26] associated with f(a,b,).

    The compensated compactness approach strongly depends on entropy inequalities satisfied by the scheme Eq (5.3). To prepare for the derivation of suitable uniform estimates, we multiply the scheme Eq (5.3) by η(ψn+1i1/2), where η is a smooth convex entropy function, and utilize the Taylor expansion

    η(ψn+1i1/2)(ψn+1i1/2ψni1/2)=η(ψn+1i1/2)η(ψni1/2)+12η(ξn+1/2i1/2)(ψn+1i1/2ψni1/2)2,

    where ξn+1/2i1/2 is an intermediate value between ψni1/2 and ψn+1i1/2. This yields

    η(ψn+1i1/2)η(ψni1/2)+12η(ξn+1/2i1/2)(ψn+1i1/2ψni1/2)2=λη(ψn+1i1/2)ΔF(ϕni,ψni)=λη(ψni1/2)ΔF(ϕni,ψni)λ(η(ψn+1i1/2)η(ψni1/2))ΔF(ϕni,ψni). (5.5)

    We now define the functions ˆf and ˇf as the partial derivatives

    ˆf(a,b,u):=uF(a,b,u,v)=q++uG(a,b,u,v)0,ˇf(a,b,v):=vF(a,b,u,v)=q+(vG(a,b,u,v)a˜vϕs(b)0.

    The dependence of uF(a,b,u,v) and vF(a,b,u,v) on u and v only, respectively, is crucial for the subsequent analysis. We define the functions

    ˆF(a,b,u):=u0ˆf(a,b,s)ds,ˇF(a,b,v):=v0ˇf(a,b,s)ds

    and note that

    F(a,b,u,v)=ˆF(a,b,u)+ˇF(a,b,v). (5.6)

    Next, we define

    ˆQ(a,b,ψ):=ψ0η(u)ˆf(a,b,u)du,ˇQ(a,b,ψ):=ψ0η(v)ˇf(a,b,v)dv,Q(a,b,ψ1,ψ2):=ˆQ(a,b,ψ1)+ˇQ(a,b,ψ2). (5.7)

    The function Q is a consistent numerical entropy flux for the scheme Eq (5.3) for the entropy function η since

    Q(a,a,ψ,ψ)=ψ0η(u)(ˆf(a,a,u)+ˇf(a,a,u))du=ψ0η(u)uF(a,a,u,u)du=ψ0η(u)˜F(a,u)du=Q(a,ψ).

    Furthermore, integration by parts yields

    ˆQ(a,b,ψ)ˆQ(a,b,˜ψ)=η(ψ)(ˆF(a,b,ψ)ˆF(a,b,˜ψ))ψ˜ψη(u)(ˆF(a,b,u)ˆF(a,b,˜ψ))du, (5.8)
    ˇQ(a,b,ψ)ˇQ(a,b,˜ψ)=η(ψ)(ˇF(a,b,ψ)ˇF(a,b,˜ψ))ψ˜ψη(u)(ˇF(a,b,u)ˇF(a,b,˜ψ))du (5.9)
    =η(˜ψ)(ˇF(a,b,ψ)ˇF(a,b,˜ψ))ψ˜ψη(u)(ˇF(a,b,u)ˇF(a,b,ψ))du. (5.10)

    Now, denoting by Δϕ and Δψ difference operators that act on both ϕ- and ψ-arguments only, respectively (leaving the two others unchanged), we observe that

    ΔF(ϕni,ψni)=ΔψF(ϕni,ψni)+Δϕ+F(ϕni1,ψni1). (5.11)

    In light of Eq (5.8) and Eq (5.10),

    η(ψni1/2)ΔψF(ϕni,ψni)=ˆQ(ϕni,ψni1/2)ˆQ(ϕni,ψni3/2)+ˇQ(ϕni,ψni+1/2)ˇQ(ϕni,ψni1/2)(η(ψni1/2)(ˆF(ϕni,ψni1/2)ˆF(ϕni,ψni3/2))ψni1/2ψni3/2η(u)(ˆF(ϕni,u)ˆF(ϕni,ψni3/2))du)(η(ψni1/2)(ˇF(ϕni,ψni+1/2)ˇF(ϕni,ψni1/2))+ψni1/2ψni+1/2η(u)(ˇF(ϕni,u)ˇF(ϕni,ψni+1/2))du)+η(ψni1/2)(ˆF(ϕni,ψni1/2)ˆF(ϕni,ψni3/2)+ˇF(ϕni,ψni+1/2)ˇF(ϕni,ψni1/2))=Q(ϕni,ψni1/2,ψni+1/2)Q(ϕni,ψni3/2,ψni1/2)+Θni1/2=ΔψQ(ϕni,ψni)+Θni1/2, (5.12)

    where the notation for evaluations and differences for Q is the same as for F and Θni1/2:=ˆΘni1+ˇΘni, where

    ˆΘni1:=ψni1/2ψni3/2η(u)(ˆF(ϕni,u)ˆF(ϕni,ψni3/2))du,ˇΘni:=ψni1/2ψni+1/2η(u)(ˇF(ϕni,u)ˇF(ϕni,ψni+1/2))du.

    Since ˆF is increasing and ˇF is decreasing in the respective third argument, there holds ˆΘni1,ˇΘni0 and therefore Θni1/20. Furthermore, we notice that

    η(ψni1/2)Δϕ+F(ϕni1,ψni1)=Δϕ+(η(ψni1/2)F(ϕni1,ψni1)). (5.13)

    From Eq (5.11) we obtain by taking into account Eq (5.12) and Eq (5.13)

    η(ψni1/2)ΔF(ϕni,ψni)=ΔψQ(ϕni,ψni)+η(ψni1/2)Δϕ+F(ϕni1,ψni1)+Θni1/2=ΔQ(ϕni,ψni)Δϕ+Q(ϕni1,ψni1)+η(ψni1/2)Δϕ+F(ϕni1,ψni1)+Θni1/2=ΔQ(ϕni,ψni)+Δϕ+(η(ψni1/2)F(ϕni1,ψni1)Q(ϕni1,ψni1))+Θni1/2.

    Consequently, Eq (5.5) can be written as

    η(ψn+1i1/2)η(ψni1/2)+12η(ξn+1/2i1/2)(ψn+1i1/2ψni1/2)2+λΘni1/2=λΔQ(ϕni,ψni)λ(η(ψn+1i1/2)η(ψni1/2))ΔF(ϕni,ψni)λΔϕ+(η(ψni1/2)F(ϕni1,ψni1)Q(ϕni1,ψni1)). (5.14)

    Multiplying Eq (5.14) by Δz and summing over (n,i)I1, where

    Ik:={(n,i)n=0,,NTk,iZ},

    we get

    ΔziZη(ψNi1/2)ΔziZη(ψ0i1/2)+Δz2I1η(ξn+1/2i1/2)(ψn+1i1/2ψni1/2)2+λΔzI1Θni1/2=λΔzI1ΔQ(ϕni,ψni)λΔzI1(η(ψn+1i1/2)η(ψni1/2))ΔF(ϕni,ψni)λΔzI1Δϕ+(η(ψni1/2)F(ϕni1,ψni1)Q(ϕni1,ψni1)),

    which implies the inequality

    ΔziZη(ψNi1/2)+Δz2I1η(ξn+1/2i1/2)(ψn+1i1/2ψni1/2)2+λΔzI1Θni1/2ΔziZη(ψ0i1/2)+2ηLΔzΔtI11Δz|ΔF(ϕni,ψni)|+CΔzΔtI11Δz|Δϕ+(η(ψni1/2)F(ϕni1,ψni1)Q(ϕni1,ψni1))|.

    The last term on the right-hand side is uniformly bounded since ϕΔz has bounded variation. Now let us choose η(v)=v2 and take into account [34] that there exists a constant CF such that

    ˆΘni11CF(ΔψˆF(ϕni,ψi1/2))2,ˇΘni1CF(Δψ+ˇF(ϕni,ψi1/2))2.

    Noticing that Lemma 4.5, applied to the present scheme (ψ-Scheme 3), implies the bound on the discrete divergence of the numerical flux

    ΔziZ1Δz|ΔFi(ϕn,ψn)|C, (5.15)

    we obtain from Eq (5.14)

    ΔziZη(ψNi1/2)+ΔzI1(ψn+1i1/2ψni1/2)2+λCFΔzI1((ΔψˆF(ϕni,ψi1/2))2+(Δψ+ˇF(ϕni,ψni1/2))2)ΔziZ(ψ0i1/2)2+CT. (5.16)

    Inequality Eq (5.16) implies the following estimate.

    Lemma 5.2. Consider numerical approximations produced by Scheme 3. There exists a constant C7 that is independent of (Δz,Δt) such that

    ΔtΔzNT1n=0iZ(ψn+1i1/2ψni1/2)2C7Δz. (5.17)

    Proof. The estimate for the "time variation" of ψΔz, Eq (5.17), follows immediately from Eq (5.16) if we consider that its right-hand side is uniformly bounded.

    Before proceeding, we prove the following lemma that is crucial for the subsequent analysis. For ease of notation we define the difference operators Δ(3)± and Δ(4)± that only act on the third or fourth argument of a function, respectively.

    Lemma 5.3. Consider numerical approximations produced by Scheme 3. There exist constants C8 and C9 that are independent of (Δz,Δt) such that for all i,

    |ΔQ(ϕni,ψni)|C8|(Δ(3)+Δ(4))F(ϕni,ψni)|+C9(|Δϕni1/2|+|Δ+ϕni1/2|). (5.18)

    Proof. We note that

    ΔQ(ϕni,ψni)=ΔψQ(ϕni,ψni)+Δϕ+Q(ϕni1,ψni1). (5.19)

    We first discuss

    ΔψQ(ϕni,ψni)=Δ(3)ˆQ(ϕni,ψni1/2)+Δ(3)ˇQ(ϕni,ψni+1/2).

    From Eq (5.8) we get

    |Δ(3)ˆQ(ϕni,ψni1/2)|=|η(ψni1/2)Δ(3)F(ϕni,ψni)ψni1/2ψni3/2η(u)(ˆF(ϕni,u)ˆF(ϕni,ψni1/2))du||η(ψni1/2)||Δ(3)F(ϕni,ψni)|+|ψni1/2ψni3/2η(u)du||Δ(3)F(ϕni,ψni)|3ηL(0,1)|Δ(3)F(ϕni,ψni)|

    and analogously

    |Δ(3)ˇQ(ϕni,ψni+1/2)|3ηL(0,1)|Δ(4)F(ϕni,ψni)|,

    hence

    |ΔψQ(ϕni,ψni)|3ηL(0,1)|(Δ(3)+Δ(4))F(ϕni,ψni)|. (5.20)

    On the other hand, we take into account that

    Δϕ+Q(ϕni1,ψni1)=Δϕ+ˆQ(ϕni1,ψni3/2)+Δϕ+ˇQ(ϕni1,ψni1/2).

    Now

    Δϕ+ˆQ(ϕni1,ψni3/2)=ψni3/20η(u)(ˆf(ϕni,u)ˆf(ϕni1,u))du=[η(u)(ˆF(ϕni,u)ˆF(ϕni1,u))]u=ψni3/2u=0ψni3/20η(u)(ˆF(ϕni,u)ˆF(ϕni1,u))du (5.21)

    and analogously

    Δϕ+ˇQ(ϕni1,ψni1/2)=ψni1/20η(u)(ˇf(ϕni,v)ˇf(ϕni1,v))dv=[η(v)(ˇF(ϕni,v)ˇF(ϕni1,v))]v=ψni1/2v=0ψni1/20η(v)(ˇF(ϕni,v)ˇF(ϕni1,v))dv.

    Consequently,

    |Δϕ+ˆQ(ϕni1,ψni3/2)|3ηL(0,1)max0uψni3/2|ˆF(ϕni,u)ˆF(ϕni1,u)|, (5.22)

    and by analogous reasoning for ˇQ,

    |Δϕ+ˇQ(ϕni1,ψni1/2)|3ηL(0,1)max0vψni1/2|ˇF(ϕni,v)ˇF(ϕni1,v)|. (5.23)

    To estimate the right-hand side of Eq (5.22), we recall that

    ˆF(ϕni,u)ˆF(ϕni1,u)=u0(((fni))+(s)((fni1))+(s))ds=:ˆDni1/2.

    We now assume that σ=1 and use Eq (3.9). To estimate ˆDni1/2, we will appeal to the trivial but useful inequality |(αx)(βy)||αβ|+|xy|. We proceed by discussing all possible cases of the location of u in relation to the extrema ˆψni and ˆψni1 of fni and fni1, respectively, and assume that σ=1.

    1. Assume that uˆψniˆψni1. Since ˆψniψnmax,i and ˆψni1ψnmax,i1, in this case ˆDni1/2=0 if ψnmax,i=ψnmax,i1 and otherwise

    |ˆDni1/2|=|fni(u)fni1(u)|=u|˜V(u/ψnmax,i)˜V(u/ψnmax,i1)|u2ψnmax,iψnmax,i1˜V|ψnmax,iψnmax,i1|˜V|ψnmax,iψnmax,i1|.

    Noticing that

    |ψnmax,iψnmax,i1|=|ϕni1/2ϕni+1/2ϕni3/2ϕni1/2||Δϕni1/2|+|Δ+ϕni1/2|, (5.24)

    we conclude that

    |ˆDni1/2|˜V(|Δϕni1/2|+|Δ+ϕni1/2|).

    2. If ˆψni<u<ˆψni1 then

    |ˆDni1/2|=|fni(ˆψni)fni1(u)||fni(ˆψni)fni1(ˆψni1)|+|fni1(ˆψni1)fni1(u)|.

    Since fni(ˆψni)=ψnmax,iω˜V(ω) for all i, we conclude that

    |fni(ˆψni)fni1(ˆψni1)|ω˜V|ψnmax,iψnmax,i1|˜V(|Δϕni1/2|+|Δ+ϕni1/2|). (5.25)

    On the other hand, in the present case

    |fni1(ˆψni1)fni1(u)|=fni1(ˆψni1)fni1(u)fni1(ˆψni1)fni1(ˆψni).

    Since for s[0,ˆψni1] there holds (fni1)(s)(fni1)(0)=˜V(0), we get

    |fni1(ˆψni1)fni1(u)|=ˆψni1u(fni1)(s)dsˆψni1ˆψni(fni1)(s)ds˜V(0)(ˆψni1ˆψni).

    Lemma 3.1 (a) implies that

    |ˆψniˆψni1|=ω|ψnmax,iψnmax,i1|, (5.26)

    hence

    |ˆDni1/2|2˜V(|Δϕni1/2|+|Δ+ϕni1/2|).

    The same estimate holds for ˆψni1<u<ˆψni.

    3. If uˆψniˆψni1 then utilizing Eq (5.25) we get

    |ˆDni1/2|=|fni(ˆψni)fni1(ˆψni1)|˜V(|Δϕni1/2|+|Δ+ϕni1/2|).

    Combining all possible cases we deduce that

    |ˆDni1/2|(˜V+2˜V)(|Δϕni1/2|+|Δ+ϕni1/2|). (5.27)

    Next, we deal with Eq (5.23), recalling that (see Eq (5.4))

    ˇF(ϕni,v)ˇF(ϕni1,v)=ˇDni1/2+Eni1/2,

    where we define

    ˇDni1/2:=v0(((fni))(s)((fni1))(s))ds,Eni1/2:=(ϕni1/2˜vϕs(ϕni+1/2)ϕni3/2˜vϕs(ϕni1/2))v.

    The discussion of all possible cases of the location of v in relation to ˆψni and ˆψni1 gives rise to the following cases for the estimation of ˇDni1/2.

    1. If vˆψniˆψni1 or vψnmax,iψnmax,i1 then ˇDni1/2=0.

    2. To handle the case ˆψniˆψni1vˆψniˆψni1 we assume that ˆψni<ˆψni1 and ˆψnivˆψni1. Then

    |ˇDni1/2|=|fni(v)fni(ˆψni)|=|v^ψni((fni))(s)ds|maxˆψnisˆψni1|((fni))(s)||ˆψniˆψni1||(fni)(ψninfl,i)||ˆψniˆψni1|.

    By Lemma 3.1 (ⅱ),

    so does not depend on and we conclude that

    Applying the argument of Eq (5.24) and Eq (5.26) yields

    The same inequality is also valid if and .

    3. Finally, assume that . In this case

    Taking into account that , we get

    (5.28)

    If , then , hence Eq (5.28) means that

    Suppose now that

    (5.29)

    Since we know that , the inequality can only be satisfied if

    On the other hand,

    so Eq (5.29) is only possible when , which means .

    If instead of Eq (5.29) we have

    (5.30)

    then implies

    The remainder of the estimate is based on Eq (5.24). Since , we conclude that if Eq (5.29) holds, then

    In combination with Eq (5.28) we obtain in this case

    (5.31)

    Next, suppose that instead of Eq (5.29) or Eq (5.30) there holds

    then the discussion of Eq (5.29) can be applied again and we get that this ordering is only feasible if all terms are equal and zero, and therefore . On the other hand, let us assume that

    In this case,

    It remains to treat the case . We then have , and analogously to the derivation of Eq (5.31) we get

    Collecting all estimates for , we see that

    (5.32)

    Furthermore, we obtain

    (5.33)

    Combining the estimates Eq (5.27), Eq (5.32) and Eq (5.33), we obtain from Eq (5.22) and Eq (5.23) the bounds

    (5.34)

    and therefore

    with constants and . Combining the last inequality with Eq (5.19) and Eq (5.20) we arrive at the desired estimate Eq (5.18).

    From Eq (5.18), and considering that for all and , we obtain

    Summing over we get

    Multiplying this inequality by and taking into account Eq (5.16) and the uniform bound on we have proved the following lemma.

    Lemma 5.4. Consider numerical approximations produced by Scheme 3. There exists a constant that is independent of and such that the following estimate holds:

    In part following the proofs of Lemmas 5.5 and 5.9 in [33] and in particular that of Lemma 3.4 in [18] we now prove the compactness result.

    Lemma 5.5. Assume that is generated by the scheme Eq (5.3) (-Scheme 3), and that is the unique entropy solution of Eq (1.8a), Eq (1.2a) on . Furthermore, we denote by the Kružkov entropy pair Eq (5.2), and the distribution

    Then the sequence belongs to a compact subset of .

    The proof of Lemma 5.5 essentially follows the proof of [18,Lemma 3.4], with some slight modifications due the particular definition of the numerical flux in the present case. The proof is presented in Appendix A.

    Since strongly in , we obtain that there exists a constant such that

    hence the sequence , where we define

    is compact in . Now, by Lemma 5.1 there exists a subsequence (which we do not relabel) and a function such that

    (5.35)

    Theorem 5.1. Assume that the maps and are the limit functions of and of as (the latter one being defined by Eq (5.35), that is, we consider Scheme 3). Then is a weak solution of the initial-value problem Eq (1.8), Eq (1.2) in the sense of Definition 1.1.

    The proof follows that of [18,Lemma 3.5], again with slight modifications. We refer to Appendix A for details.

    To simplify error estimations we utilize a grid with the property that the boundaries of the tank agree with the boundaries of a cell (see Figure 2 (left)). Since an exact solution is frequently difficult to obtain, we use an approximate reference solution obtained with a large number cells against which the error of other simulated solutions with is measured. The error is estimated on a fixed interval slightly larger than the column of height so that the outflow volume fractions are included. We define the coarsest grid of cells with and place the column between and . This corresponds to Figure 2 (left) with . We define the length of the interval of error estimation as .

    To estimate the convergence order, we simulate with cells, , where the integer defines the number of cells of the reference solution. Then we define , and the factor of refinement from cells to as . We note that for all .

    We will now measure the error between the piecewise constant numerical solution obtained by cells (we skip the index for a moment) and the reference solution obtained with cells on the grid refined by a factor . The refined grid satisfies and we have . The corresponding numerical solutions on the refined grids are denoted by (skipping the time index ) , , etc., where are defined by . The refined cells are numbered such that the first cell for above contains the value . Then . This means that the cells within contain the values , and analogously for ; see Figure 2 (right).

    Note that the location of the spatial discontinuities and will coincide with a cell boundary for any mesh considered in the refinement process while the locations of the inlets , etc. will be chosen in such a way that each of them lies inside a cell for the finest grid; hence, they do this for all the coarser meshes. In this way, the numerical fluxes at cell boundaries are well defined.

    We compute the estimated error at a time point and define

    The -difference between two numerical solutions computed on grids with cell sizes and is calculated as follows for :

    with the summands defined by

    The approximate relative error for in the interval is then defined as

    We define analogously and hence, the total relative error can be defined as

    and the observed convergence order between two discretizations and as

    For smooth solutions and a constant (see Eq (1.7)), we also use an alternative way [7] of calculating approximate errors and convergence orders in which a reference solution is not needed. One can use cubic interpolation to compute the quantities from a grid with cells, , with an integer, taking into consideration that . Then, is given by

    The alternative approximate relative -error for can then be calculated as

    We can define and analogously along with the alternative total approximate -error and convergence order

    For the first example, in Section 6.3, we use a smooth solution away from spatial discontinuities, to estimate the order of convergence of the numerical scheme. For this example, we use , for and ; hence, and .

    In Sections 6.4 and 6.5, we exemplify counter- and co-current flows of the primary and secondary disperse phases, respectively. For these two examples, we use , and . We set three inlets , and dividing the tank into four equal parts each with the height , where is used. These three inlets are defined so that they lie inside a cell for any mesh size considered. A fixed quantity of the is introduced through inlet , a fixed quantity of the secondary disperse phase through inlet and some wash water through inlet .

    Tables 1 and 2 show the estimated errors and convergence orders for the three scenarios studied. In the calculations of the alternative approximate error and convergence order in Table 1, we use .

    Table 1.  Simulation of a smooth solution (Section 6.3). Total estimated relative -error , alternative relative -error , estimated convergence order , and its alternative counterpart , calculated with and .
    500 ——
    1000
    2000
    4000
    8000
    16000 ——

     | Show Table
    DownLoad: CSV
    Table 2.  Approximate total relative -errors and convergence orders calculated between consecutive values of , with (a) for Application 1 (counter-current flow) at simulated time , (b) for Application 2 (co-current flow) at simulated time .

     | Show Table
    DownLoad: CSV

    We consider a vessel with a constant cross-sectional area of , and we set all inlet and outlet volumetric flows to zero, i.e, . (Under these assumptions, the scheme reduces to Scheme 3 for Model 3.) For the velocity functions and , given by Eq (2.10) and Eq (2.11), respectively, we use the parameters , , and , and consider (counter-current flow). The initial datum is a sinusoidal function for both phases with support in the interval .

    We simulate a short time, until , before the first discontinuity appears; see the first row in Figure 4 where is used. In the second row, we compare two approximate solutions obtained with a coarse mesh with and a finer one with . Table 1 shows the estimated errors and convergence orders. Both and assume values close to one as increases, as expected, confirming that the scheme is first-order accurate for smooth solution.

    Figure 4.  Simulation of a smooth solution (Section 6.3). First row: Time evolution of the volume fractions of the primary disperse phase (left) and the secondary disperse phase (right) from to . Second row: Approximate solutions at time computed with (left) and (right).

    The remainder of examples refer to Model 1 discretized by Scheme 1. We first illustrate that the "crossing condition" is satisfied as mentioned in Section 4.1. For this we use the constant and simulate a tank that initially contains only water, i.e., for all . At we start pumping aggregates, solids, fluid and wash water with , , , , and . We choose the volumetric flows , so that the volumetric flows in the tank are positive in all zones but not in zone 1. Three inlets , and divide the tank into four zones of equal height.

    Figure 5 shows the graphs of the flux functions on both sides of each discontinuity in , for this case three inlets and two outlets . We see that the fluxes (defined in Eq (4.9)) intersect when , and , and do not intersect in at neither nor . Figure 6 shows the simulation results during .

    Figure 5.  Illustration of the crossing condition (Section 6.4). The crossing condition is satisfied at each of the five spatial discontinuities.
    Figure 6.  Simulation of the example in Section 6.4 during .

    We consider now a complete tank with ; hence, the primary disperse phase will move upwards and the secondary disperse phase downwards with respect to the volume average velocity of the mixture. A straightforward interpretation of this scenario is the flotation process used in the mineral industry to recover valuable minerals from crushed ore; see the model in [8,9]. In that model, the primary disperse phase consists of aggregates, which are air bubbles fully loaded with hydrophobic minerals, and the secondary disperse phase is the tailings, consisting of hydrophilic particles suspended in the fluid that do not attach to air bubbles. We consider three inlets , and , dividing the tank into four regions with equal height. At , only gas is fed, at only wash water, while at a slurry of solids and water is fed into the column.

    The cross-sectional area is discontinuous (cf. Figure 2) due to a centered pipe from the top down to that introduces material into the tank. It is given by

    These values correspond to the reflux flotation cell studied in [21].

    We consider that the column is initially filled only with fluid, hence for all , when we start pumping aggregates and solids with concentrations , , , , and , along with fluid and/or wash water. We choose , so that the mixture flows in zones 2 and 3 are positive, i.e., directed upwards: in zone 2 and in zone 3.

    Figure 7 shows the time evolution of the volume fractions of and . It can be seen that the aggregates rise fast to the top, while the solids are travelling both up and down the vessel, leaving through the effluent and the underflow.

    Figure 7.  Application 1: Counter-current flows. Time evolution of the volume fraction profiles of the primary disperse phase (left) and secondary disperse phase (right) from time to seen from two different angles (first and second rows).

    At time , we change the volumetric flow from to . After this change, the solids settle and we obtain a steady state. We mention that this is not a desired steady state in the mining industry (the capacity of the device is not fully used); see [9] for more examples. Table 2 (a) shows the estimated errors and convergence orders for this simulation. As in the smooth example in Section 6.3, the convergence orders tend to one as increases.

    For the last example, we consider , i.e., both the primary and secondary disperse phases have a density smaller than that of the fluid and therefore move upwards relative to the mixture. This scenario could be a flotation process with two buoyant phases differing in density and possibly also in size. We consider here the same flotation column as in Application 1 and choose , , , and so that we have two buoyant phases with different (upwards-directed) velocities relative to the mixture. As in the previous example, only the primary disperse phase is fed into the tank at and only the secondary at . The column is initially filled with only fluid at time , hence for all , when we start pumping both phases with the following volume fractions: , , , , and . We choose the volumetric flows , so that the volumetric flows in the tank are positive in all zones with the exception of zone 1.

    Figure 8 shows the time evolution of the volume fractions of both phases. It can be seen that, for times , the primary disperse phase leaves the tank through both the underflow and effluent outlets, while the secondary disperse phase quickly rises to the top part of the tank and leaves it just through the effluent outlet.

    Figure 8.  Application 2: Co-current flow. Time evolution of the volume fraction profiles of primary disperse phase (left) and secondary disperse phase (right) from time to seen from two different angles (first and second rows).

    At , we change the volumetric flow of the inlet from to , maintaining the other volumetric flows. As a consequence we can see that the primary disperse phase rises and leaves zone 1, exiting the tank only through the effluent while the secondary disperse phase maintains the same behaviour as before and is present only above the feed level . Table 2 (b) shows the estimated errors and convergence orders for this simulation, which have the same behaviour as the ones in the numerical examples in Sections 6.3 and 6.4.

    The present study outlines a numerical method for a triangular system of two PDEs, whose flux functions have several spatial discontinuities due to in- and outflows of a one-dimensional tank with possibly varying cross-sectional area. The triangular structure is utilized in the following way in the numerical scheme. The numerical update formula corresponding to the first scalar equation contains, for the nonlinear term, a numerical flux where the the volume fraction in the left cell is multiplied with the velocity computed in the right cell; see [6]. The update formula for the second equation uses the Engquist-Osher numerical flux for the term modelling the nonlinear relative flux of the secondary disperse phase, chosen in a particular way since this flux also depends on the primary disperse phase volume fraction. The other terms of the second update formula are also chosen in such a way that the entire scheme is proved to be monotone under the CFL condition Eq (3.2). We prove that the numerically obtained volume fractions satisfy the invariant-region property that they stay between zero and one, as is physically expected.

    The numerical scheme is applied to simulate the hydrodynamic movement of simultaneously rising aggregates (air bubbles with attached hydrophobic particles) and settling hydrophilic particles in the fluid under in- and outflows of a flotation column. As a demonstration of the capabilities of the numerical method, three different settings are simulated. The convergence order of the numerical method is estimated. As expected, in regions where the solution is smooth, the order is one. The first-order scheme proposed in this paper could be improved to achieve second-order accuracy, for instance, by techniques of variable extrapolation [6,14].

    In [9], the authors proposed a staggered scheme to compute numerical solutions for a flotation column, following the approach of Karlsen et al. [32,33]. Although the staggered scheme worked for a single inlet for a mixture of aggregates and solids, we have, in the case of several feed inlets, found it difficult to get consistent numerical solutions with respect to different mesh sizes.

    We are currently [12] extending the model and numerical scheme to the explicit description of drainage of liquid from the foam forming at the top of a flotation column. This phenomenon gives rise to a model similar to Eq (1.1) but with an additional degenerating diffusion term. The numerical solution of the resulting system of non-linear convection-diffusion equations calls for semi-implicit discretizations to alleviate the severe restrictions in the CFL condition due to the diffusion term.

    R.B. is supported by ANID (Chile) through Anillo project ANID/PIA/ACT210030; Centro de Modelamiento Matemático (BASAL projects ACE210010 and FB210005); CRHIAM, project ANID/Fondap/15130015; and Fondecyt project 1210610. S.D. is supported by the Swedish Research Council (Vetenskapsrådet, 2019-04601). M.C.M. is supported by grant PID2020-117211GB-I00 funded by MCIN/AEI/10.13039/501100011033 and project CIAICO/2021/227. Y.V. is supported by IFARHU-SENACYT (Panama) and project UCO1866.

    The authors declare there is no conflict of interest.

    Proof of Lemma 5.5. Following [18], we work with smooth entropies instead of , so we denote by a smooth convex approximation to , so that and , and . Moreover, if is the entropy flux associated with , then there also holds as . Then we split as , where we define

    If denotes a test function with compact support, then as in [18], one verifies that

    hence is compact in . By an integration by parts we get

    so we may finally write

    (A.1)

    We define the cell average

    (A.2)

    Replacing the integral in the first term of the right-hand side of Eq (A.1) by produces the following error, where we follow the derivation of Eq (3.27) in [18]:

    (see Eq (5.17)). By similar arguments we obtain the bound

    Consequently, and further following [18], we have shown that

    We now utilize the "scheme for ", Eq (5.14), to rewrite the term in curled brackets as , where we define

    (A.3)

    Thus, where

    and , , and are defined analogously. In view of Lemma 5.2, we get

    and therefore . Appealing to the divergence bound of the numerical flux Eq (5.15) and taking into account the bound on it also follows that , and therefore .

    Finally, to deal with we consider first and let , and denote the entropy and numerical entropy fluxes calculated from Eq (5.1) and Eq (5.7), respectively, where . Since is consistent with ,

    (cf. Eq (5.9)). By using the monotonicity of with respect to its -argument we get

    so in the limit ,

    (A.4)

    Noticing that

    we get from Eq (A.3)

    By a summation by parts and applying Eq (A.4) we get

    (A.5)

    We now write

    where

    From Eq (5.34), and considering in that bound, we deduce there exists a constant such that , therefore there exists (another) constant such that

    (A.6)

    Consequently, from Eq (A.5) we deduce that

    From Eq (5.16) we infer that there exists a constant such that

    Noticing that also

    we conclude that there exists a constant such that

    Next, we deal with . Applying again a summation by parts, we get

    The definition of (see Eq (5.7)) yields

    By a computation similar to Eq (5.21) we get

    where

    The discussion of is similar to that of above, and appealing to Eq (5.34) we see that

    On the other hand, Eq (A.6) implies that

    Thus

    and we deduce that can be bounded in a similar way as . In particular,

    and we conclude that also

    so is compact in . Thus, the sequence , and therefore also the sequence belong to a compact subset of .

    Proof of Theorem 5.1. We only need to verify that is a weak solution of Eq (1.8b), that is, that Eq (1.10) holds. To this end, we choose a test function , recall the definition Eq (A.2) of cell averages , multiply the -scheme Eq (5.3) by , sum over and , and apply a summation by parts to obtain an identity , where

    By exactly following the estimates of terms and in the proof of [18,Lemma 3.5] and appealing to the bounded convergence theorem we may prove that

    (A.7)

    The treatment of differs from that of the term in [18,Lemma 3.5] since here the numerical flux depends on four arguments (not three, as in [18]). We here get

    where we define

    The term can be estimated by choosing a constant such that for and noting that

    (A.8)

    Furthermore, in light of Eq (5.6) the difference arising in can be written as

    Utilizing the estimate Eq (5.34) with yields that there exists a constant such that

    hence

    (A.9)

    To estimate , we utilize Eq (5.16). Then

    (A.10)

    From Eq (A.8), Eq (A.9) and Eq (A.10) and appealing to the strong convergence of and we deduce that

    (A.11)

    The limits Eq (A.7) and Eq (A.11) imply that the limit is a weak solution.



    [1] Adimurthi, G. D. Veerappa Gowda, S. Mishra, Optimal entropy solutions for conservation laws with discontinuous flux, J. Hyperbolic Differ. Equ., 2 (2005), 1–56. https://doi.org/10.1142/S0219891605000361 doi: 10.1142/S0219891605000361
    [2] B. Andreianov, K. H. Karlsen, N. H. Risebro, A theory of -dissipative solvers for scalar conservation laws with discontinuous flux, Arch. Ration. Mech. Anal., 201, (2011), 1–60. https://doi.org/10.1007/s00205-010-0389-4 doi: 10.1007/s00205-010-0389-4
    [3] B. Andreianov, C. Donadello, S. S. Ghoshal, U. Razafison, On the attainable set for a class of triangular systems of conservation laws, J. Evol. Equ., 15 (2015), 503–532. https://doi.org/10.1007/s00028-014-0267-x doi: 10.1007/s00028-014-0267-x
    [4] B. Andreianov, D. Mitrović, Entropy conditions for scalar conservation laws with discontinuous flux revisited, Ann. Inst. H. Poincaré C Anal. Non Linéaire, 32 (2015), 1307–1335. https://doi.org/10.1016/j.anihpc.2014.08.002 doi: 10.1016/j.anihpc.2014.08.002
    [5] A. Bressan, G. Guerra, W. Shen, Vanishing viscosity solutions for conservation laws with regulated flux, J. Differ. Equ., 266 (2018), 312–351. https://doi.org/10.1016/j.jde.2018.07.044 doi: 10.1016/j.jde.2018.07.044
    [6] R. Bürger, A. García, K. H. Karlsen, J.D. Towers, A family of numerical schemes for kinematic flows with discontinuous flux, J. Eng. Math., 60 (2008), 387–425. https://doi.org/10.1007/s10665-007-9148-4 doi: 10.1007/s10665-007-9148-4
    [7] S. Boscarino, R. Bürger, P. Mulet, G. Russo, L. M. Villada, Linearly implicit IMEX Runge-Kutta methods for a class of degenerate convection-diffusion problems, SIAM J. Sci. Comput., 37 (2015), B305–B331. https://doi.org/10.1137/140967544 doi: 10.1137/140967544
    [8] R. Bürger, S. Diehl, M. C. Martí, A conservation law with multiply discontinuous flux modelling a flotation column, Netw. Heterog. Media, 13 (2018), 339–371. https://doi.org/10.3934/nhm.2018015 doi: 10.3934/nhm.2018015
    [9] R. Bürger, S. Diehl, M. C. Martí, A system of conservation laws with discontinuous flux modelling flotation with sedimentation, IMA J. Appl. Math., 84 (2019), 930–973. https://doi.org/10.1093/imamat/hxz021 doi: 10.1093/imamat/hxz021
    [10] R. Bürger, S. Diehl, M. C. Martí, Y. Vásquez, Flotation with sedimentation: Steady states and numerical simulation of transient operation, Minerals Eng., 157 (2020), 106419. https://doi.org/10.1016/j.mineng.2020.106419 doi: 10.1016/j.mineng.2020.106419
    [11] R. Bürger, S. Diehl, M.C. Martí, Y. Vásquez, Simulation and control of dissolved air flotation and column froth flotation with simultaneous sedimentation, Water Sci. Tech., 81 (2020), 1723–1732. https://doi.org/10.2166/wst.2020.258 doi: 10.2166/wst.2020.258
    [12] R. Bürger, S. Diehl, M. C. Martí, Y. Vásquez, A degenerating convection-diffusion system modelling froth flotation with drainage, arXiv: 2202.04679, [Preprint], (2022), [cited 2022 Nov 22 ]. Available from: https://doi.org/10.48550/arXiv.2202.04679
    [13] R. Bürger, K. H. Karlsen, N. H. Risebro, J. D. Towers, Well-posedness in and convergence of a difference scheme for continuous sedimentation in ideal clarifier-thickener units, Numer. Math., 97 (2004), 25–65. https://doi.org/10.1007/s00211-003-0503-8 doi: 10.1007/s00211-003-0503-8
    [14] R. Bürger, K. H. Karlsen, H. Torres, J. D. Towers, Second-order schemes for conservation laws with discontinuous flux modelling clarifier-thickener units, Numer. Math., 116 (2010), 579–617. https://doi.org/10.1007/s00211-010-0325-4 doi: 10.1007/s00211-010-0325-4
    [15] R. Bürger, K. H. Karlsen, J. D. Towers, An Engquist-Osher-type scheme for conservation laws with discontinuous flux adapted to flux connections, SIAM J. Numer. Anal., 47 (2009), 1684–1712. https://doi.org/10.1007/s00211-010-0325-4 doi: 10.1007/s00211-010-0325-4
    [16] R. Bürger, K. H. Karlsen, J. D. Towers, On some difference schemes and entropy conditions for a class of multi-species kinematic flow models with discontinuous flux, Netw. Heterog. Media, 5 (2010), 461–485. https://doi.org/10.3934/nhm.2010.5.461 doi: 10.3934/nhm.2010.5.461
    [17] G. M. Coclite, K. H. Karlsen, S. Mishra, N. H. Risebro, Convergence of vanishing viscosity approximations of triangular systems of multi-dimensional conservation laws, Boll. Unione Mat. Ital., 2 (2009), 275–284.
    [18] G. M. Coclite, S. Mishra, N. H. Risebro, Convergence of an Engquist-Osher scheme for a multi-dimensional triangular system of conservation laws, Math. Comp., 79 (2010), 71–94. https://doi.org/10.1090/S0025-5718-09-02251-0 doi: 10.1090/S0025-5718-09-02251-0
    [19] M. G. Crandall, L. Tartar, Some relations between nonexpansive and order preserving mappings, Proc. Amer. Math. Soc., 78 (1980), 385–390. https://doi.org/10.1090/S0002-9939-1980-0553381-X doi: 10.1090/S0002-9939-1980-0553381-X
    [20] V. G. Danilov, D. Mitrovic, Delta shock wave formation in the case of triangular hyperbolic system of conservation laws. J. Differ. Equ., 245 (2008), 3704–3734. https://doi.org/10.1090/S0025-5718-09-02251-0 doi: 10.1090/S0025-5718-09-02251-0
    [21] J. E. Dickinson, K. P. Galvin, Fluidized bed desliming in fine particle flotation, Part Ⅰ, Chem. Eng. Sci., 108 (2014), 283–298. https://doi.org/10.1016/j.ces.2013.11.006 doi: 10.1016/j.ces.2013.11.006
    [22] S. Diehl, On scalar conservation laws with point source and discontinuous flux function, SIAM J. Math. Anal., 26 (1995), 1425–1451. https://doi.org/10.1137/S0036141093242533 doi: 10.1137/S0036141093242533
    [23] S. Diehl, A conservation law with point source and discontinuous flux function modelling continuous sedimentation, SIAM J. Appl. Math., 56 (1996), 388–419. https://doi.org/10.1137/S0036139994242425 doi: 10.1137/S0036139994242425
    [24] S. Diehl, Operating charts for continuous sedimentation Ⅰ: control of steady states, J. Eng. Math., 41 (2001), 117–144. https://doi.org/10.1023/A:1011959425670 doi: 10.1023/A:1011959425670
    [25] S. Diehl, The solids-flux theory—confirmation and extension by using partial differential equations, Water Res. 42 (2008), 4976–4988. https://doi.org/10.1016/j.watres.2008.09.005 doi: 10.1016/j.watres.2008.09.005
    [26] B. Engquist, S. Osher, One-sided difference approximations for nonlinear conservation laws, Math. Comp., 36 (1981), 321–351. https://doi.org/10.1090/S0025-5718-1981-0606500-X doi: 10.1090/S0025-5718-1981-0606500-X
    [27] J. A. Finch, G. S. Dobby, Column Flotation, London: Pergamon Press, 1990.
    [28] K. P. Galvin, J. E. Dickinson, Fluidized bed desliming in fine particle flotation Part Ⅱ: Flotation of a model feed, Chem. Eng. Sci., 108 (2014), 299–309. https://doi.org/10.1016/j.ces.2013.11.027 doi: 10.1016/j.ces.2013.11.027
    [29] T. Gimse, N. H. Risebro, Riemann problems with a discontinuous flux function. Third International Conference on Hyperbolic Problems, Theory, Numerical Methods and Applications, 1 (1991), 488–502.
    [30] H. Holden, N. H. Risebro, Front Tracking for Hyperbolic Conservation Laws, 2nd edition, Berlin: Springer Verlag, 2015.
    [31] E. L. Isaacson, J. B. Temple, Analysis of a singular hyperbolic system of conservation laws, J. Differ. Equ., 65 (1986), 250–268.
    [32] K. H. Karlsen, S. Mishra, N. H. Risebro, Semi-Godunov schemes for general triangular systems of conservation laws, J. Eng. Math., 60 (2008), 337–349. https://doi.org/10.1007/s10665-007-9163-5 doi: 10.1007/s10665-007-9163-5
    [33] K. H. Karlsen, S. Mishra, N. H. Risebro, Convergence of finite volume schemes for triangular systems of conservation laws, Numer. Math., 111 (2009), 559–589. https://doi.org/10.1007/s00211-008-0199-x doi: 10.1007/s00211-008-0199-x
    [34] K. H. Karlsen, N. H. Risebro, Convergence of finite difference schemes for viscous and inviscid conservation laws with rough coefficients, ESAIM: Math. Model. Numer. Anal., 35 (2001), 239–269. https://doi.org/10.1051/m2an:2001114 doi: 10.1051/m2an:2001114
    [35] K. H. Karlsen, N. H. Risebro, J. D. Towers, stability for entropy solutions of nonlinear degenerate parabolic convection-diffusion equations with discontinuous coefficients, Skr. K. Nor. Vidensk. Selsk., 3 (2003), 1–49.
    [36] K. H. Karlsen, J. D. Towers, Convergence of a Godunov scheme for conservation laws with a discontinuous flux lacking the crossing condition, J. Hyperbolic Differ. Equ., 14 (2017), 671–701. https://doi.org/10.1142/S0219891617500229 doi: 10.1142/S0219891617500229
    [37] S. N. Kružkov, First order quasi-linear equations in several independent variables, Math. USSR Sb., 10 (1970), 217–243.
    [38] G. J. Kynch, A theory of sedimentation, Trans. Faraday Soc., 48 (1952), 166–176. https://doi.org/10.1039/tf9524800166 doi: 10.1039/tf9524800166
    [39] D. Mitrovic, New entropy conditions for scalar conservation laws with discontinuous flux, Discrete Contin. Dyn. Syst., 30 (2011), 1191–1210. https://doi.org/10.3934/dcds.2011.30.1191 doi: 10.3934/dcds.2011.30.1191
    [40] D. Mitrovic, V. Bojkovic, V. G. Danilov, Linearization of the Riemann problem for a triangular system of conservation laws and delta shock wave formation process, Math. Meth. Appl. Sci., 33 (2010), 904–921. https://doi.org/10.1002/mma.1226 doi: 10.1002/mma.1226
    [41] R. Pal, J. H. Masliyah, Flow characterization of a flotation column, Canad. J. Chem. Eng., 67 (1989), 916–923. https://doi.org/10.1002/cjce.5450670608 doi: 10.1002/cjce.5450670608
    [42] R. Pal, J. H. Masliyah, Oil recovery from oil in water emulsions using a flotation column, Canad. J. Chem. Eng., 68 (1990), 959–967. https://doi.org/10.1002/cjce.5450680611 doi: 10.1002/cjce.5450680611
    [43] E. Yu. Panov, Existence and strong pre-compactness properties for entropy solutions of a first-order quasilinear equation with discontinuous flux, Arch. Rational Mech. Anal., 195 (2010), 643–673. https://doi.org/10.1007/s00205-009-0217-x doi: 10.1007/s00205-009-0217-x
    [44] E. Yu. Panov, Erratum to: Existence and strong pre-compactness properties for entropy solutions of a first-order quasilinear equation with discontinuous flux, Arch. Rational Mech. Anal., 196 (2010), 1077–1078. https://doi.org/10.1007/s00205-009-0217-x doi: 10.1007/s00205-009-0217-x
    [45] P. Quintanilla, S. J. Neethling, P. R. Brito-Parada, Modelling for froth flotation control: A review, Min. Eng., 162 (2021), 106718. https://doi.org/10.1016/j.mineng.2020.106718 doi: 10.1016/j.mineng.2020.106718
    [46] J. F. Richardson, W. N. Zaki, Sedimentation and fluidisation: Part Ⅰ, Trans. Instn. Chem. Engrs. (London), 32 (1954), 35–53.
    [47] M. D. Rosini, Systems of conservation laws with discontinuous fluxes and applications to traffic, Ann. Univ. Mariae Curie-Skłodowska Sect. A, 73 (2019), 135–173.
    [48] C. O. R. Sarrico, A distributional product approach to -shock wave solutions for a generalized pressureless gas dynamic system, Int. J. Math., 25, (2014), 145007.
    [49] W. Shen, On the Cauchy problems for polymer flooding with gravitation, J. Differ. Equ., 261 (2016), 627–653. https://doi.org/10.1016/j.jde.2016.03.020 doi: 10.1016/j.jde.2016.03.020
    [50] P. Stevenson, P. S. Fennell, K. P. Galvin, On the drift-flux analysis of flotation and foam fractionation processes, Canad. J. Chem. Eng., 86 (2008), 635–642. https://doi.org/10.1002/cjce.20076 doi: 10.1002/cjce.20076
    [51] J. Vandenberghe, J. Chung, Z. Xu, J. Masliyah, Drift flux modelling for a two-phase system in a flotation column, Canad. J. Chem. Eng., 83 (2005), 169–176.
    [52] Y. Vásquez, Conservation Laws with Discontinuous Flux Modeling Flotation Columns, (Spanish), Doctoral Thesis of Universidad de Concepción, Concepción, 2022.
  • This article has been cited by:

    1. Raimund Bürger, Stefan Diehl, M Carmen Martí, Yolanda Vásquez, A degenerating convection–diffusion system modelling froth flotation with drainage, 2022, 87, 0272-4960, 1151, 10.1093/imamat/hxac033
    2. Li Feng, Yunjuan Jin, Yinzheng Sun, A class of piecewise constant Radon measure solutions to Riemann problems of compressible Euler equations with discontinuous fluxes: pressureless flow versus Chaplygin gas, 2024, 75, 0044-2275, 10.1007/s00033-024-02353-1
    3. Stefan Diehl, Jaime Manríquez, Catherine J. Paul, Tage Rosenqvist, A convection-diffusion-reaction system with discontinuous flux modelling biofilm growth in slow sand filters, 2025, 137, 0307904X, 115675, 10.1016/j.apm.2024.115675
    4. Raimund Bürger, Julio Careaga, Stefan Diehl, Romel Pineda, Numerical schemes for a moving-boundary convection-diffusion-reaction model of sequencing batch reactors, 2023, 57, 2822-7840, 2931, 10.1051/m2an/2023068
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2435) PDF downloads(119) Cited by(4)

Figures and Tables

Figures(8)  /  Tables(2)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog