Loading [MathJax]/jax/output/SVG/jax.js
Research article

A monotone finite volume scheme for linear drift-diffusion and pure drift equations on one-dimensional graphs

  • Received: 28 October 2024 Revised: 28 April 2025 Accepted: 23 May 2025 Published: 12 June 2025
  • We propose numerical schemes for the approximate solution of problems defined on the edges of a one-dimensional graph. In particular, we consider linear transport and a drift-diffusion equations, and discretize them by extending finite volume schemes with upwind flux to domains presenting bifurcation nodes with an arbitrary number of incoming and outgoing edges, and implicit time discretization. We show that the discrete problems admit positive unique solutions, and we test the methods on the intricate geometry of electrical treeing.

    Citation: Beatrice Crippa, Anna Scotti, Andrea Villa. A monotone finite volume scheme for linear drift-diffusion and pure drift equations on one-dimensional graphs[J]. Networks and Heterogeneous Media, 2025, 20(2): 670-700. doi: 10.3934/nhm.2025029

    Related Papers:

    [1] Patrick Henning, Mario Ohlberger . The heterogeneous multiscale finite element method for advection-diffusion problems with rapidly oscillating coefficients and large expected drift. Networks and Heterogeneous Media, 2010, 5(4): 711-744. doi: 10.3934/nhm.2010.5.711
    [2] Dimitra Antonopoulou, Georgia Karali . A nonlinear partial differential equation for the volume preserving mean curvature flow. Networks and Heterogeneous Media, 2013, 8(1): 9-22. doi: 10.3934/nhm.2013.8.9
    [3] Stefan Berres, Ricardo Ruiz-Baier, Hartmut Schwandt, Elmer M. Tory . An adaptive finite-volume method for a model of two-phase pedestrian flow. Networks and Heterogeneous Media, 2011, 6(3): 401-423. doi: 10.3934/nhm.2011.6.401
    [4] Karoline Disser, Matthias Liero . On gradient structures for Markov chains and the passage to Wasserstein gradient flows. Networks and Heterogeneous Media, 2015, 10(2): 233-253. doi: 10.3934/nhm.2015.10.233
    [5] François James, Nicolas Vauchelet . One-dimensional aggregation equation after blow up: Existence, uniqueness and numerical simulation. Networks and Heterogeneous Media, 2016, 11(1): 163-180. doi: 10.3934/nhm.2016.11.163
    [6] Markus Musch, Ulrik Skre Fjordholm, Nils Henrik Risebro . Well-posedness theory for nonlinear scalar conservation laws on networks. Networks and Heterogeneous Media, 2022, 17(1): 101-128. doi: 10.3934/nhm.2021025
    [7] Caihong Gu, Yanbin Tang . Global solution to the Cauchy problem of fractional drift diffusion system with power-law nonlinearity. Networks and Heterogeneous Media, 2023, 18(1): 109-139. doi: 10.3934/nhm.2023005
    [8] Piermarco Cannarsa, Genni Fragnelli, Dario Rocchetti . Null controllability of degenerate parabolic operators with drift. Networks and Heterogeneous Media, 2007, 2(4): 695-715. doi: 10.3934/nhm.2007.2.695
    [9] Chunlin Du, Yu Zhang, Haolan Qu, Haowen Guo, Xinbo Zou . Acceleration of solving drift-diffusion equations enabled by estimation of initial value at nonequilibrium. Networks and Heterogeneous Media, 2024, 19(1): 456-474. doi: 10.3934/nhm.2024020
    [10] Panpan Xu, Yongbin Ge, Lin Zhang . High-order finite difference approximation of the Keller-Segel model with additional self- and cross-diffusion terms and a logistic source. Networks and Heterogeneous Media, 2023, 18(4): 1471-1492. doi: 10.3934/nhm.2023065
  • We propose numerical schemes for the approximate solution of problems defined on the edges of a one-dimensional graph. In particular, we consider linear transport and a drift-diffusion equations, and discretize them by extending finite volume schemes with upwind flux to domains presenting bifurcation nodes with an arbitrary number of incoming and outgoing edges, and implicit time discretization. We show that the discrete problems admit positive unique solutions, and we test the methods on the intricate geometry of electrical treeing.



    The present work is focused on the investigation of a finite volume-based numerical scheme for the solution of time-dependent linear advection and advection–diffusion partial differential equations (PDEs), defined on one-dimensional graphs. Our target application is the simulation of charges' movement inside defects embedded in the insulation systems of power networks. In fact, the reliability of an electric system is highly influenced by internal propagation of defects in its insulating components, a phenomenon known as treeing [1,2]. Electrical treeing structures are characterized by highly branched geometries [3,4], evolving over time due to partial discharges [5] in cold plasmas, where electrons and ions are free to move and cause avalanche effects [6]. While several semi empirical models exist in the literature [7,8], to the best of the authors' knowledge, models based on first principles remain limited to either simple configurations or early-stage treeing structures [9,10], as fully developed geometries are so complex that even the creation of a coherent mesh is almost impossible to achieve. For these reasons, the use of simplified, one-dimensional (1D) approximations of treeing as a graph is a more viable approach. Previously developed methods [11,12,13,14], introduced for the three-dimensional (3D) problem, will be adapted to simulate, on a 1D graph, the evolution of the volume concentration of charged particles inside the electrical treeing, modeled as drift-diffusion equations, and of the charge concentrations on the surface separating the treeing and the external dielectric domain, which instead are modeled as linear transport equations. A crucial property that the numerical methods for such models must satisfy is the positivity of the solution, which is necessary in view of the future coupling of the described problem with the chemical effects [2] to account for partial discharges and the electron avalanche effect. The one-dimensional domain can be described as a network, made of interconnected nodes, with differential equations defined on its edges and proper conditions at the intersections. A complete introduction on network structures can be found in [15,16].

    Further applications of problems on networks, with possibly non-linear drift, can be found in traffic flow on road networks [17,18], gas flow in pipe networks [19,20], supply chain models [21], water networks [22], and neuron models [23]. For these problems, ad-hoc solvers have been proposed [24,25,26,27], based on Godunov schemes [28], discrete velocities kinetic methods [29], finite differences [30] and domain decomposition [31]. Moreover, an extension of finite element methods to differential equations on quantum graphs [32] is discussed in [33].

    The particular structure of one-dimensional graphs introduces some complications at the nodes. In fact, the concentration and transport velocity are not defined there, since the functions on metric graphs are defined on the edges. Moreover, as junctions are represented by 0-dimensional nodes, the velocity direction is also ambiguous in such points. As a consequence, the definition of fluxes is not a trivial extension of fluxes on one-dimensional domains, as has been pointed out, for instance, in the case of nonlinear transport, where specific Riemann solvers [18] are introduced to correctly partition fluxes at junctions. A deep investigation and axiomatic definition of advection operators on graphs is presented in [34].

    In this paper, we focus on the solution of problems with a linear advection term, without any additional assumption on the advection velocities on each branch of the domain. Since, in the application of our interest, the equations represent the evolution of charge densities, we look for a conservative and monotone numerical scheme. For higher-dimensional problems, not defined on graphs, many numerical methods have been adapted to ensure such properties: Xu and Zikatanov [35] discussed the positivity of finite elements for convection-dominated equations, Lipnikov et al. [36] presented a monotone finite volume scheme on polygonal grids, while a conservative finite difference scheme is introduced by Hundsdorfer et al. [37], with flux being limited to ensure the positivity of the solution.

    Finally, we are interested in coupling, in the future, these one-dimensional equations with the electrostatic problem, defined both on the graph and on an external three-dimensional insulator where the electrical treeing is formed [38], requiring high computational effort for the numerical solution. Thus, we also look for an efficient numerical method for the solution of the linear transport and advection-diffusion equations. The coupling between the electrostatic equation, which describes the potential, and linear advection–diffusion equations, which model the density of charged species, is commonly known as the drift-diffusion problem [39,40,41]. This has been widely studied in the field of semiconductors and is typically solved via the Scharfetter–Gummel method [42], resulting in a modification of the advection–diffusion equations that introduces the exponential dependence of the diffusion coefficient on the electric field. However, in our case, the restrictive assumptions of constant mobility are not applicable. Similar mixed-dimensional coupled problems with transport on one-dimensional structures are also typical of flows in two-dimensional (2D) fractured porous media [43] and reservoir simulations [44], root–soil interaction [45], and drug delivery through microcirculation [46].

    We decided to investigate the properties of a numerical method that combines an extension of the finite volume scheme with upwind flux for space discretization and the implicit Euler method for the time discretization. This choice is motivated by the need for a conservative and monotone scheme to solve the model of densities of charged species, and by the efficiency of the low-order schemes, which do not introduce too much computational cost in the solution of the complete coupled problem. Moreover, the considered schemes provide stability when solving equations with a dominant transport term. A similar approach is applied in [43] to mixed-dimensional fluid flow and transport in porous media. Although this is a typical approach for the discretization of linear advection and diffusion problems, to the best of the authors' knowledge, there is no rigorous proof in the literature of the existence and uniqueness of the discrete solutions on one-dimensional graphs, or of monotonicity and consistency in the presence of bifurcation nodes. Instead, convergence results are discussed by Badwaik and Ruf [47] for upwind finite volumes applied to problems with jump discontinuities in the flux. We will show that the same convergence rate is preserved when these discontinuities are given by the jumps of the transport velocities at bifurcation nodes. The method presented in this paper can be regarded as an extension of these well established methods to graph domains: Indeed, in our case, the graph does not necessarily derive from a topological reduction of higher-dimensional thin domains. Therefore, we have to discuss the issue of defining fluxes at 0D nodes without relying on their physical meaning as small intersection regions.

    We start from a general formulation of linear transport and advection–diffusion problems on one-dimensional graphs, without introducing specific physics-based assumptions regarding the transport velocities. We only consider the flux conservation conditions at the internal nodes, which is necessary for the well-posedness of the problems [48,49]. We then semidiscretize the problem in time using an Implicit Euler scheme, and discuss an extension of the well-known finite volume method with two-point flux approximation (TPFA) and upwind flux in presence of bifurcation points or junctions between two segments with different direction. Finally, we test the discussed schemes on increasingly complex geometries, verifying that the well-known convergence properties of upwind and TPFA finite volumes on a 1D straight line are preserved in presence of bifurcation nodes, with the proposed cell-centered flux approximation. The last test case illustrates an example of the application of the discusses methods to the solution of the drift-diffusion problem on the electrical treeing, where the transport velocity is piecewise constant, consistent with the approximation of the electric field presented in [38].

    The original contributions of this paper consist of the analysis of the linear transport and advection–diffusion problems on a one-dimensional graph in a very general form, which introduces some complications in the extension of the widely discussed finite volume methods. In particular, to overcome the absence of a univocal definition of the flux at the junctions, due to the eventual discontinuities of transport velocities and concentrations and to the different directions of the velocities on the corresponding incoming and outgoing edges, we introduce cell-centered flux approximations and replace the unknown concentration at the internal nodes on the basis the flux conservation conditions, leading to the introduction of weights in the transport flux approximation. Based on the M-matrix structure of the transport discrete linear system and on the connectivity of the graph domain, reflected in the structure of the discrete linear system of diffusion, we prove the existence and uniqueness of the numerical solutions and the positivity of the scheme, meaning that, starting from a non-negative initial condition, at every time step, the solution is non-negative. Moreover, we also discuss the consistency of the numerical fluxes at the graph nodes, in the presence of bifurcations with an arbitrary number of incoming and outgoing edges.

    The rest of this manuscript is organized as follows: In Section 2, we present the structure of the domain and define the model problems on it; the numerical methods are introduced in Section 3. The existence, uniqueness, and positivity of their solution, together with consistency of the numerical fluxes, is investigated in Section 4. Finally, in Section 5, we present some tests of these methods on increasingly complex geometries, from a straight line to graphs with only one intersection node to an application to the geometry of a typical electrical tree.

    We will focus on time-dependent problems on a finite time interval [0,T], TR, and on a spatial domain represented by a one-dimensional, connected, simple, and directed graph Λ=(E,N), consisting of a set of edges E={ek}nek=1 and a set of nodes N={Υk}nΥk=0. Each edge is parametrized by its arc length parameter lk and normalized with a map πk : [0,1]ek. We suppose that the parametrization is opposite to the direction of the advective flux, in agreement with the previous literature, and we say that an edge ek is incoming to a node Υ if πk(0)=Υ and outgoing if πk(1)=Υ, namely if the flux on ek is directed towards Υ or away from it, respectively. We denote by E+Υ and EΥ the sets of incoming and outgoing edges of a node ΥN

    E+Υ:={ejE : πj(1)=Υ},
    EΥ:={ejE : πj(0)=Υ}.

    We use EΥ=E+ΥEΥ to denote the set of all the edges intersecting at Υ. In an analogous way, we can define the sets E+k and Ek, k=1,,ne, of incoming edges into the inflow endpoint, and outgoing edges from the outflow node of an edge ek, respectively

    E+k:={ejE : πj(0)=πk(1)},
    Ek:={ejEk : πj(1)=πk(0)}.

    On such a graph we can define the spatial derivative s on each edge as the derivative on each edge with respect to the parametrization.

    We consider the following pure transport problem:

    Find u : Λ×[0,T]R such that

    ut+s(cu)=f,on Λ×(0,T]; (2.1)

    with the transport velocity c: Λ[0,+) and the source term given by a function f: Λ×[0,T]R, and a drift-diffusion equation.

    Find u : Λ×[0,T]R such that

    uts(νus)+s(cu)=f,on Λ×(0,T], (2.2)

    with the coefficients c: Λ[0,+) and ν: Λ[0,+), and the source term f: Λ×[0,T]R. The set of nodes N can be partitioned into two non overlapping sets, namely N=NiNb, NiNb=, where we call Ni the set of internal nodes, i.e., nodes shared by at least two segments, and Nb the set of boundary nodes (i.e., connected to only one edge of the graph).

    On all the nodes in Nb, we must impose boundary conditions for problems (2.1) and (2.2). For each problem, we use Ns to denote the set of Dirichlet nodes, and by Ne the set of end nodes, where we impose Neumann conditions

    {u=ˉu,on every ΥNs,duds=0,on every ΥNe. (2.3)

    We assume that the nodes where we impose Dirichlet conditions are sources, meaning they have no incoming edges, and those with Neumann (or Robin) conditions are sinks, with no outgoing edges. In the advection problem (2.1), we can also have sink nodes where no boundary condition is imposed, if they are outflow nodes, where outflow conditions are not necessary.

    Moreover, in both cases, we allow one node to be shared among more than two edges, and we call such a node a bifurcation. We denote the set of bifurcations as NbifNi. In each internal node, there must be at least one incoming and one outgoing edge, in order to avoid mass accumulation there.

    We use fk to denote the restriction f|ek of a function f: ΛR to the edge ek, k=1,,ne. In particular, we can write the velocity and concentration on each edge as ck=c|ek and uk=u|ek, respectively. Then, the conservation of fluxes at a node ΥN is:

    ekE+Υukck=ejEΥujcj,

    and the complete formulation of the transport problem (2.1) on the network reads as follows:

    {ukt+s(ckuk)=fk,k{1,,ne},in[0,1]×(0,T],(2.4a)ekE+Υckuk=ejEΥcjuj,ΥNi,in(0,T],(2.4b)uk(πk(1),t)=ˉuk(t),k{1,,ne}:πk(1)Ns,t(0,T],(2.4c)uk(x,0)=u0(x),x[0,1],k{1,,ne},(2.4d)

    while the complete formulation of the drift-diffusion Eq (2.2) problem on the network is:

    {ukts(νkuks)+s(ckuk)=fk,k{1,,ne},in[0,1]×(0,T],(2.5a)ekE+Υckuk=ejEΥcjuj,ΥNi,in(0,T],(2.5b)uk(πk(1),t)=ˉuk(t),k{1,,ne}:πk(1)Ns,t(0,T],(2.5c)uks(πk(0),t)=0,k{1,,ne}:πk(0)Ne,t(0,T],(2.5d)uk(x,0)=u0(x),x[0,1],k{1,,ne}.(2.5e)

    All the hypotheses introduced above are necessary for stating the well-posedness of these continuous Eqs on a graph. Indeed, according to [48], the drift-diffusion problems (2.5a)–(2.5e) admits a unique solution if the coefficients ν and c belong to the space Cα,α/2(Λ×[0,T]) of Holder continuous functions for some α, while the existence of a unique solution of the transport problems (2.4a)–(2.4d) is ensured by Banasiak and Namayanja [49].

    Remark 1. In [49], the possibility of dealing with homogeneous Neumann conditions on an open graph is discussed. In this case, we would have the injection of some quantity from the shource nodes, whose concentration is represented by the unknown u, and no outflow from the sinks, so we need a workaround in order to make the problem well-posed and avoid mass accumulation on the vertices. The authors propose to extend the graph by adding one extra vertex for each sink, connected to it by two edges going in opposite directions, as to create a cycle (see Figure 1), and they prove that these additional vertices do not influence the behavior of the solution on the original graph, the flow on each appended subgraph is asymototically periodic in time, and the edges incoming to the sinks will be eventually depleted as t, or in a finite time if they are in the acyclic part of the graph. However, as we will see in Section 4, this solution is not required in the discrete setting.

    Figure 1.  A simple one-dimensional domain with three edges and one bifurcation node Υ2. If we impose homogeneous Neumann conditions on the set of sinks Ne={Υ3, Υ4}, the graph should be extended by also including the nodes Υ5 and Υ6 and the additional dashed edges, according to [49].

    We start by discretizing the two model problems (2.1) and (2.2) in time by introducing a set of equispaced instants {tl : l=0,,N} on [0,T], such that t0=0, tN=T and tl+1=tl+Δt, l{0,,N1}. We denote by ul the solution u at time tl, l=0,,N and apply an implicit Euler discretization scheme on both problems

    ul+1ulΔt+dds(ul+1 c)=fl+1,on Λ, l{0,,N1}, (3.1a)
    ul+1ulΔtdds(νdul+1ds)+dds(ul+1 c)=fl+1,on Λ, l{0,,N1}. (3.1b)

    The choice of an implicit solver is due to the final aim of defining methods for the simulation of thermal plasma without restrictions on the time step length, due to stability issues [12,50].

    We complete the discretization of Eqs (3.1a) and (3.1b) by applying a finite volume (FV) method on the 1D domain Λ, represented by a set of connected segments with possible branches.

    On the basis of the concepts discussed in the previous section, we introduce a numerical scheme for the approximate solution of the aforementioned PDEs. We start by the discretization of a linear transport problem in Section 3.1.1, followed by the diffusion equations in Section 3.1.2, which introduces the solver for drift-diffusion problems, presented in Section 3.1.3.

    Consider time-dependent transport equations in the form of Eq (2.1), semi-discretized in time as Eq (3.1a). For simplicity, as a partition for Λ, we consider the set of segments representing its edges {ek}nek=1, and we complete the discretization of the equation by integrating Eq (3.1a) on each segment:

    ekul+1ulΔt=ekdds(ul+1c),k{1,,ne}, l{0,,N1}. (3.2)

    Notice that a further partition of the edges can be trivially introduced by considering more internal nodes on each edge.

    On each segment, we approximate the value of u as a constant, given by its integral mean. From now on, with a little abuse of notation, we will use ulk to denote the constant approximation of the solution on each segment ek at time tl: ul|ekulk, k{1,,ne}, l{0,,N}. We choose to approximate the right-hand integral as the sum of numerical fluxes at the endpoints of each segment, with an upwind scheme. Integrating the first derivative of the flux F(u)=uc on the edge ek, we obtain the difference between its value on the end nodes Υ2k=πk(0) and Υ1k=πk(1). Since the numerical solution may have jump discontinuities on the nodes separating two neighboring edges, we have to define a numerical flux ˜F to set its value there. According to the upwind method, we set the flux on the first end node Υ1k to be equal to the upstream flux ˜F(uΥ)=ckuΥ, with the concentration u evaluated in the incoming neighbors of ek, the second Υ2k to be equal to the upstream flux ˜F(uk)=ckuk on ek, thus obtaining the following approximation:

    ekdds(ul+1c)ck(ul+1kul+1Υ),k{1,,ne}, l{0,,N}, (3.3)

    where ul+1Υ indicates the concentration u at time tl+1 from the incoming neighbors of ek through the vertex Υ=πk(1). We define it as a function of the concentrations on the neighboring edges of ek, based on the flux conservation condition at the junction node Υ.

    We can substitute Eq (3.3) into Eq (3.2) and ul with its constant approximation ulk on each segment to obtain the following expression:

    |ek|Δtul+1k+ck(ul+1kul+1Υ)=|ek|Δtulk,k{1,,ne},l{0,,N}. (3.4)

    On the node Υ=πk(1), we can impose the conservation of fluxes and determine the portion of mass transported from the incoming edges ei to the outgoing edges ek. The numerical fluxes ˜F+Υ and ˜FΥ, entering and leaving the node are given by the sum of fluxes on its incoming and outgoing edges, respectively:

    ˜F+Υ=ejE+Υ˜Fj=ejE+Υcjuj,˜FΥ=ejEΥ˜Fj=ejEΥcjuΥ. (3.5)

    If we balance the fluxes Eq (3.5), we obtain:

    ul+1Υ=ejE+Υ(cjekEΥckul+1j). (3.6)

    Then, the value of u at a node Υ is given by a weighted sum of its value on the incoming edges, where the weights are given by

    wj=cjelEΥcl,j : ejE+Υ, ΥNNe.

    A similar weighted flux is also introduced in other works, such as [51], which discusses finite volumes on one-dimensional domains with bifurcation nodes.

    If we substitute the expression of ul+1Υ obtained by this relation in Eq (3.4), we end up with the following numerical scheme:

    (|ek|Δt+ck)ul+1kckeiE+kwiul+1i=|ek|Δtulk,k{1,,ne}, l{0,,N}, (3.7)

    which can be written in matrix form as:

    Mtrul+1=gl,l{0,,N1}, (3.8)

    where gl=[|e1|Δtul1,,|ene|Δtulne]TRne is a vector depending on the unknowns at the previous time tl, and MRne×ne is a square matrix with entries given by:

    (Mtr)ki={|ek|Δt+ck,if i=k, k=1,,ne,ckwi,if ik, i,k=1,,ne, eiE+k,0,otherwise. (3.9)

    Let us introduce a numerical scheme for the diffusion term, which will then be combined with the method presented in the previous section.

    Consider the following steady equation:

    dds(νduds)=0,on Λ, (3.10)

    completed by the Dirichlet and Neumann boundary conditions (2.3) introduced in Section 2.

    For clarity of exposition, we consider the simple case where Λ is made of three segments, one having vertices numbered as Υ0 and Υ1, and the other two branching from vertex Υ1, as in Figure 2. We denote A as the vertex with index 0, I as the vertex with index 1 and B and C as the endpoints of the two edges of the branch. Then, the segments composing Λ are e1=AI, e2=IB, and e3=IC.

    Figure 2.  Simple one-dimensional domain with three edges and one bifurcation node Υ1. The set of edges is E={e1=AI,e2=IB,e3=IC}. We set Dirichlet boundary conditions on A and Neumann boundary conditions on the set of end nodes Ne={B,C}. The only bifurcation node is I.

    On this domain, the problems (3.10) and (2.3) read as follows:

    {dds(νduds)=0,on Λ,u(A)=ˉu,duds(B)=duds(C)=0.

    If we now integrate the first equation of this system on each segment applying the fundamental theorem of calculus, we obtain:

    ekdds(νduds)=ek0,ν(πk(0))duds(πk(0))ν(πk(1))duds(πk(1))=0.

    Let us denote νk,j as the value of the diffusion coefficient defined on edge ej, evaluated at the node Υ=πk(j), j=0,1. Furthermore, we can approximate the first derivative of u at a node by finite differences, similar to the TPFA approach [52]. In particular, considering two adjacent and collinear segments, as in Figure 3, we approximate its value on the node between ek and ek+1 as

    dudsTk,k+1(uk+1uk),k=1,,ne1,
    Figure 3.  Segment discretized by contiguous elements without bifurcations. The set of edges is E={e1=AB, e2=BC} and there are no bifurcations. We impose Dirichlet boundary conditions on the source node Υ0 and homogeneous Neumann conditions on the sink Υ2.

    where Tk,k+1 denotes the inverse of the distance between the centers of ek and ek+1, and in this case, it is simply

    Tk,k+1=(|ek|+|ek+1|2)1. (3.11)

    If the two segments are not collinear, the approximation of the first derivative remains the same and the distance T1k,j between the centers of consecutive elements is considered as the distance on the graph, as in Eq (3.11). If a node belongs to one edge only, then it must be either a source node ΥNs or an end node ΥNe. In the first case, we have

    duds(Υ)Tk,1(ukˉuk),ΥNs, (3.12)

    where we denote Tk,1=2|ek| as the inverse of distance from the center of the segment ek to the source node. Otherwise, if it is an end node, we impose homogeneous Neumann boundary conditions and we can directly substitute

    duds(Υ)=0,ΥNe. (3.13)

    Finally, in our case, a node may be shared by more than two edges. Hence, we will introduce the following generalization.

    Consider the node I in the graph of Figure 2, connected to three edges. Here, we separately consider the three contributions and then balance the flux at the intersection Υ1=I. If we assume that the orientations of the edges e1, e2, and e3 are AI, IB, and IC, respectively, then Υ1=π1(0)=π2(1)=π3(1). In this case, the only incoming edge is e1, while the outgoing ones are e2 and e3. If we assume that there is no mass accumulation at the nodes, then the fluxes are balanced, meaning that the sum of incoming and outgoing fluxes at this point is null

    ν1(Υ1)du1ds(Υ1)=ν2(Υ1)du2ds(Υ1)+ν3(Υ1)du3ds(Υ1). (3.14)

    We can consider the bifurcation nodes as boundary points of subgraphs of Λ not containing bifurcations. On these additional boundary points, we impose Dirichlet boundary conditions:

    uk(Υ)=u(Υ),ΥNbif. (3.15)

    Clearly, we do not know the value of the function u at these nodes, but we can compute it by making use of the flux balance condition (3.14). This means that we need to consider an additional set of discrete unknowns {une+j}nbifj=1, where nbif=|Nbif| is the number of bifurcation nodes in Λ, such that une+j=u(Υj), j=1,,nbif.

    The approximation of derivatives on ΥNbif is given by:

    dukds(Υj){Tk,ne+j(ukune+j),if Υj=πk(1),Tk,ne+j(une+juk),if Υj=πk(0). (3.16)

    Here, we denote Tk,ne+j=2|ek| as the inverse of the distance between the center of the element ek and the bifurcation node Υj along the graph.

    Combining Eqs (3.12), (3.13), and (3.16), we obtain a linear system of the form

    Adiffu=0, (3.17)

    where Adiff is an (ne+nbif)×(ne+nbif) matrix with the following block structure:

    Adiff=[BCCTD], (3.18)

    with BRne×ne, CRne×nbif, and DRnbif×nbif diagonal matrix. The entries of these matrices are:

    Ck,j={νk,j|ek|Tk,j,if πk(1)=ΥjNbif or πk(0)=ΥjNbif,0,otherwise,k=1,,ne,j=1,,nbif;Dk,j={nei=1Ci,j,if k=j,0,otherwise,k,j=1,,nbif;Bk,j={νk,j|ek|Tk,j,if ejEk,v= ejekNbif,0,otherwisek,j=1,,ne,jk;Bk,k={nej=1, jkBk,jnbifj=ne+1Ck,j,if πk(1)Ns,nej=1, jkBk,jnbifj=ne+1Ck,j+νk,1Tk,1,otherwisek=1,,ne. (3.19)

    The vector uRne+nbif introduced in Eq (3.17) contains the unknowns, ordered as the ne values of u at the edges ek, k=1,,ne, ˆu=[u1,,une]T, and at the bifurcation points uB=[une+1,,une+nbif]T

    u=[ˆuuB].

    By combining Eqs (3.8) and (3.17) we can discretize the drift-diffusion Eq (2.2), which contains a system of ordinary differential equations in the unknowns ˆu and uB. Note that in the discretization of the transport term, the only considered unknown is ˆu, representing the approximate values of u on the edges. Therefore, we need to augment the corresponding matrix in order to be able to sum it to the diffusion part

    M=[I+Mtr000]+ΔtAdiff, (3.20)

    where MtrRne×ne and AdiffR(ne+nbif)×(ne+nbif) are defined in Eqs (3.9) and (3.18), respectively, while the vector of unknowns and the right-hand side are defined as in Section 3.1.2, Eq (3.17).

    In this section, we focus on proving some properties of the matrices defined in Eqs (3.9), (3.18), and (3.20), proving the existence and uniqueness of the solutions and the positivity of the associated linear systems. We will show these results on the basis of some properties of Z-matrices and M-matrices [53], recalled for the readers' convenience

    Definition 1. A matrix ARN×N is a Z-matrix if AZN×N, with

    ZN×N:={MRN×N : Mij0 i,j suchthat ij, 1i,jN}.

    Definition 2. A matrix MRN×N is an M-matrix if  sR, BRN×N such that M=sIB, with Bij0, 1i,jN, and sρ(B), where ρ(B) denotes the spectral radius of B.

    Let us start by analyzing the discrete transport equation introduced in Section 3.1.1. Observe that the matrix Mtr, defined in Eq (3.9), has all non-negative diagonal entries and non-positive extra-diagonal entries. As a consequence, the following property holds.

    Property 1. The matrix Mtr defined in Eq (3.9) is a Z-matrix.

    Moreover, in the following, we prove that Mtr has all positive column sums. These properties allow us to prove that it is invertible and that, in Eq (3.8), g0u0. Consequently, the system Mtru=g admits a unique solution, which is positive if the source term is. In fact, we have

    Property 2. The matrix Mtr defined in Eq (3.9) has positive column sums.

    Proof. If we fix any column i{1,,ne} of Mtr, the sum of its elements is given by

    nek=1(Mtr)ki=(Mtr)ii+nek=1ki(Mtr)ki=|ei|Δt+cikEΥckwi=|ei|Δt+cicikEΥckkEΥck=|ei|Δt>0,

    where Υ=πk(0).

    This result allows us to prove the following theorems:

    Theorem 1. Let Mtr be the ne×ne matrix defined by Eq (3.9). Then, Mtr is invertible and the system Mtru=g is positive.

    Proof. Since Mtr is a Z-matrix (Property 2), so is its transpose MTtr is. According to Theorem 1 of [53], this property is necessary and sufficient to state that MTtr is a nonsingular M-matrix, which is also equivalent to MTtr being inverse-positive, i.e., (MTtr)1, ((MTtr)1)ij0, and M1tr0.

    By the definition of inverse, we have I=MTtr(MTtr)1, where I denotes the ne×ne identity matrix. If we transpose both sides of the equations, we get

    I=(MTtr(MTtr)1)T=((MTtr)1)TMtr.

    This is only possible if Mtr is invertible and its inverse is M1tr=((MTtr)1)T. Then, Mtr is also inverse-positive, i.e., (M1tr)ij0,i,j{1,,ne}, and, as a consequence,

    uk=nei=1(M1tr)kifi0,k{1,,ne},

    because they are linear combinations of non-negative quantities.

    Finally, we can prove the consistency of the proposed numerical flux on the graph nodes.

    Theorem 2. The flux approximation Eq (3.3) is consistent at every internal graph node.

    Proof. Denote the exact flux at an internal node ΥNi as FΥ=(cu)(sΥ), where sΥΛ is the position of the node on the graph and consider the corresponding incoming and outgoing numerical fluxes ˜F+Υ and ˜FΥ defined in Eq (3.5).

    We introduce the evaluation of the numerical fluxes with respect to the exact solution as follows:

    ˜F,+Υ=eiE+Υ(cu)(si),
    ˜F,Υ=eiEΥc(si)u(sΥ),

    and the approximate flux at the bifurcation node should satisfy FΥ=F,+Υ=F,Υ, by balance of fluxes.

    We then obtain the following expression for the exact solution at the node:

    u(sΥ)=eiE+Υc(si)u(si)eiEΥc(si),

    which can be substituted in the definition of the exact flux, to obtain:

    FΥ=c(sΥ)eiE+Υc(si)u(si)eiEΥc(si)=c(sΥ)eiEΥc(si)FΥ.

    Finally, the truncation error can be computed as follows:

    τΥ=FΥFΥ=(c(sΥ)eiEΥc(si)1)FΥ=(c(sΥ)eiEΥc(si))u(si).

    We have no hypotheses on the continuity of the speed c at nodes, so it may present jumps at those points. We need to introduce the following compatibility condition, setting the value of the velocity at the nodes as:

    c(sΥ)=eiEΥc(si), (4.1)

    where si indicates the position of the midpoint of the edge ei on Λ. Notice that, if the node Υ has only one outgoing edge, the velocity at the node is set equal to the outgoing velocity.

    By this definition of the speed at the nodes, we get τΥ=0 and we can conclude that the the method is strongly consistent at the graph nodes.

    Remark 2. The proposed numerical scheme does not take the outflow conditions on sinks into account, but it can be adapted to enforce boundary conditions, such as Robin or homogeneous Neumann conditions, by imposing them on the edges connected to the sinks. The approximation error on these edges should decrease with the dimension of the space discretization, and consistency should still hold.

    Let us now analyze the diffusion matrix Adiff, defined by Eq (3.18). This is a block (ne+nbif)×(ne+nbif) matrix, with the following properties:

    Property 3. The matrix Adiff defined in Eq (3.18) is a Z-matrix;

    Property 4. The matrix Adiff defined in Eq (3.18) is symmetric.

    These two properties follow from the definition of the matrix. Moreover, we can prove the following result.

    Property 5. The row sums of the elements of the matrix Adiff defined in Eq (3.18) are always zero but for a number n=|Ns|>0 of rows, and the same holds for the columns.

    Proof. Let us start with the last nbif rows, consisting of the blocks

    [CTD],

    where the only non-null element of D on every row is on its diagonal and is equal, by definition (Eq (3.19)), to the opposite of the column sum of C. In this case,

    ne+nbifk=1(Adiff)jk=nek=1CTjk+nbifk=1Djk=nek=1Ckj+Dkk=nek=1Ckjnek=1Ckj=0,j=ne+1,,ne+nbif.

    Among the first ne rows, we find the non-null row sums, corresponding to the elements on whose end points the Dirichlet boundary conditions are imposed. Indeed, by definition (Eq (3.19)), the diagonal elements of Adiff are defined as the opposite of the sum of the non diagonal entries of each row, plus an additive term in the case of Dirichlet conditions imposed on an end node:

    ne+nbifk=1(Adiff)jk=(Adiff)kk+ne+nbifk=1kj(Adiff)jk=Bkk+ne+nbifk=1kj(Adiff)jk=={ne+nbifk=1kj(Adiff)jk+ne+nbifk=1kj(Adiff)jk=0,if πk(1)Ns,ne+nbifk=1kj(Adiff)jk+νk,1Tk,1+ne+nbifk=1kj(Adiff)jk=νk,1Tk,1>0,otherwise.

    Since the only rows with positive sums are the ones corresponding to the edges with Dirichlet boundary conditions on one end point, the number of rows of Adiff with a non-null sum is equivalent to the number of such edges |Ns|, which was assumed to be non-null.

    Finally, since Adiff is symmetric (Property 4), the same holds for the columns.

    A consequence of Property 5 concerns the diagonal dominance of the matrix Adiff.

    Property 6. The matrix Adiff defined in Eq (3.18) is diagonally dominant on every row and strictly diagonally dominant on a number n=|Ns| of rows.

    These rows correspond to edges where Dirichlet boundary conditions are imposed.

    In order to prove the existence and uniqueness of the solution and the positivity of the discrete diffusion problem, we also need to introduce the SC property, and show that it is satisfied by Adiff.

    Property 7. The (ne+nbif)×(ne+nbif) square matrix Adiff defined in Eq (3.18) satisfies the SC property (Definition 6.2.7 of [54]).

    p,q{1,,ne+nbif}, with pq, there is a sequence of distinct integers {ki}mi=1{1,,ne+nbif}, such that k1=p, km=q and (Adiff)k1k2, (Adiff)k2k3, , (Adiff)km1km are all nonzeros.

    Proof. Let us analyze the non-zero entries of the matrix: (Adiff)ik0k=i or where k and i represent neighboring edges, or where k and i represent a bifurcation and one of the edges connected to it.

    As a consequence, wehave k1,,km{1,,ne+nbif} distinct indices such that (Adiff)k1k2, (Adiff)k2k3,, (Adiff)km1km 0  a path on the 1D domain Λ connecting the edge or bifurcation k1 to the edge or bifurcation km without passing more than once through the same edge or bifurcation. This condition is true for every couple of edges/bifurcation points because the domain Λ is a completely connected graph by construction.

    Thus, we can conclude that Adiff satisfies the SC property.

    Thanks to these properties, we can prove the positivity of the numerical scheme proposed in Section 3.1.2 and the existence and uniqueness of the solution to the discrete problem.

    Theorem 3. Let Adiff be the (ne+nbif)×(ne+nbif) matrix defined in Eq (3.19). Then, Adiff is invertible and the problem (3.17) is positive.

    Proof. Since the matrix Adiff satisfies the SC Property 7, and Property 6 ensures that it is also diagonally dominant on each row and strictly diagonally dominant on at least one row, then it is invertible (see Corollary 6.2.9 of [54]).

    Moreover, these conditions are also satisfied by Adiff+˜D, where ˜D denotes a general diagonal (ne+nbif)×(ne+nbif) matrix ˜D with all positive diagonal entries. Thus, Adiff+˜D is nonsingular ˜D, by Corollary 6.2.9 of [54]. Consequently, since Adiff is a Z-matrix (Property 3), this implies (Theorem 1 of [53]) that Adiff is a nonsingular M-matrix and also inverse-positive.

    This means that (A1diff)ik0 i,k=1,,ne+nbif, and A1diff0. Then, u=A1diffg0, g0.

    Let us now explicitly write the system for the vectors ˆu and uB of unknowns on the edges and on the bifurcation points, respectively, starting from the discrete problem (3.17) in matrix form

    {Bˆu+CuB=f,CTˆu+DuB=0.

    Since D is invertible (a diagonal matrix with all positive entries), we can solve the second equation for uB and obtain the following system:

    {(BCD1CT)ˆu=f,uB=D1CTˆu.

    Observe that BCD1CT is nothing but the Schur complement of D, which we denote by A/D:

    {(A/D)ˆu=f,uB=D1CTˆu.

    Since A and D are nonsingular, we know that A/D is nonsingular as well (Theorem 1.2 of [55]) and that det(A/D)=det(A)det(D), by Schur's formula.

    Our goal is now to prove that also the Schur complement A/D is invertible and positive. To this end, we need to recall the definition of the inertia of an Hermitian matrix and an associated result (see [55]).

    Definition 3. We call the inertia of an Hermitian matrix A the triplet In(A):=(p(A),q(A),z(A)), where p(A) is the number of positive eigenvalues of A, q(A) is the number of negative eigenvalues of A, and z(A) is the multiplicity of the 0 eigenvalue.

    Theorem 4. Theorem 1.6 of [55]:

    Let A be an Hermitian matrix and A11 a nonsingular principal submatrix of A. In this case,

    In(A)=In(A11)+In(A/A11).

    We can now prove the following result.

    Theorem 5. Let A be the (ne+nbif)×(ne+nbif) matrix defined in Eq (3.18) and A/D be the Schur complement of its south east nbif×nbif block D. Then, A/D is invertible and the system (A/D)u=f is positive.

    Proof. The matrices A and D are both positive definite, because they have positive inverse and a diagonal with strictly positive entries. Then, In(A)=(ne×nbif,0,0), In(D)=(nbif,0,0), and In(A/D)=In(A)In(D)=(ne,0,0), by the previous theorem. Then, A/D only has positive eigenvalues, and therefore, it is positive definite.

    Let us now inspect the entries of the matrix CD1CT.

    The non-diagonal ones are positive because they are given the sums of non-negative quantities:

    (CD1CT)ik=nbifh=1CihCkhDhh0 i,k=1,,ne, ik,

    and the same holds for diagonal ones:

    (CD1CT)ii=nbifh=1C2ihDhh{>0if i is connected to a bifurcation,0otherwise.

    Indeed, (D1CT)ik=1Dii(CT)ik=CkiDii0 i=1,,nbif, k=1,,ne and Cih0 i=1,,nbif, k=1,,ne.

    In this case,

    (BCD1CT)ik=Bik(CD1CT)ik0,i,k=1,,ne, ik,

    because Bik0 and (CD1CT)0, i,k=1,,ne, ik, and

    (BCD1CT)ii=Bii(CD1CT)iinej=1jiBi,jnbifj=1Ci,jnbifj=1C2i,jDjj==nej=1jiBi,jnbifj=1CijDjjCi,jDjj>0,i=1,,ne,

    because nej=1, jiBi,j0, i{1,,ne}, DjjCi,j>0, Cij<0, Djj>0, i{1,,ne}, and j{1,,nbif}.

    Consequently, A/D=BCD1CT is a Z-matrix and has positive diagonal entries. Thus, the following conditions are equivalent (Theorem 1 of [53]):

    A/D has all positive eigenvalues;

    A/D is a nonsingular M-matrix;

    A/D is inverse-positive.

    As a consequence, A/D is invertible. Since it is inverse-positive, u=(A/D)1f0, i.e., the system (A/D)u=f is positive.

    Finally, we can prove that the flux approximation is consistent at the graph nodes by following a similar procedure as in Theorem 7.1 of [52] for the upwind flux with TPFA on the diffusion term.

    Theorem 6. If the diffusion coefficient is globally continuous on Λ and νkC1(ek), k{1,,Ne}, then the flux approximation Eq (3.1.2) is consistent at every internal node of Λ.

    Proof. Define the numerical flux entering a node Υ from the incoming edge ek as

    ˜F+Υ,k=νk,Υ2|ek|(uΥuk), (4.2)

    where uΥ and uk are the approximate solutions at the node Υ and on ek, respectively, and νk,Υ is the diffusion constant evaluated at the node.

    Similarly, let

    ˜FΥ,k=νk,Υ2|ek|(ukuΥ)

    be the numerical flux going from the node Υ into the outgoing edge ek.

    Then the total incoming and outgoing fluxes at the node Υ are ˜F+Υ=ekE+Υ˜F+Υ,k and ˜FΥ=ekEΥ˜FΥ,k.

    Let us now introduce the exact flux of the velocity on an edge ek evaluated at the node on sΥΛ as follows:

    FΥ,k=(νdu|ekds)(sΥ),

    and the evaluation of the numerical flux with respect to the exact solution u is

    FΥ=ekE+ΥF,+k=ekEΥF,k, (4.3)

    where F,+k and F,k are the evaluations of the numerical fluxes on ek for incoming and outgoing edges, respectively, with respect to the exact solution

    F,+k=νk,Υ2|ek|(u|ek(sΥ)u|ek(sk)), (4.4)
    F,k=νk,Υ2|ek|(u|ek(sk)u|ek(sΥ)).

    We want to show that the truncation error τΥ:FΥFΥ is controlled by the maximum length of the segments discretizing Λ, h=maxkE|ek|, i.e.,

    CR+ such that |τΥ|Ch. (4.5)

    We start by observing that since νC1(ek), ekEΥ, then uC2(ek); therefore, we can write the second-order Taylor expansion of u on each edge as u(sk)=u(sΥ)+|ek|2duds(sk)+Rk, with RkCk|ek|Ckh, Ck>0. We then substitute it into Eqs (4.4) and (4.2), with k{1,,Ne} CkR+ such that

    {F,+k=ωkFΥ+R+k, |R+k|Ckh, ekE+Υ,F,j=ωjFΥ+Rj, |Rj|Cjh, ejEΥ, (4.6)

    where ωkFΥ denotes the fraction of the flux FΥ through the node Υ coming from each edge ek, and is given by:

    ωk=νk,Υν(sΥ), ekEΥ.

    As a consequence of Eqs (4.3) and (4.6), we have

    ekE+kωkFΥ+R+=ekEkωkFΥ+R,

    with R+=ekE+ΥR+k and R=ekEΥRk.

    Thus, we obtain

    FΥ=RR+kEΥωk=RR+kEΥνk,Υν(sΥ).

    Assuming that the diffusion coefficient ν is continuous, its evaluation on each edge ek at the common end point must be the same as at the node, i.e., νk,Υ=ν(sΥ), ekEΥ. Then, denoting |EΥ| as the number of segments intersecting at the node Υ, we have ekEΥνk,Υ=|EΥ|ν(sΥ) and ekEΥωk=|EΥ|. Consequently, FΥ=RR+|EΥ|, and the truncation error is given by

    |τΥ|=|FΥFΥ|=|(|E+Υ|1)FΥ+R+|=|(|E+Υ|1)RR+|EΥ|+R+|.

    By the triangular inequality and recalling that every internal node has at least one incoming and one outgoing edge, i.e., |EΥ||E+Υ|1

    |τΥ|(|E+Υ||EΥ|+1|EΥ|)(|R|+|R+|)+|R+|2(|R|+|R+|)+|R+|(2C+3C+)h=Ch,

    thus proving that Eq (4.5) holds.

    In Section 3.1.3, we introduced a numerical scheme for advection–diffusion Eq (2.2) and obtained a linear system of equations Mu=f. The matrix M is given by a sum of block matrices resulting from discretization of the pure transport and diffusion problems, analyzed above. Then, on the basis of the results presented for the two separate problems, we can prove the same results for this problem as well. In particular, since M is given by the sum of the two Z-matrices Mtr and Mdiff,

    Property 8. The matrix M defined in Eq (3.20) is a Z-matrix.

    Thanks to the construction of M, its resulting structure is very similar to those of Mtr and Adiff, and thus we can prove the desired properties similarly to what we have done for the two separate problems.

    Since the set of nonzero elements in Mtr is a subset of the non-zero elements of the top left block B of Adiff, the sum of these two matrices still satisfies the SC property. Moreover, transposing them does not affect their pattern and, as a consequence, also MT satisfies the SC property.

    Property 9. The matrix M defined in Eq (3.20) and its transpose MT satisfy the SC property (Definition 6.2.7 of [54]).

    Moreover, as a consequence of Properties 2, 5, and 4, we know that M has non-negative column sums.

    Property 10. The row sums of the elements of the matrix M defined in Eq (3.20) are always zero but for a number n=|Ns|>0 of rows.

    As a consequence of Properties 7, 9, and 10, we can show, as in the proof of Theorem 3, that the transpose MT of the matrix M is inverse-positive and, therefore, so is M. Then, the following result holds:

    Theorem 7. The matrix M defined in Eq (3.20) is non-singular and the associated discrete drift-diffusion scheme is positive.

    Finally, we can apply the same procedure as that adopted for the diffusion equation in the previous section, and explicitly write the system (3.20) for the unknown ˆu as follows:

    {(M/Dt)ˆu=f,uB=D1CTˆu,

    where we have used Dt to denote the bottom right block of the drift-diffusion matrix, Dt=ΔtD, and M/Dt as its Schur complement. Following the same steps as in the proof of Theorem 5, we can prove the following result:

    Theorem 8. Let M be the (ne+ndiff)×(ne+ndiff) matrix defined by Eq (3.20) and M/Dt is the Schur complement of its south-east diagonal block Dt. Then, M/Dt is invertible and the system (3.20) is positive.

    In this section, we show some applications of the proposed methods to different geometries. We first check that the approximation errors on one-dimensional straight lines respect the well-known convergence results of finite volumes, then we study the convergence on domains with bifurcations and finally apply the methods to a drift-diffusion problem on the complex geometry of an electrical tree.

    The error is computed as the L1 norm of the difference between the approximated and the exact solutions at the final time T, normalized with respect to the L1 norm of the exact solution:

    error=uex(,T)˜uNtL1([0,1])uex(,T)L1([0,1]),

    where ˜ul denotes the approximate solution at each discrete time tl as a piecewise constant function, taking values ˜ulk on ek, k=1,,ne, and l=1,,Nt.

    The discretization presented in Section 3.1.1 of the transport equation on a line is nothing but the finite volume upwind scheme; see [56].

    We apply it to solve the homogeneous transport problem in Eq (2.1) with the null source term f=0 and an advection speed c=0.5, on the space-time domain [0,1]×[0,1], with the initial condition u0(x)=0, x[0,1], and the inflow boundary condition u(0,t)=1, t(0,1].

    The exact solution of this problem consists of the propagation of a rectangular wave of speed c

    uex(x,t)={1,x[0,ct], t(0,1],0,x(ct,1], t(0,1]. (TC1A)

    In the plots in Figure 4, we can observe that the convergence rate in space and time of 12, which is the expected value for drift equations with discontinuities [47], is recovered by our method.

    Figure 4.  TC1A – Convergence test for the transport equation on a line with a discontinuous exact solution. The normalized error in norm L1([0,1]) is computed at time t=1 corresponding to space and time meshes with different dimensions.

    Convergence of order 1 is expected when the solution is continuous, and we can appreciate it in Figure 5, where we have applied the proposed method to an homogeneous transport equation with the advection speed c=0.5 with the initial condition

    u0(x)=sin(πx),x[0,1],
    Figure 5.  TC1B – Convergence test for the transport equation on a line with a continuous exact solution. The normalized error in norm L1([0,1]) is computed at time t=1 corresponding to space and time meshes with different dimensions.

    whose exact solution is

    uex(x,t)=u0(xct),x[0,1], t[0,1]. (TC1B)

    We have imposed a Dirichlet condition on the inflow boundary accordingly: u(0,t)=u0(ct), t[0,1].

    In both cases, we observe that the error in space saturates for small values of h if the time grid is not refined enough and vice versa. This is due to the predominance of the error made on the coarsest grid. In the discontinuous (TC1A) case (Figure 4), the two plots have almost the same behavior and both errors saturate for comparable mesh spacings, while in the continuous case (Figure 5), the error in space needs a finer time discretization to avoid saturation, indicating the predominance of the error committed by the discretization in space over that in time in the problem (TC1B). This is due to the diffusivity of the implicit time scheme, which suffers in presence of discontinuities.

    The method is applied to the solution of Eq (2.2) with c=0, a constant diffusion coefficient ν=2 and the initial condition u0(x)=sin(πx), where x[0,1]. We want to approximate the exact solution

    uex(x,t)=sin(πx)e2π2t,x[0,1],t[0,1], (TC2)

    and impose Dirichlet boundary conditions at the end points accordingly: uex(0,t)=uex(1,t)=0, t(0,1].

    In Figure 6, we can observe that the theoretical orders of convergence for the implicit Euler scheme and FV-TPFA (1 and 2 respectively) are verified. The error saturates for very small sizes of the spatial mesh, due to the predominance of the error induced by the time discretization.

    Figure 6.  TC2 – Convergence test for the diffusion equation on a line. The normalized L1 error is computed at time t=1 corresponding to space and time meshes with different dimensions.

    The first branched domain we consider is represented in Figure 7, made of one bifurcation node and three edges of equal length |ek|=L=2, k=1,2,3, each further partitioned into segments of length h for spatial discretization. We set the advection speed c directed as the respective edges and constant on each of them

    ck=c|ek={10,k=2,5,k=1,3.
    Figure 7.  Simple one-dimensional domain with three edges and one bifurcation node Υ1. The set of edges is E={e1=IA,e2=BI,e3=IC}. The only source node is B and the set of end nodes is Ne={A,C}. There is one bifurcation node I.

    If we set the initial condition u0=0 on the whole domain Λ and the inflow boundary condition uB=1 on the vertex Υ2, we are able to compute the exact solution at each time step using the method of characteristics. On e2, we have the propagation of a rectangular wave with amplitude 1 and speed c2, which reaches the bifurcation point for t=t=|e2|c2=210=0.2

    uex(x,t)={1,if xc2t and 0tt, or t>t,0,if x>c2t, 0tt.

    Let u be the inflow condition on the edges e1 and e3 through the node I. We then obtain the propagation of a rectangular wave of amplitude u=c2c1+c3u2 and speed ck on each ek, k=1,3, starting at time t=t. Finally, the exact solution on ek, k=1,3, is the following:

    uex(x,t)={1,if xckt, t>t,0,otherwise.

    Thus, we expect a piecewise constant solution on the edges of the extended graph. However, in Figure 8, representing the numerical solution at three different times for h=102, we can observe that the jump discontinuities are not exactly captured. For instance, the solution at time t=0.2 should be u=2 on e2 and u=0 on e1 and e3, but a little dissipation is shown in proximity of the bifurcation, due to the diffusivity of the upwind scheme.

    Figure 8.  TC3 – Numerical solution at the initial time t=0, at t=0.2 when the edge e2 is full, and at the final time t=1. Black lines represent the domain Λ.

    Since, on each branch, the numerical method is a finite volume scheme with upwind flux, we expect to observe similar convergence properties as in the analogous test case (TC1B) presented in Section 5.1 with a discontinuous solution. Indeed, Figure 9 shows the convergence with order 12 both in time and space of the normalized L1 error.

    Figure 9.  TC3 – Convergence test for the transport equation on a graph with one bifurcation node. The normalized L1 error is computed at time t=1 corresponding to space and time meshes with different dimensions.

    The second test case on a graph is the solution of diffusion Eq (3.10) on the domain represented in Figure 10 with four edges and one bifurcation node.

    Figure 10.  Simple one-dimensional domain with four edges and one bifurcation node Υ1. The set of edges is E={e1=IA,e2=BI,e3=IC,e4=ID}. All the boundary nodes A, B, C, D are source nodes, and the only bifurcation node is I.

    We set diffusion coefficient ν=4 and impose homogeneous Dirichlet boundary conditions on the set of boundary nodes Nb={Υ0,Υ2,Υ3,Υ4} and the initial condition u0(s)=cos(π2s). This problem was constructed by starting from the following exact solution:

    uex(s,t)=cos(π2s)eπ2t.

    Due to the symmetry of the domain and of the solution, we expect the numerical solver to behave as it does on two separate straight lines and to have similar performances as in the test case (TC2) presented in Section 5.2. Indeed, in Figure 11, we can see that the orders of convergence 1 and 2 for time and space, respectively, are shown also in presence of a bifurcation. Finally, the numerical solution at three different times is represented in Figure 12 and shows a good approximation of the initial datum and of the expected decrease in the peak of the solution.

    Figure 11.  TC4 – Convergence test for the diffusion equation on a graph with one bifurcation node. The L1 error is computed at time t=1 corresponding to space and time meshes with different dimensions.
    Figure 12.  TC4 – Numerical solution at the initial time t=0, at t=0.1, and at the final time t=1. Black lines represent the domain Λ.

    Finally, we solve the drift-diffusion problem (2.2) on a ramified domain representing the skeleton of a typical electrical treeing, experimentally obtained by X-ray computed tomography [3] of an existing defect in an electrical cable, resulting in a graph composed of 12,544 edges. Such an extended defect can hardly be discretized in 3D due to its geometrical complexity. Thus, we need an ad hoc numerical solver to reduce it to PDEs on a 1D graph of problems defined on a 3D electrical treeing.

    This applied test case is derived by a geometrical reduction to one dimension of the model describing the movement of charge densities across the 3D electrical treeing, where the unknown u represents the integral mean of the charge concentration on transversal sections of the 3D domain. We consider a constant electron diffusion coefficient ν=0.5×106m2μs, while the value of the transport speed on each edge ek of the graph is given by ck=μ|Ek|, k=1,,Ne, where μ=105m2kV μs denotes the electronic mobility of the material comprising the domain and |Ek| is the magnitude of the electric field along ek, k=1,,Ne. We assume that the electric field is different on each edge of the graph, respecting the compatibility condition (4.1). The electric field on the edges connected to the inflow node is equal to 50 kV, and it is split equally among its outgoing neighbors, at each bifurcation. Figure 13 shows the value of the velocity corresponding to each edge of the graph.

    Figure 13.  TC5 – Drift velocities on the graph edges. The maximum value is on the edges connected to the inflow node, and it decreases towards the outflow node until the minimum value of order 1051.

    We assume that the set of source nodes consists of a single node, given by the root of the tree, where we impose the inflow electron concentration through the Dirichlet condition u=1001m2, and simulate over the time interval (0,500s], starting from the initial condition u0=01m2. The choice of such a long time interval for the simulation is for illustrative purposes only, due to the low transport velocity as we approach the periphery of the tree.

    For the time discretization, we consider a time step dt=0.1s, while the edges of the space mesh coincide with the edges of the graph, which are not further partitioned. The final result is shown in Figure 14, where the concentration imposed as the inflow condition on top of the domain is spread across the whole graph. Not all the edges are filled because, as displayed in Figure 13, the transport velocity becomes close to zero at the middle of the graph. However, we expect that, in the long term, the graph will all be filled, i.e., the concentration on all the edges will eventually be the same as those prescribed by the inflow condition.

    Figure 14.  TC5 – Numerical solution at the first time step t=0.1s and at the final time t=500s on the domain of an electrical treeing.

    We were able to simulate this problem in approximately 2 hours in parallel on six cores on a laptop with 16 GiB RAM and 11th Gen Intel(R) Core(TM) i7-11. In this test case, we were able to deal with a very complex and extended graph, solving the problem on a relatively coarse time grid and without further partitioning the edges of the graph, obtaining a positive solution at each time step, in reasonable time.

    We have defined two numerical schemes based on finite volumes for the approximate solution of a linear transport equation and a drift-diffusion problem on one-dimensional graph domains, with implicit time discretization. We have proven the existence and uniqueness of the solution to both discrete problems and that, starting from a non-negative initial condition and source term, the solutions at each time step remain non-negative. Moreover, we have shown that the numerical fluxes are consistent at the graph nodes. The positivity of the approximate solution at each time step is necessary for the problems of our interest, concerning plasma, where the presented methods are coupled with chemical solvers, which cannot deal with negative concentrations. Moreover, starting from general problems with few hypotheses and exploiting the properties of graph domains, we have shown that the numerical fluxes are consistent at graph nodes, provided that a compatibility condition is satisfied by the transport velocities. This condition is physically meaningful for different applications and thus is not restrictive.

    The two methods are extensions of a finite volume upwind scheme for the transport equation and TPFA for the diffusion part. We have observed that the convergence results of the proposed methods on a one-dimensional domain consisting in a straight line coincide with the theoretical results for these two schemes, and that this property is also extended to graph domains with bifurcation nodes.

    Finally, we have applied the two numerical methods to the solution of a transport equation and a drift-diffusion problem on the geometry of an electrical tree, and obtained a solution in a very short computational time, which was our starting motivation. The methods present good performance on real complex geometries with a moderate computational cost, despite the high number of branches. Moreover, we have made a parallel implementation of each solver, allowing us to exploit resources of high-performance computers. Unfortunately, both the upwind and implicit Euler methods are diffusive, but our choice fell on these schemes because of the necessity of a moderate computational cost. In fact, the possibility of adopting large time spacings Δt, due to the stability of the time numerical scheme, helps in reducing the computational time.

    We aim to incorporate the presented problems, describing the evolution of charge densities in non-thermal plasma, into a more comprehensive framework modelling the whole phenomenon of partial discharges and the evolution of electrical treeing. The goal is to overcome the employment of the semi empirical schemes used so far in the literature on the topic, keeping a controlled computational cost [50,57]. However, a possible extension could be the application of higher-order methods, such as Discontinuous Galerkin in space and Crank–Nicholson, which is still unconditionally stable, in time.

    Conceptualization, A. S. and A. V.; methodology, B. C., A. S. and A. V.; formal analysis, B. C.; writing-original draft preparation, B. C.; writing-review and editing, A. S. and A. V.. All authors have read and agreed to the published version of the manuscript.

    The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

    This work has been financed by the Research Found for the Italian Electrical System under the Contract Agreement between RSE (Ricerca Sul Sistema Energetico) and the Ministry of Economic Development. The mesh of the treeing structure is courtesy of Prof. R. Schurch, Universidad Técnica Federico Santa María Valparaíso, Chile.

    The authors declare there is no conflict of interest.



    [1] S. Bahadoorsingh, S. M. Rowland, The role of power quality in electrical treeing of epoxy resin, in 2007 Annual Report-Conference on Electrical Insulation and Dielectric Phenomena, IEEE, Vancouver, BC, Canada, (2007), 221–224. https://doi.org/10.1109/CEIDP.2007.4451500
    [2] G. Buccella, A. Villa, D. Ceresoli, R. Schurch, L. Barbieri, R. Malgesini, et al., A computational modelling of carbon layer formation on treeing branches, Modell. Simul. Mater. Sci. Eng., 31 (2023), 035001. https://doi.org/10.1088/1361-651X/acac44 doi: 10.1088/1361-651X/acac44
    [3] R. Schurch, S. M. Rowland, R. S. Bradley, P. J. Withers, Imaging and analysis techniques for electrical trees using X-ray computed tomography, IEEE Trans. Dielectr. Electr. Insul., 21 (2014), 53–63. https://doi.org/10.1109/TDEI.2013.003911 doi: 10.1109/TDEI.2013.003911
    [4] R. Schurch, J. Ardila-Rey, J. Montana, A. Angulo, S. M. Rowland, I. Iddrissu, et al., 3D characterization of electrical tree structures, IEEE Trans. Dielectr. Electr. Insul., 26 (2019), 220–228. https://doi.org/10.1109/TDEI.2018.007486 doi: 10.1109/TDEI.2018.007486
    [5] W. J. K. Raymond, H. A. Illias, H. Mokhlis, Partial discharge classifications: Review of recent progress, Measurement, 68 (2015), 164–181. https://doi.org/10.1016/j.measurement.2015.02.032 doi: 10.1016/j.measurement.2015.02.032
    [6] S. K. Pankaj, K. M. Keener, Cold plasma: Background, applications and current trends, Curr. Opin. Food Sci., 16 (2017), 49–52. https://doi.org/10.1016/j.cofs.2017.07.008 doi: 10.1016/j.cofs.2017.07.008
    [7] Y. Nyanteh, L. Graber, C. Edrington, S. Srivastava, D. Cartes, Overview of simulation models for partial discharge and electrical treeing to determine feasibility for estimation of remaining life of machine insulation systems, in 2011 Electrical Insulation Conference (EIC), IEEE, (2011), 327–332. https://doi.org/10.1109/EIC.2011.5996172
    [8] J. Jow, W. K. Lee, G. S. Cieloszyk, Stochastic simulation of water treeing in heterogeneous media using a field enhancement equation, in Proceedings of Conference on Electrical Insulation and Dielectric Phenomena-CEIDP'96, IEEE, 2 (1996), 758–761. https://doi.org/10.1109/CEIDP.1996.564619
    [9] G. Callender, P. L. Lewin, Plasma dynamic simulations of partial discharges within electrical tree structures, in 2019 IEEE Electrical Insulation Conference (EIC), IEEE, (2019), 340–343. https://doi.org/10.1109/EIC43217.2019.9046581
    [10] A. Villa, R. Schurch, L. Barbieri, G. Buccella, R. Malgesini, D. Palladini, Towards the plasma-polymer simulation in treeing branches, in 2022 IEEE 4th International Conference on Dielectrics (ICD), IEEE, (2022), 171–174. https://doi.org/10.1109/ICD53806.2022.9863222
    [11] A. Villa, L. Barbieri, M. Gondola, A. R. Leon-Garzon, R. Malgesini, A PDE-based partial discharge simulator, J. Comput. Phys., 345 (2017), 687–705. https://doi.org/10.1016/j.jcp.2017.05.045 doi: 10.1016/j.jcp.2017.05.045
    [12] A. Villa, L. Barbieri, M. Gondola, A. R. Leon-Garzon, R. Malgesini, An implicit three-dimensional fractional step method for the simulation of the corona phenomenon, Appl. Math. Comput., 311 (2017), 85–99. https://doi.org/10.1016/j.amc.2017.04.037 doi: 10.1016/j.amc.2017.04.037
    [13] A. Villa, L. Barbieri, M. Gondola, A. R. Leon-Garzon, R. Malgesini, Stability of the discretization of the electron avalanche phenomenon, J. Comput. Phys., 296 (2015), 369–381. https://doi.org/10.1016/j.jcp.2015.05.013 doi: 10.1016/j.jcp.2015.05.013
    [14] A. Villa, L. Barbieri, G. Marco, R. Malgesini, A. R. Leon-Garzon, Simulation of the AC corona phenomenon with experimental validation, J. Phys. D: Appl. Phys., 50 (2017), 435201. https://doi.org/10.1088/1361-6463/aa84f0 doi: 10.1088/1361-6463/aa84f0
    [15] E. Estradam, The Structure of Complex Networks: Theory and Applications, American Chemical Society, 2012.
    [16] M. Newman, Networks, Oxford University Press, 2018.
    [17] A. Bressan, S. Čanić, M. Garavello, M. Herty, B. Piccoli, Flows on networks: Recent results and perspectives, EMS Surv. Math. Sci., 1 (2014), 47–111. https://doi.org/10.4171/EMSS/2 doi: 10.4171/EMSS/2
    [18] B. Piccoli, M. Garavello, Traffic flow on networks, Am. Inst. Math. Sci., 2006.
    [19] S. Göttlich, M. Herty, A. Klar, Network models for supply chains, Commun. Math. Sci., 3 (2005), 545–559.
    [20] M. Herty, J. Mohring, V. Sachers, A new model for gas flow in pipe networks, Math. Methods Appl. Sci., 33 (2010), 845–855. https://doi.org/10.1002/mma.1197 doi: 10.1002/mma.1197
    [21] M. Herty, A.l Klar, B. Piccoli, Existence of solutions for supply chain models based on partial differential equations, SIAM J. Math. Anal., 39 (2007), 160–173. https://doi.org/10.1137/060659478 doi: 10.1137/060659478
    [22] J. Hild, G. Leugering, Real-time control of urban drainage systems, Math. Optim. Water Networks, (2012), 129–150.
    [23] W. A. M. Wybo, D. Boccalini, B. Torben-Nielsen, M. O. Gewaltig, A sparse reformulation of the green's function formalism allows efficient simulations of morphological neuron models, Neural Comput., 27 (2015), 2587–2622. https://doi.org/10.1162/NECO_a_00788 doi: 10.1162/NECO_a_00788
    [24] G. Bretti, R. Natalini, B. Piccoli, Numerical approximations of a traffic flow model on networks, Networks Heterogen. Media, 1 (2005), 57–84. https://doi.org/10.3934/nhm.2006.1.57 doi: 10.3934/nhm.2006.1.57
    [25] P. Goatin, E. Rossi, Comparative study of macroscopic traffic flow models at road junctions, Networks Heterogen. Media, 15 (2020), 261–279. https://doi.org/hal-02474650
    [26] V. Gyrya, A. Zlotnik, An explicit staggered-grid method for numerical simulation of large-scale natural gas pipeline networks, Appl. Math. Modell., 65 (2019), 34–51. https://doi.org/10.1016/j.apm.2018.07.051 doi: 10.1016/j.apm.2018.07.051
    [27] G. Leugering, Domain decomposition of an optimal control problem for semi-linear elliptic equations on metric graphs with application to gas networks, Appl. Math., 8 (2017), 1074. https://doi.org/10.4236/am.2017.88082 doi: 10.4236/am.2017.88082
    [28] S. K. Godunov, I. Bohachevsky, Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, Math. Sbornik, 47 (1959), 271–306.
    [29] D. Aregba-Driollet, V. Milišić, Kinetic approximation of a boundary value problem for conservation laws, Numer. Math., 97 (2004), 595–633. https://doi.org/10.1007/s00211-003-0514-5 doi: 10.1007/s00211-003-0514-5
    [30] J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2004.
    [31] J. E. Lagnese, G. Leugering, Domain in Decomposition Methods in Optimal Control of Partial Differential Equations, Springer Science & Business Media, 2004.
    [32] G. Berkolaiko, P. Kuchment, Introduction to Quantum Graphs, American Mathematical Society, 2013.
    [33] M. Benzi, G. H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numer., 14 (2005), 1–137. https://doi.org/10.1017/S0962492904000212 doi: 10.1017/S0962492904000212
    [34] M. Benzi, F. Durastante, F. Zigliotto, Modelling advection on distance-weighted directed networks, preprint, arXiv: 2410.11352.
    [35] J. C. Xu, L. Zikatanov, A monotone finite element scheme for convection-diffusion equations, Math. Comput., 68 (1999), 1429–1446.
    [36] K. Lipnikov, D. Svyatskiy, Y. Vassilevski, A monotone finite volume method for advection–diffusion equations on unstructured polygonal meshes, J. Comput. Phys., 229 (2010), 4017–4032. https://doi.org/10.1016/j.jcp.2010.01.035 doi: 10.1016/j.jcp.2010.01.035
    [37] W. Hundsdorfer, B. Koren, J. G. Verwer, A positive finite-difference advection scheme, J. Comput. Phys., 117 (1995), 35–46. https://doi.org/10.1006/jcph.1995.1042 doi: 10.1006/jcph.1995.1042
    [38] B. Crippa, A. Scotti, A. Villa, A mixed-dimensional model for the electrostatic problem on coupled domains, J. Comput. Phys., (2025), 114015. https://doi.org/10.1016/j.jcp.2025.114015
    [39] W. F. Fang, K. Ito, Global solutions of the time-dependent drift-diffusion semiconductor equations, J. Differ. Equ., 123 (1995), 523–566.
    [40] A. Jüngel, Drift-diffusion equations, Transp. Equ. Semicond., (2009), 1–29. https://doi.org/10.1007/978-3-540-89526-8_5
    [41] P. A. Markowich, C. A. Ringhofer, C. Schmeiser, Semiconductor Equations, Springer Science & Business Media, 2012.
    [42] D. L. Scharfetter, H. K. Gummel, Large-signal analysis of a silicon read diode oscillator, IEEE Trans. Electron Devices, 16 (1969), 64–77. https://doi.org/10.1109/T-ED.1969.16566 doi: 10.1109/T-ED.1969.16566
    [43] I. Berre, F. Doster, E. Keilegavlen, Flow in fractured porous media: A review of conceptual models and discretization approaches, Transp. Porous Media, 130 (2019), 215–236.
    [44] M. Karimi-Fard, L. J. Durlofsky, K. Aziz, An efficient discrete-fracture model applicable for general-purpose reservoir simulators, SPE J., 9 (2004), 227–236. https://doi.org/10.2118/88812-PA doi: 10.2118/88812-PA
    [45] T. Koch, Projection-based resolved interface 1D-3D mixed-dimension method for embedded tubular network systems, Comput. Math. Appl., 109 (2022), 15–29. https://doi.org/10.1016/j.camwa.2022.01.021 doi: 10.1016/j.camwa.2022.01.021
    [46] L. Cattaneo, P. Zunino, A computational model of drug delivery through microcirculation to compare different tumor treatments, Int. J. Numer. Methods Biomed. Eng., 30 (2014), 1347–1371. https://doi.org/10.1002/cnm.2661 doi: 10.1002/cnm.2661
    [47] J. Badwaik, A. M. Ruf, Convergence rates of monotone schemes for conservation laws with discontinuous flux, SIAM J. Numer. Anal., 58 (2020), 607–629. https://doi.org/10.1137/19M1283276 doi: 10.1137/19M1283276
    [48] J. von Below, Classical solvability of linear parabolic equations on networks, J. Differ. Equ., 72 (1988), 316–337. https://doi.org/10.1016/0022-0396(88)90158-1 doi: 10.1016/0022-0396(88)90158-1
    [49] J. Banasiak, P. Namayanja, Asymptotic behaviour of flows on reducible networks, Networks Heterogen. Media, 9 (2014). https://doi.org/10.3934/nhm.2014.9.197
    [50] A. Villa, R. Schurch, L. Barbieri, R. Malgesini, G. Buccella, An uncoupled implementation of the local mean energy plasma model, J. Comput. Phys., 447 (2021), 110674. https://doi.org/10.1016/j.jcp.2021.110674 doi: 10.1016/j.jcp.2021.110674
    [51] T. H. Sandve, I. Berre, J. M. Nordbotten, An efficient multi-point flux approximation method for Discrete Fracture–Matrix simulations, J. Comput. Phys., 231 (2012), 3784–3800. https://doi.org/10.1016/j.jcp.2012.01.023 doi: 10.1016/j.jcp.2012.01.023
    [52] R. Eymard, T. Gallouët, R. Herbin, Finite volume methods, Handb. Numer. Anal., 7 (2000), 713–1018. https://doi.org/10.1016/S1570-8659(00)07005-8
    [53] R. J. Plemmons, M-matrix characterizations. I—nonsingular M-matrices, Linear Algebra Appl., 18 (1977), 175–188. https://doi.org/10.1016/0024-3795(77)90073-8 doi: 10.1016/0024-3795(77)90073-8
    [54] R. A. Horn, C. R. Johnson, Matrix Analysis, Cambridge University Press, 2012.
    [55] F. Z. Zhang, The Schur Complement and its Applications, Springer Science & Business Media, 2006.
    [56] R. J. LeVeque, Numerical Methods for Conservation Laws, Springer Birkhäuser Basel, 1992.
    [57] A. Villa, R. Schurch, G. Buccella, L. Barbieri, C. Laurano, R. Malgesini, et al., Simulation of surface-plasma interaction with high surface conductivity, J. Comput. Phys., 456 (2022), 111029.
  • Reader Comments
  • © 2025 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(169) PDF downloads(17) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog