Mixed Mesh Finite Volume Method for 1D Hyperbolic Systems with Application to Plug-Flow Heat Exchangers

Abstract

We present a finite volume method formulated on a mixed Eulerian-Lagrangian mesh for highly advective 1D hyperbolic systems altogether with its application to plug-flow heat exchanger modeling/simulation. Advection of sharp moving fronts is an important problem in fluid dynamics, and even a simple transport equation cannot be solved precisely by having a finite number of nodes/elements/volumes. Finite volume methods are known to introduce numerical diffusion, and there exist a wide variety of schemes to minimize its occurrence; the most recent being adaptive grid methods such as moving mesh methods or adaptive mesh refinement methods. We present a solution method for a class of hyperbolic systems with one nonzero time-dependent characteristic velocity. This property allows us to rigorously define a finite volume method on a grid that is continuously moving by the characteristic velocity (Lagrangian grid) along a static Eulerian grid. The advective flux of the flowing field is, by this approach, removed from cell-to-cell interactions, and the ability to advect sharp fronts is therefore enhanced. The price to pay is a fixed velocity-dependent time sampling and a time delay in the solution. For these reasons, the method is best suited for systems with a dominating advection component. We illustrate the method’s properties on an illustrative advection-decay equation example and a 1D plug flow heat exchanger. Such heat exchanger model can then serve as a convection-accurate dynamic model in estimation and control algorithms for which it was developed.

1. Introduction

Numerical computation is a fundamental tool for the study of complex phenomena described by sets of partial differential equations (PDE). Being able to simulate these equations gives the opportunity to answer not only research questions but to accelerate technological and knowledge advances in general. To name a few, numerical solutions of PDEs are important in fluid dynamics, magnetohydrodynamics, particle physics, automatic control, optimization, etc. In this paper, we deal with a class of problems arising in computational fluid dynamics, in particular, the advection-(diffusion)-reaction (ADR) equation. The ADR equation represents a variety of physical, chemical, or biological processes, where the examples include water pollution analysis, chemical reactors models, heat exchanger models, population growth, and others.

There exists countless amount of solution methods for such a system. Coarsely, one can choose from finite difference, finite element, or finite volume methods (FVM). Picking the FVM, which differs from the latter by engaging the divergence theorem, the main ingredients of the solutions are a spatial reconstruction scheme, numerical flux selection, grid (adaptation) selection, and a time integration method. Analysis of numerical methods for hyperbolic systems is given in [1].

A breakthrough paper by Godunov [2] presents a first-order upwind method and discusses the monotonicity of the solution schemes. Regarding spatial accuracy, his method is of the first order. A second-order scheme is referred to as MUSCL according to its originating paper [3]. This was also the first second-order total variation diminishing (TVD) [4] scheme. Such a property is there obtained by the use of a proper flux limiter function, e.g., minmod or superbee limiters [5]. There are, in fact infinitely many TVD flux limiter functions; they only need to fit the TVD region [6]. Arbitrary order, piecewise polynomial, essentially non-oscillatory scheme (ENO) for spatial reconstruction was presented in [7] and further developed in [8,9]. The most used nowadays, however, is the weighted essentially non-oscillatory scheme (WENO) first presented in [10] and further developed and analyzed in [11].

Numerical flux is driven by the Riemann problem solution, which can be either exact, e.g., [2], or approximate as in [12,13]. The latest development in the field focuses on adaptive mesh techniques. Moving mesh methods continuously reposition a finite number of grid points such that some measure is equalized over the grid; a review can be found in [14]. Another technique, the adaptive mesh refinement method, is, however, the most popular. Between the recent well-cited papers are [15,16], for example. Time integration is either multi-step, Runge–Kutta typically, or one-step, where the discontinuous Galerkin methods gain popularity lately.

This article presents a mixed-mesh method for a particular class of hyperbolic systems—systems with two characteristics, where one is non-zero and time-dependent. Traditional FVM may add significant artificial numerical diffusion at the low number of states/nodes; the presented mixed-mesh method eliminates this drawback. Sharp edges are convected undistorted even in low-dimensional situations, which is beneficial for control-oriented models for example. The model is given in the state-space representation for the convenience of the control community. Full implementation in Matlab is published as an open-source [17]. An exact discretization of the continuous model is prepared for a follow-up article.

The paper is organized as follows. First, the problem formulation is stated and a brief resume of FVM is given. The mixed-mesh finite volume method (MMFVM) is defined in the next section An application to a simple advection equation, together with a heat exchanger application, is presented at the end of the article.

1.1. Problem Formulation

Consider a general continuity equation in one spatial dimension

$${\partial }_t\mathit{q}+{\partial }_x\mathit{f}=\mathit{s},$$

(1)

where the flux $\mathit{f}\left(x,t\right)$ is of the form

$$\mathit{f}=\mathit{V}\mathit{q}+\mathit{D}\nabla \mathit{q}$$

and the variables of interest (states) $\mathit{q}\left(x,t\right):\mathbb{R}\times \mathbb{R}\to {\mathbb{R}}^n$ are defined on a spatial domain $x\in \Omega$ and time domain $t\in \mathcal{J}$, $\mathit{V}\left(t\right)$ is a matrix of time-dependent transport velocities, $\mathit{D}$ is a diffusion coefficient matrix, $\mathit{s}\left(\mathit{q},x,t\right)\in {\mathbb{R}}^n$ is a source/sink term and ${\partial }_t=\partial /\partial t,{\partial }_x=\partial /\partial x$ are shorthand notations for partial differential operators with respect to time and space (Euler notation), respectively.

Under the assumption of an incompressible fluid, spatially invariant transport velocity and constant diffusion coefficients, we get a common advection-diffusion-reaction equation

$${\partial }_t\mathit{q}+\mathit{V}{\partial }_x\mathit{q}+\mathit{D}{\partial }_x^2\mathit{q}=\mathit{s}.$$

(2)

Focusing further on highly advective systems (with the Péclet number much larger than 10) the diffusion term is marginal and can be omitted [18]. Consider now the initial boundary value problem (IBVP) for the system (3), where ${\mathit{q}}_0$$\left(x\right)$ is an initial spatial profile at time $t=t_0$ and boundary conditions at $x=0$ are $\mathit{b}\left(t\right)$ (transport velocities are considered non-negative in the positive direction of x),

$$\begin{array}{ccc}\hfill {\partial }_t\mathit{q}+\mathit{V}{\partial }_x\mathit{q}& =& \mathit{s}\hfill \end{array}$$

(3a)

$$\begin{array}{ccc}\hfill \mathit{q}\left(x,t_0\right)& =& {\mathit{q}}_0\left(x\right)\hfill \end{array}$$

(3b)

$$\begin{array}{ccc}\hfill \mathit{q}\left(0,t\right)& =& \mathit{b}\left(t\right).\hfill \end{array}$$

(3c)

The goal is to obtain a transient solution for the variables $\mathit{q}\left(x,t\right)$ on a defined time interval $\mathcal{J}\equiv 〈t_0,T〉$. In particular, to obtain a numerical solution that preserves the wave behavior of the system without introducing numerical diffusion.

1.2. Finite Volume Method

Finite volume method (FVM) is one of the numerical solution techniques used for solving fluid dynamics problems. The core of the method lies in the usage of the divergence theorem on cells of finite volume and then on an approximation of the cell boundary fluxes from cell averages (a cell-centered method is assumed). More specifically, define a grid of N nodes, equidistant for simplicity, with a node to node distance $\Delta x=1/N$. For any $\mathit{h}:\mathbb{R}\times \mathbb{R}\to {\mathbb{R}}^n$ define a notation of an average value over one grid length,

$${\overline{\left(\mathit{h}\right)}}_x\left(t\right)=\frac{1}{\Delta x}{\int }_{c\left(x\right)}\mathit{h}\left(s,t\right)\mathrm{d}s,$$

(4)

where the interval $c\left(x\right)=〈x-\Delta x/2,x+\Delta x/2〉$. Now, applying the cell average to (3) and using the divergence theorem results in a system of ordinary differential equations for all cell-average variables

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathit{q}}_i+\frac{\mathit{V}}{\Delta x}\left({\mathit{q}}_{i-1/2}-{\mathit{q}}_{i+1/2}\right)={\mathit{s}}_i,$$

(5)

where ${\mathit{q}}_i={\overline{\left(\mathit{q}\right)}}_{x_i}$ is a cell average about position $x_i=\left(i-1/2\right)\Delta x,i=1,\cdots ,N$, ${\mathit{q}}_{i\pm 1/2}=\mathit{q}\left(x_i\pm \Delta x/2,t\right)$ denote left and right face values, and ${\mathit{s}}_i={\overline{\left(\mathit{s}\right)}}_{x_i}$ is a cell source term average. The face value approximation depends on a numerical scheme used. In the first-order upwind scheme, for example, are the face values simply ${\mathit{q}}_{i-1/2}={\mathit{q}}_{i-1},{\mathit{q}}_{i+1/2}={\mathit{q}}_i$. The first-order upwind approximation is monotone (no non-physical oscillations) but introduces a considerable numerical diffusion. Higher-order schemes bring higher accuracy, where the solution is smooth, but by the order barrier theorem, see [2] cannot be monotone. The remedy is to approximate the face values by a low order scheme, where the solution has high gradients, and by a high order scheme, where the solution is smooth. The switch between schemes is performed according to a slope/flux limiter function, and the method coined the term MUSCL [3]. Numerous limiters were introduced over the past five decades, starting with Van Leer [19]. An important class of slope limiters is the one forbidding an appearance of any new local extrema-the total variation diminishing class [4]. In the MUSCL scheme, the face values are approximated, component-wise, as

$$\begin{array}{ccc}\hfill q_{i-1/2}& =& q_{i-1}+\frac{\varphi \left(r_{i-1}\right)}{2}\left(q_{i-1}-q_{i-2}\right)\hfill \\ \hfill q_{i+1/2}& =& q_i+\frac{\varphi \left(r_i\right)}{2}\left(q_i-q_{i-1}\right),\hfill \end{array}$$

where the slope limiter function is, for example, the one by Van Leer [19]

$$\varphi \left(r\right)=\frac{r+\left|r\right|}{1+\left|r\right|}$$

and the smoothness monitor is

$$r_i=\frac{q_i-q_{i-1}}{q_{i+1}-q_i}.$$

(6)

See [20,21] for more schemes and slope limiters.

Time evolution of (5) is obtained using standard integration techniques such as RK4; special care, however, has to be taken concerning the sampling period. To contain the physical domain of dependence in the numerical domain of dependence (or in other words avoid interaction of adjacent Riemann problems), the sampling period must satisfy the Courant–Friedrichs–Lewy (CFL) condition [22]

$$max\left(\left|\mathrm{eig}\left(\mathit{V}\right)\right|\right)\Delta t\le \Delta x.$$

(7)

2. Mixed Mesh Finite Volume Method

The challenge of the numerical approximation of hyperbolic transport equations is to obtain high accuracy of the solution in both discontinuous and smooth regions [23]. Highly advective systems also present a challenge for the finite volume method. The cells are considered coherent (ideally mixed), and as a consequence, the flux at one boundary immediately influences fluxes at all the other boundaries. The numerical propagation speed is therefore always infinite for FVM. Pure discontinuities, as well, always dissolve into smooth step changes. The finer the grid, the less the phenomenon is observed. An exact solution, advecting discontinuities sharp, is obtained for an infinite number of cells [24].

To overcome this fundamental limitation, we propose a modified approach.

Remark 1.

Where coherent cells of finite volumes are used, an exact finite transport speed is only preserved by moving the grid by this velocity.

A Mixed-mesh Finite Volume Method (MMFVM) based on the observation of Remark 1 is presented here. The underlying principle is simple: advect (move) cells of the flowing field and simultaneously solve FVM for both the moving and static cells until the moving grid passes the static grid by one grid length; sample the state and reinitialize all moving cells to their starting positions and assign appropriate cell averages. The idea was first introduced in preceding conference papers of the authors [25,26] and will be generalized here. We believe a similar approach was used in [27,28], but a description of the spatially discretized system and its analysis are missing therein.

Definition 1 (Monocharacteristic hyperbolic system).

A system of the form (1), where the flux Jacobian $\mathit{J}\equiv \partial \mathit{f}/\partial \mathit{q}$ can be decomposed to

$$\mathit{J}=v\mathit{R}\mathit{D}{\mathit{R}}^{-1},\mathit{D}=\mathrm{diag}\left(d_1,\dots ,d_n\right),d_i\in \left\{0,1\right\}\forall i$$

(8)

with $v\left(t\right)\in \mathbb{R}$ being (for this paper) a space-independent velocity and $\mathit{R}$ being a square decomposition matrix, is a monocharacteristic hyperbolic system, i.e., the system has a monocharacteristic hyperbolic property.

Remark 2.

As a hyperbolic system has a diagonalizable flux Jacobian with real eigenvalues a system is monocharacteristic hyperbolic when the Jacobian has just one nonzero eigenvalue, i.e., when a subset of the vector variable$\mathit{q}$advects with a transport (characteristic) velocity v, and the rest of $\mathit{q}$ is stationary.

Definition 2 (Normal monocharacteristic form).

A normal monocharacteristic form is

$${\partial }_t\mathit{q}+v\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\partial }_x\mathit{q}=\mathit{s}$$

(9)

or equivalently

$$\begin{array}{ccc}\hfill {\partial }_t\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}\\ {\mathit{q}}_{\mathrm{s}}\end{array}\right]+v\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\partial }_x\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}\\ {\mathit{q}}_{\mathrm{s}}\end{array}\right]& =& \left[\begin{array}{c}{\mathit{s}}_{\mathrm{a}}\\ {\mathit{s}}_{\mathrm{s}}\end{array}\right],\hfill \end{array}$$

(10)

where $\mathit{q}={\left[\begin{array}{cc}{\mathit{q}}_{\mathrm{a}}^{\mathsf{T}}& {\mathit{q}}_{\mathrm{s}}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}$, $\mathit{s}={\left[\begin{array}{cc}{\mathit{s}}_{\mathrm{a}}^{\mathsf{T}}& {\mathit{s}}_{\mathrm{s}}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}$ and ${\mathit{q}}_{\mathrm{a}}\left(x,t\right):\mathbb{R}\times \mathbb{R}\to {\mathbb{R}}^{n_{\mathrm{a}}}$ is a vector of advected states, ${\mathit{q}}_{\mathrm{s}}\left(x,t\right):\mathbb{R}\times \mathbb{R}\to {\mathbb{R}}^{n_{\mathrm{s}}}$ is a vector of stationary states and${\mathit{s}}_{\mathrm{a}}\left(\mathit{q},\mathit{u},t\right)\in {\mathbb{R}}^{n_{\mathrm{a}}},{\mathit{s}}_{\mathrm{s}}\left(\mathit{q},\mathit{u},t\right)\in {\mathbb{R}}^{n_{\mathrm{s}}}$are respective spatially independent sink/source terms. $\mathit{I},\mathbf{0}$ are square unit and zero matrices of appropriate sizes and $\left(·\right)$ denotes rectangular zero matrix of coresponding size.

Lemma 1.

A system (1) with a monocharacteristic property (8) has a normal monocharacteristic form (9), (10).

Proof. 

We will prove the statement by transformation construction. The system (1) can be decomposed using the chain rule as

$${\partial }_t\mathit{q}+{\left.\frac{\partial \mathit{f}}{\partial \mathit{q}}\right|}_{x=\mathrm{const}.}{\partial }_x\mathit{q}=\mathit{s}-{\left.\frac{\partial }{\partial x}\mathit{f}\right|}_{\mathit{q}=\mathrm{const}.},$$

where $\mathit{J}\equiv {\left.\partial \mathit{f}/\partial \mathit{q}\right|}_{x=\mathrm{const}.}$ is the flux Jacobian and $\overline{\mathit{s}}\equiv \mathit{s}-{\left.\partial /\partial x\mathit{f}\right|}_{\mathit{q}=\mathrm{const}.}$ is an equivalent sink/source term. As the system is monocharacteristic hyperbolic, $\mathit{J}$ can be decomposed to form (8). Defining now a permutation matrix $\mathit{P}$ such that

$${\mathit{P}}^{\mathsf{T}}\mathit{D}\mathit{P}=\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right],$$

the Jacobian is

$$\mathit{J}=v\mathit{R}\mathit{P}\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\left(\mathit{R}\mathit{P}\right)}^{-1}.$$

So (1) can be rewritten as

$$\begin{array}{ccc}\hfill {\partial }_t\mathit{q}+v\mathit{R}\mathit{P}\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\left(\mathit{R}\mathit{P}\right)}^{-1}{\partial }_x\mathit{q}& =& \overline{\mathit{s}}\hfill \\ \hfill {\left(\mathit{R}\mathit{P}\right)}^{-1}{\partial }_t\mathit{q}+v\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\left(\mathit{R}\mathit{P}\right)}^{-1}{\partial }_x\mathit{q}& =& {\left(\mathit{R}\mathit{P}\right)}^{-1}\overline{\mathit{s}}.\hfill \end{array}$$

Now, define a transformation $T\left(x\right)={\left(\mathit{R}\mathit{P}\right)}^{-1}x$, then $\tilde{\mathit{q}}\equiv T\left(\mathit{q}\right)$, $\tilde{\mathit{s}}\equiv T\left(\overline{\mathit{s}}\right)$ and

$${\partial }_t\tilde{\mathit{q}}+v\left[\begin{array}{cc}\mathit{I}& ·\\ ·& \mathbf{0}\end{array}\right]{\partial }_x\tilde{\mathit{q}}=\tilde{\mathit{s}}$$

is in a normal monocharacteristic form (9). □

Method Formulation

Consider a problem (1) with a monocharacteristic property (8) without a loss of generality in the monocharacteristic form (10)

$$\begin{array}{ccc}\hfill {\partial }_t{\mathit{q}}_{\mathrm{a}}+v{\partial }_x{\mathit{q}}_{\mathrm{a}}& =& {\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}},{\mathit{q}}_{\mathrm{s}},\mathit{u}\right)\hfill \end{array}$$

(11a)

$$\begin{array}{ccc}\hfill {\partial }_t{\mathit{q}}_{\mathrm{s}}& =& {\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}},{\mathit{q}}_{\mathrm{s}},\mathit{u}\right)\hfill \end{array}$$

(11b)

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{a}}\left(0,t\right)& =& \mathit{b}\left(t\right)\hfill \end{array}$$

(11c)

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{a}}\left(x,0\right)& =& {\mathit{q}}_{\mathrm{a},0}\left(x\right)\hfill \end{array}$$

(11d)

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{s}}\left(x,0\right)& =& {\mathit{q}}_{\mathrm{s},0}\left(x\right),\hfill \end{array}$$

(11e)

where $\mathit{u}$$\left(t\right)$ are external inputs, initial spatial profiles for advected states and static states are ${\mathit{q}}_{\mathrm{a},0}\left(x\right)$, ${\mathit{q}}_{\mathrm{s},0}\left(x\right)$, respectively and $v\left(t\right)\ge 0$. Define an equidistant grid of N elements with a grid spacing $\Delta x$. The stationary grid ${\Xi }_{\mathrm{s}}$ resembles the one of FVM

$${\Xi }_{\mathrm{s}}=\left\{\left.x_{\mathrm{s}}^i\right|i\in \mathbb{Z},0\le i\le N,\Delta x=\frac{1}{N},x_{\mathrm{s}}^i=\left(i-\frac{1}{2}\right)\Delta x\right\},$$

(12)

but the advected grid differs importantly. The basic idea behind the MMFVM is to allow the grid of advected states to flow along the stationary states with the characteristic velocity v.

Definition 3 (Relative cell shift).

A relative stationary-to-advected grid shift p is a distance a tracer would travel with a velocity v in a time interval $〈t_0,t_0+\tau 〉$ taken relative to one grid length $\Delta x$

$$p\left(t\right)=\frac{1}{\Delta x}{\int }_{t_0}^{t_0+\tau }v\left(s\right)\mathrm{d}s,$$

(13)

where the integration start $t_0$ will be specified later on.

Definition 4 (Advected grid).

Grid positions of the advected states are

$${\Xi }_{\mathrm{a}}=\left\{\left.x_{\mathrm{a}}^i\right|i\in \mathbb{Z},0\le i\le N,x_{\mathrm{a}}^i=x_{\mathrm{s}}^i+p\Delta x\right\},$$

(14)

where p is given by (13).

Remark 3.

Note that there is always at least one upwind “ghost” cell for the advected states at position $x_{\mathrm{a}}^0=\left(p-1/2\right)\Delta x$, where a boundary condition is applied. Other ghost cells are added depending on the boundary conditions. See Figure 1 for illustration. Note that there is also a ghost cell at position $x_{\mathrm{s}}^0$ that can be used to apply boundary conditions for static states, but this cell in this paper does not play any role and is only added for future use and simplification of the notation.

Theorem 1 (MMFVM step).

Time evolution of a problem (11) discretized on grids ${\Xi }_{\mathrm{s}}$, ${\Xi }_{\mathrm{a}}$, is on a time interval $〈t_0,t_0+\tau 〉$ given by

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}^i\\ {\mathit{q}}_{\mathrm{s}}^i\end{array}\right]& =& \left(1-p\right)\left[\begin{array}{c}{\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\\ {\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\end{array}\right]\hfill \end{array}$$

$$\begin{array}{ccc}& & +p\left[\begin{array}{c}{\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^{i+1},\mathit{u}\right)\\ {\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^{i-1},{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\end{array}\right]\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}p& =& \frac{v}{\Delta x},\hfill \end{array}$$

(16)

where ${\mathit{q}}_{\mathrm{a}}^i$, ${\mathit{q}}_{\mathrm{s}}^i$, $i=0,\dots ,N$, denote centered cell average values of the advective and static grid, respectively, and ${\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^0,{\mathit{q}}_{\mathrm{s}}^0,\mathit{u}\right)=\mathbf{0}$, ${\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^N,{\mathit{q}}_{\mathrm{s}}^{N+1},\mathit{u}\right)=\mathbf{0}$, ${\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^0,{\mathit{q}}_{\mathrm{s}}^0,\mathit{u}\right)=\mathbf{0}$, ${\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^{-1},{\mathit{q}}_{\mathrm{s}}^0,\mathit{u}\right)=\mathbf{0}$. The length of the time interval τ is given by the relative cell shift to integrate from $p\left(t_0\right)=0$ to $p\left(t_0+\tau \right)=1$ or equivalently

$${\int }_{t_0}^{t_0+\tau }v\left(t\right)\mathrm{d}t=\Delta x.$$

(17)

The initial conditions are

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{a}}^i\left(t_0\right)& =& \overline{{\left({\mathit{q}}_{\mathrm{a}}\left(x,t_0\right)\right)}_{x_{\mathrm{s}}^i}}\hfill \end{array}$$

(18a)

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{s}}^i\left(t_0\right)& =& \overline{{\left({\mathit{q}}_{\mathrm{s}}\left(x,t_0\right)\right)}_{x_{\mathrm{s}}^i}}\hfill \end{array}$$

(18b)

$$\begin{array}{ccc}\hfill p\left(t_0\right)& =& 0,\hfill \end{array}$$

(18c)

for $i=1,\dots ,N$. The boundary conditions are applied to the upwind ghost cell${\mathit{q}}_{\mathrm{a}}^0\left(t_0\right)=\mathit{b}\left(t_0\right)$ at the beginning of the integration.

Proof. 

Direct. Spatial discretization of (11), by applying cell average operator (4) about positions $x_{\mathrm{s}}^i,x_{\mathrm{a}}^i$, gives

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}+{\overline{\left(v{\partial }_x{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& {\overline{\left({\mathit{s}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& {\overline{\left({\mathit{s}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}\hfill \end{array}$$

(20)

and to prove Theorem 1 we will deal with all the terms in (19) and (20) separately. Remember that $x_{\mathrm{a}}^i\left(t\right)=x_{\mathrm{s}}^i+p\left(t\right)\Delta x$.

The time differential term (19) reads as

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \frac{1}{\Delta x}{\int }_{x_{\mathrm{a}}^i\left(t\right)-\frac{\Delta x}{2}}^{x_{\mathrm{a}}^i\left(t\right)+\frac{\Delta x}{2}}{\partial }_t{\mathit{q}}_{\mathrm{a}}\left(x,t\right)\mathrm{d}x,\hfill \end{array}$$

where using the Leibnitz’s integral rule the order of integration and differentiation can be switched to obtain

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \frac{1}{\Delta x}\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{x_{\mathrm{a}}^i\left(t\right)-\frac{\Delta x}{2}}^{x_{\mathrm{a}}^i\left(t\right)+\frac{\Delta x}{2}}{\mathit{q}}_{\mathrm{a}}\left(x,t\right)\mathrm{d}x-\hfill \\ & & \frac{1}{\Delta x}{\mathit{q}}_{\mathrm{a}}\left(x_{\mathrm{a}}^i\left(t\right)-\frac{\Delta x}{2},t\right)\frac{\mathrm{d}}{\mathrm{d}t}\left(x_{\mathrm{a}}^i\left(t\right)-\frac{\Delta x}{2}\right)+\hfill \\ & & \frac{1}{\Delta x}{\mathit{q}}_{\mathrm{a}}\left(x_{\mathrm{a}}^i\left(t\right)+\frac{\Delta x}{2},t\right)\frac{\mathrm{d}}{\mathrm{d}t}\left(x_{\mathrm{a}}^i\left(t\right)+\frac{\Delta x}{2}\right).\hfill \end{array}$$

Now, first calculating the integral bounds derivative,

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left(x_{\mathrm{a}}^i\left(t\right)\pm \frac{\Delta x}{2}\right)& =& \frac{\mathrm{d}}{\mathrm{d}t}\left(x_{\mathrm{s}}^i+p\left(t\right)\Delta x\right)\hfill \\ & =& \Delta x\frac{\mathrm{d}}{\mathrm{d}t}\left(p\left(t\right)\right)\hfill \\ & =& \Delta x\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{t_0}^{\mathsf{T}}\frac{v\left(s\right)}{\Delta x}\mathrm{d}s\hfill \\ & =& v\left(t\right),\hfill \end{array}$$

we obtain the first term as

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \frac{1}{\Delta x}\left(\Delta x\frac{d}{dt}{\mathit{q}}_{\mathrm{a}}^i-v\left({\mathit{q}}_{\mathrm{a}}^{i+1/2}-{\mathit{q}}_{\mathrm{a}}^{i-1/2}\right)\right)\hfill \\ & =& \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{q}}_{\mathrm{a}}^i-\frac{v}{\Delta x}\left({\mathit{q}}_{\mathrm{a}}^{i+1/2}-{\mathit{q}}_{\mathrm{a}}^{i-1/2}\right),\hfill \end{array}$$

(21)

where ${\mathit{q}}_{\mathrm{a}}^i\left(t\right)=\overline{{\left({\mathit{q}}_{\mathrm{a}}\right)}_{x_{\mathrm{a}}^i}}\left(t\right)$ is a cell average value and ${\mathit{q}}_{\mathrm{a}}^{i\pm 1/2}\left(t\right)\equiv {\mathit{q}}_{\mathrm{a}}\left(x_{\mathrm{a}}^i\pm \Delta x/2,t\right)$ denote cell face values.

The transport term (19) reads as

$$\begin{array}{ccc}\hfill {\overline{\left(v{\partial }_x{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \frac{1}{\Delta x}{\int }_{x_{\mathrm{a}}^i-\frac{\Delta x}{2}}^{x_{\mathrm{a}}^i+\frac{\Delta x}{2}}v{\partial }_x{\mathit{q}}_{\mathrm{a}}\left(x,t\right)\mathrm{d}x,\hfill \end{array}$$

which using the Ostrogradsky’s divergence theorem simplifies to

$$\begin{array}{ccc}\hfill {\overline{\left(v{\partial }_x{\mathit{q}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \frac{v}{\Delta x}\left({\mathit{q}}_{\mathrm{a}}^{i+1/2}-{\mathit{q}}_{\mathrm{a}}^{i-1/2}\right).\hfill \end{array}$$

(22)

The time differential term (20) reads as

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& \frac{1}{\Delta x}{\int }_{x_{\mathrm{s}}^i-\frac{\Delta x}{2}}^{x_{\mathrm{s}}^i+\frac{\Delta x}{2}}{\partial }_t{\mathit{q}}_{\mathrm{s}}\left(x,t\right)\mathrm{d}x\hfill \end{array}$$

which by substitution ${\mathit{q}}_{\mathrm{s}}^i\left(t\right)={\overline{\left({\mathit{q}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}\left(t\right)$ immediately gives

$$\begin{array}{ccc}\hfill {\overline{\left({\partial }_t{\mathit{q}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{q}}_{\mathrm{s}}^i.\hfill \end{array}$$

(23)

The static state source term (20) reads as

$$\begin{array}{ccc}\hfill {\overline{\left({\mathit{s}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& \frac{1}{\Delta x}{\int }_{x_{\mathrm{s}}^i-\frac{\Delta x}{2}}^{x_{\mathrm{s}}^i+\frac{\Delta x}{2}}{\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}\left(x,t\right),{\mathit{q}}_{\mathrm{s}}\left(x,t\right),\mathit{u}\left(t\right)\right)\mathrm{d}x.\hfill \end{array}$$

Given ${\mathit{s}}_{\mathrm{s}}$ is spatially independent and inspecting Figure 1 we observe, that the integral may be separated as

$$\begin{array}{ccc}\hfill {\overline{\left({\mathit{s}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& \frac{1}{\Delta x}\left({\int }_{x_{\mathrm{s}}^i-\frac{\Delta x}{2}}^{x_{\mathrm{s}}^i-\frac{\Delta x}{2}+p\left(t\right)\Delta x}{\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}\left(x,t\right),{\mathit{q}}_{\mathrm{s}}\left(x,t\right),\mathit{u}\left(t\right)\right)\mathrm{d}x\right.\hfill \\ & & \left.+{\int }_{x_{\mathrm{s}}^i-\frac{\Delta x}{2}+p\left(t\right)\Delta x}^{x_{\mathrm{s}}^i+\frac{\Delta x}{2}}{\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}\left(x,t\right),{\mathit{q}}_{\mathrm{s}}\left(x,t\right),\mathit{u}\left(t\right)\right)\mathrm{d}x\right),\hfill \end{array}$$

which, under a piecewise constant spatial profile approximation of ${\mathit{q}}_{\mathrm{a}}$, ${\mathit{q}}_{\mathrm{s}}$ is

$$\begin{array}{ccc}\hfill {\overline{\left({\mathit{s}}_{\mathrm{s}}\right)}}_{x_{\mathrm{s}}^i}& =& \left(1-p\right){\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \\ & & +p{\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^{i-1},{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \end{array}$$

(24)

and finally following the same reasoning as above we obtain the advected state source term (19)

$$\begin{array}{ccc}\hfill {\overline{\left({\mathit{s}}_{\mathrm{a}}\right)}}_{x_{\mathrm{a}}^i}& =& \left(1-p\right){\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \\ & & +p{\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^{i+1},\mathit{u}\right).\hfill \end{array}$$

(25)

The source term formulas are only valid for $0\le p\left(t\right)\le 1$, which corresponds with (17). Note that higher order spatial reconstructions are possible following that appropriate calculation of the source terms (24), (25) above is performed.

Now, by substitution of (21) to (25) into (19), (20) we get

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{q}}_{\mathrm{a}}^i-\frac{v}{\Delta x}\left({\mathit{q}}_{\mathrm{a}}^{i+1/2}-{\mathit{q}}_{\mathrm{a}}^{i-1/2}\right)\\ \hfill +\frac{v}{\Delta x}\left({\mathit{q}}_{\mathrm{a}}^{i+1/2}-{\mathit{q}}_{\mathrm{a}}^{i-MM1/2}\right)& =& \left(1-p\right){\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \\ & & +p{\mathit{s}}_{\mathrm{a}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^{i+1},\mathit{u}\right)\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{q}}_{\mathrm{s}}^i& =& \left(1-p\right){\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^i,{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \\ & & +p{\mathit{s}}_{\mathrm{s}}\left({\mathit{q}}_{\mathrm{a}}^{i-1},{\mathit{q}}_{\mathrm{s}}^i,\mathit{u}\right)\hfill \end{array}$$

moreover, we see that the velocity-dependent terms cancel out, and by aggregation of states, we obtain (15). Differentiation of the relative cell shift Equation (13),

$$\frac{\mathrm{d}}{\mathrm{d}t}p\left(t\right)=\frac{v\left(t\right)}{\Delta x}$$

with $p\left(t_0\right)=0$ gives (16) and concludes the statement. □

Corollary 1 (MMFVM step for systems with a linear source term).

If a monocharacteristic hyperbolic system of Theorem 1 has a linear source term

$$\mathit{s}\left({\mathit{q}}_{\mathrm{a}},{\mathit{q}}_{\mathrm{s}},\mathit{u}\right)=\left[\begin{array}{cc}{\mathit{A}}_{\mathrm{aa}}& {\mathit{A}}_{\mathrm{as}}\\ {\mathit{A}}_{\mathrm{sa}}& {\mathit{A}}_{\mathrm{ss}}\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}\\ {\mathit{q}}_{\mathrm{s}}\end{array}\right]+\left[\begin{array}{c}{\mathit{B}}_{\mathrm{a}}\\ {\mathit{B}}_{\mathrm{s}}\end{array}\right]\mathit{u},$$

(26)

then (15) simplifies to

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}^i\\ {\mathit{q}}_{\mathrm{s}}^i\end{array}\right]& =& \left[\begin{array}{cc}{\mathit{A}}_{\mathrm{aa}}& {\mathit{A}}_{\mathrm{as}}\\ {\mathit{A}}_{\mathrm{sa}}& {\mathit{A}}_{\mathrm{ss}}\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_{\mathrm{a}}^i\\ {\mathit{q}}_{\mathrm{s}}^i\end{array}\right]+\hfill \\ & & +p\left[\begin{array}{cc}\mathbf{0}& {\mathit{A}}_{\mathrm{as}}\\ {\mathit{A}}_{\mathrm{sa}}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}\left({\mathit{q}}_{\mathrm{a}}^{i-1}-{\mathit{q}}_{\mathrm{a}}^i\right)\\ \left({\mathit{q}}_{\mathrm{s}}^{i+1}-{\mathit{q}}_{\mathrm{s}}^i\right)\end{array}\right]+\hfill \\ & & +\left[\begin{array}{c}{\mathit{B}}_{\mathrm{a}}\\ {\mathit{B}}_{\mathrm{s}}\end{array}\right]\mathit{u}\hfill \end{array}$$

(27)

and the whole system of equations can be formulated as

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathit{x}& =& \left({\mathit{A}}_{\mathrm{c}}+p{\mathit{A}}_{\mathrm{cp}}\right)\mathit{x}+{\mathit{B}}_{\mathrm{c}}\mathit{u},\hfill \end{array}$$

(28)

where the state vector $\mathit{x}\in {\mathbb{R}}^{\left(N+1\right)·n}$ aggregates the advected states and stationary states as

$$\mathit{x}={\left[\begin{array}{c}{\left({\mathit{q}}_{\mathrm{a}}^0\right)}^{\mathsf{T}},\dots ,{\left({\mathit{q}}_{\mathrm{a}}^N\right)}^{\mathsf{T}},{\left({\mathit{q}}_{\mathrm{s}}^0\right)}^{\mathsf{T}},\dots ,{\left({\mathit{q}}_{\mathrm{s}}^N\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}.$$

Proof. 

By construction. Substituting (26) into (15), after minor arrangements, gives (27). Defining now ${\tilde{\mathit{A}}}_{\mathrm{c},1}\equiv \left[\begin{array}{cc}{\mathit{A}}_{\mathrm{aa}}& {\mathit{A}}_{\mathrm{as}}\\ {\mathit{A}}_{\mathrm{sa}}& {\mathit{A}}_{\mathrm{ss}}\end{array}\right]$, ${\tilde{\mathit{A}}}_{\mathrm{c},2}\equiv \left[\begin{array}{cc}\mathbf{0}& \mathbf{0}\\ {\mathit{A}}_{\mathrm{sa}}& \mathbf{0}\end{array}\right]$, ${\tilde{\mathit{A}}}_{\mathrm{c},3}\equiv \left[\begin{array}{cc}\mathbf{0}& {\mathit{A}}_{\mathrm{as}}\\ \mathbf{0}& \mathbf{0}\end{array}\right]$, ${\tilde{\mathit{B}}}_{\mathrm{c}}\equiv \left[\begin{array}{c}{\mathit{B}}_{\mathrm{a}}\\ {\mathit{B}}_{\mathrm{s}}\end{array}\right]$, $M=N+1$ and aggregating the states as

$$\tilde{\mathit{x}}={\left[\begin{array}{c}{\left[\begin{array}{cc}{\left({\mathit{q}}_{\mathrm{a}}^0\right)}^{\mathsf{T}}& {\left({\mathit{q}}_{\mathrm{s}}^0\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}},\dots ,{\left[\begin{array}{cc}{\left({\mathit{q}}_{\mathrm{a}}^{N+1}\right)}^{\mathsf{T}}& {\left({\mathit{q}}_{\mathrm{s}}^{N+1}\right)}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}}\end{array}\right]}^{\mathsf{T}},$$

we can summarize (27) for $i=0,\dots ,M$ into a formula

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\tilde{\mathit{x}}& =& \left({\mathit{I}}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},1}\right)\tilde{\mathit{x}}\hfill \\ & +& p\left({\mathit{K}}_{-1}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},2}+{\mathit{K}}_{+1}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},3}\right)\tilde{\mathit{x}}\hfill \\ & +& \left({\mathbf{1}}^M\otimes {\mathit{I}}^n\otimes {\tilde{\mathit{B}}}_{\mathrm{c}}\right)\mathit{u},\hfill \end{array}$$

(29)

where ${\mathit{K}}_{-1}^M=\left({\mathit{I}}_{-1}^M-{\mathit{I}}^M\right)$, ${\mathit{K}}_{+1}^M=\left({\mathit{I}}_{+1}^M-{\mathit{I}}^M\right)$, ${\mathit{I}}^k$ denotes an identity matrix of size ${\mathbb{R}}^{k\times k}$, ${\mathit{I}}_{-1}^k$, ${\mathit{I}}_{+1}^k$ are lower and upper shift matrices of size ${\mathbb{R}}^{k\times k}$, respectively, ${\mathbf{1}}^k\in {\mathbb{R}}^k$ is a column vector of ones and ⊗ denotes the Kronecker product. Now, define a permutation matrix $\mathit{P}$ such that $\mathit{x}={\mathit{P}}^{\mathsf{T}}\tilde{\mathit{x}}$

$$\mathit{P}=\tilde{\mathit{P}}\otimes {\mathbf{1}}^2,\tilde{\mathit{P}}=\left\{p_{i,j}\right\}=〈\begin{array}{cc}1& i\mathrm{odd},j=\frac{i+1}{2}\\ 1& i\mathrm{even},j=\frac{i+2N}{2}\\ 0& \mathrm{otherwise}\end{array},$$

(30)

then (29) becomes (28) with

$$\begin{array}{ccc}\hfill \mathit{A}{}_{\mathrm{c}}& =& {\mathit{P}}^{\mathsf{T}}\left({\mathit{I}}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},1}\right)\mathit{P}\hfill \end{array}$$

(31a)

$$\begin{array}{ccc}\hfill \mathit{A}{}_{\mathrm{cp}}& =& {\mathit{P}}^{\mathsf{T}}\left(\left({\mathit{K}}_{-1}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},2}+{\mathit{K}}_{+1}^M\otimes {\tilde{\mathit{A}}}_{\mathrm{c},3}\right)\right)\mathit{P}\hfill \end{array}$$

(31b)

$$\begin{array}{ccc}\hfill {\mathit{B}}_{\mathrm{c}}& =& {\mathit{P}}^{\mathsf{T}}\left({\mathbf{1}}^M\otimes {\mathit{I}}^n\otimes {\tilde{\mathit{B}}}_{\mathrm{c}}\right).\hfill \end{array}$$

(31c)

□

Theorem 2 (Mixed Mesh Finite Volume Method, MMFVM).

A numerical solution to a monocharacteristic hyperbolic system is obtained by a repetitive solution of the MMFVM step (Theorem 1/ Corollary 1).

There is a delay τ in the solution

$${\tau }^i\left(t\right):\left\{\begin{array}{cc}{\int }_{t-{\tau }^i\left(t\right)}^{t}v\left(s\right)\mathrm{d}s=\frac{\Delta x}{2};\hfill & x=x_{\mathrm{s}}^i,i=1,\dots ,N\hfill \\ {\int }_{t-{\tau }^i\left(t\right)}^{t}v\left(s\right)\mathrm{d}s=\Delta x;\hfill & x=x_{\mathrm{s}}^i,i=N+1\hfill \end{array}\right.$$

(32)

caused by the fact that there is an additional ${\left(N+1\right)}^{\mathrm{st}}$ cell in the advected grid.

Proof. 

During the MMFVM step, the cell average values evolve by the governing PDE, as proved in the Theorem 1. Consider now a ghost cell ${\mathit{q}}_{\mathrm{a}}^0$ initialized by an appropriate boundary condition. It takes $N+1$ MMFVM steps to move this cell along all the stationary cells. Only after $N+1$ MMFVM steps, the cell has interacted with all the stationary cells and its value stops evolving.

The solution delay of the cell inside the spatial domain ($x=x_{\mathrm{s}}^i,i=1,\dots ,N$) is caused by a length the MMFVM solution has to travel further opposed to a true analytic solution. More precisely, the true analytic residence time ${\varsigma }_{\mathrm{T}}^i\left(x,t\right)$ is by the method of characteristics given as

$$\begin{array}{ccc}\hfill {\varsigma }_{\mathrm{T}}^i\left(x_{\mathrm{s}}^i,t\right):{\int }_{t-{\varsigma }_{\mathrm{T}}^i}^{t}v\left(s\right)\mathrm{d}s& =& x_{\mathrm{s}}^i=\left(i-\frac{1}{2}\right)\Delta x.\hfill \end{array}$$

(33)

The MMFVM solution, however, has to travel a longer path (see Figure 1) and its residence time ${\varsigma }_{\mathrm{A}}^i\left(x,t\right)$ is given by

$${\varsigma }_{\mathrm{A}}^i\left(x_{\mathrm{s}}^i,t\right):{\int }_{t-{\varsigma }_{\mathrm{A}}^i}^{t}v\left(s\right)\mathrm{d}s=x_{\mathrm{s}}^i+\frac{\Delta x}{2}=i\Delta x,$$

which can be, due to linearity of the integral operator, also written as

$${\int }_{t-{\varsigma }_{\mathrm{A}}^i}^{t-{\varsigma }_{\mathrm{A}}^i+{\varsigma }_{\mathrm{T}}^i}v\left(s\right)\mathrm{d}s+{\int }_{t-{\varsigma }_{\mathrm{A}}^i+{\varsigma }_{\mathrm{T}}^i}^{t}v\left(s\right)\mathrm{d}s=\left(i-\frac{1}{2}\right)\Delta x+\frac{\Delta x}{2}.$$

(34)

Defining a solution delay as ${\tau }^i\left(t\right)={\varsigma }_{\mathrm{A}}^i\left(t\right)-{\varsigma }_{\mathrm{T}}^i\left(t\right)$ and identifying (33) in (34) we conclude the first part of (32)

$${\int }_{t-{\tau }^i\left(t\right)}^{t}v\left(s\right)\mathrm{d}s=\frac{\Delta x}{2},i=1,\dots ,N.$$

Similarly, for $i=N+1$, the analytic residence time ${\varsigma }_{\mathrm{T}}^{N+1}$ is given by

$${\int }_{t-{\varsigma }_{\mathrm{T}}^{N+1}}^{t}v\left(s\right)\mathrm{d}s=1.$$

The MMFVM residence time ${\varsigma }_{\mathrm{A}}^{N+1}$ is given by

$$\begin{array}{ccc}\hfill {\int }_{t-{\varsigma }_{\mathrm{A}}^{N+1}}^{t}v\left(s\right)\mathrm{d}s& =& \left(N+1-\frac{1}{2}\right)\Delta x+\frac{\Delta x}{2}\hfill \\ & =& 1+\Delta x,\hfill \end{array}$$

from which we can, by the same manipulation as above, conclude the remaining part of (32)

$${\int }_{t-{\tau }^{N+1}\left(t\right)}^{t}v\left(s\right)\mathrm{d}s=\Delta x.$$

□

Corollary 2 (State jump).

When a consecutive MMFVM step is being initialized (at $p\left(t\right)=1$ from the preceding step), then the initial averaging of states (18) over the stationary grid ${\Xi }_{\mathrm{s}}$ (instead of the fully shifted ${\Xi }_{\mathrm{a}}$ grid) effectively performs a state shift/jump

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{a}}^i\left(t^{+}\right)& =& {\mathit{q}}_{\mathrm{a}}^{i-1}\left(t\right)\hfill \end{array}$$

(35a)

$$\begin{array}{ccc}\hfill {\mathit{q}}_{\mathrm{s}}^i\left(t^{+}\right)& =& {\mathit{q}}_{\mathrm{s}}^i\left(t\right)\hfill \end{array}$$

(35b)

$$\begin{array}{ccc}\hfill p\left(t^{+}\right)& =& 0,\hfill \end{array}$$

(35c)

where $t^{+}$ denotes the just-after-jump instance at time t. Note that this relation holds only for a piecewise constant spatial profile approximation. Analogous relations can be derived for higher-order reconstructions.

Proof. 

Integration of (18) at $p=1$ gives (35). See Figure 2 for illustration.

□

Note that in the case of a simple transport equation ${\partial }_t\mathit{q}+v{\partial }_x\mathit{q}=s\left(\mathit{q}\right)$, $s\left(\mathit{q}\right)=\mathit{A}\mathit{q}+\mathit{b}$, the solution delay (32) of MMFVM for the position $x_{\mathrm{s}}^{N+1}$ (corresponding to the analytic solution at $x=1$) can be compensated by a modification to the velocity $\widehat{v}=\left(N+1/N\right)v$ and modification to the source term $\widehat{s}\left(\mathit{q}\right)=N+1/N\left(\mathit{A}\mathit{q}+\mathit{b}\right)$. This statement can be shown to be true, but as it does not apply to the more general case considered, the proof is omitted.

Remark 4.

There is a minor numerical diffusion caused by the cell homogeneity (constant spatial profile over a cell). Imagine an advected cell being halfway shifted $\left(p=MM1/2\right)$, as it interacts with its left stationary cell, its state value changes and this, by cell homogeneity, in turn, affects the interaction with its right stationary cell and vice versa. This effect is, however, minor compared to FVM, where the cell homogeneity affects the greater advective flux.

In FVM, the sampling in time is bounded by the CFL condition (and possibly other stability conditions), but otherwise can be chosen freely. Sampling instances $t_k$ in MMFVM are, however, rigidly defined by the grid refinement and the actual characteristic velocity; in particular by

$${\int }_{t_0}^{t_0+t_k}v\left(t\right)\mathrm{d}t=k\Delta x,k\in \mathbb{N}.$$

Those are the instances where the stationary and advected grid, (12) and (14), coincide, and when the last advected cell is done interacting with the stationary cell, and its value remains unchanged from thereon. The advected state boundary condition is sampled at this rate. Homogeneity assumption (piecewise constant spatial profile reconstruction) disallows free time sampling at fixed spatial positions due to passing discontinuities. Free time sampling is in general possible, but the quality of the solution strongly depends on the quality of the spatial profile reconstruction.

When the characteristic velocity tends to zero, the sampling intervals tend to infinity; this is not a desired property. To prevent this scenario, the following step is advised.

Remark 5.

When the sampling interval exceeds an allowed threshold, switch the solution method to FVM, i.e., stop the advected grid and incorporate the advective flux as is usual in FVM. It may also be beneficial to momentarily refine the grid.

Note that for low velocities the relative importance of the advective flux decreases and any numerical advection artifacts of FVM decrease. Also, as the velocity increases, the solution delay of MMFVM decreases. The MMFVM method is in a sense complementary to FVM, and its advantages prevail in highly advective scenarios.

Remark 6.

For systems with no static grid, apply the boundary condition directly to the first advected cell (and the upwind cells in the case of boundary conditions on the derivatives). MMFVM reduces to a system of ODEs with state shifts and solves the problem without transport delays in the solution.

Remark 7.

The (11) can be imagined as time and spatial evolution of an infinite number of infinitesimally small slices/points. For illustration, in a cooled water pipe, we can imagine infinitely many slices of water of width dx traveling by a defined input velocity through a pipe consisting again of infinitely many, now static, slices. For increasing N, MMFVM approximates ever so closely the continuum of the governing PDE; the solution delay of Theorem 2 diminishes and so does the numerical diffusion introduced in Remark 4.

3. Application

Two examples will be given; a simple transport-decay equation, for which an analytic solution is known and a plug-flow heat exchanger, for which comparison to a fine FVM solution will validate the solution.

3.1. Advection-Decay Test

A simple transport equation with decay

$$\begin{array}{ccc}\hfill {\partial }_{t}q\left(x,t\right)+v{\partial }_{x}q\left(x,t\right)& =& -sq\left(x,t\right)\hfill \\ \hfill q\left(0,t\right)& =& b\left(t\right)\hfill \\ \hfill q\left(x,0\right)& =& 0\hfill \end{array}$$

has by the method of characteristics an analytical solution that is only a time-shifted initial or boundary condition with exponential decay along the transport,

$$q\left(x,t\right)=\left\{\begin{array}{cc}0\hfill & t<\frac{x}{v}\hfill \\ b\left(t-\frac{x}{v}\right)e^{-\frac{s}{v}x}\hfill & t\ge \frac{x}{v}.\hfill \end{array}\right.$$

This analytic solution is compared to an FVM solution with a first order upwind scheme, an FVM solution with a superbee limiter [5] and finally to the MMFVM.

According to the Theorem 2, the MMFVM solution is obtained by repetitive solution of the MMFVM step, i.e., integrate the following system in time

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}q& =& -sq\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}t}p& =& Nv\hfill \end{array}$$

for $p\in 〈0,1〉$ and perform a state jump according to Remark 2 at $p=1,$ when $b\left(t\right)$ is assigned to the first volume state.

Figure 3 shows the results for $q\left(1,t\right)$, with a simulation setting $N=5$, $v=0.1$ and $s=-vln\left(0.6\right)$. A number of elements $N=5$ is set low to highlight low-dimensional behavior of the said methods. Time integration is performed by BS23 [29] method with event location [30] (at $p=1$). The analytic solution is $q\left(1,t\right)=0.6\times b\left(t-10\right),t\ge 10.$ Figure 3 also contains the input signal $b\left(t\right)$.

It is noticeable that the FVM solution with first-order profile approximation brings great numerical diffusion to the solution. The usage of superbee slope limiter [5] improves the behavior, but under-performs in advecting discontinuities. The MMFVM solution is capable of truly advecting the input. The sampling is, however, dictated by the advection speed and is rather sparse in the case of only a few elements ($N=5$). Solution preserves the shape of input but suffers from the solution delay as analyzed in Theorem 2 and its proof.

Note that the problem contains no static states and therefore Remark 6 applies. However, to highlight properties of the method for general scenarios, the boundary condition has been applied to a first upwind ghost cell to generate the transport delay in the solution (N+1 integration steps instead of N). The source term had to be altered as $s^{\prime }=s\frac{N}{N+1}$ to account for this modification.

3.2. Water-to-Air Heat Exchanger Model

A plug flow heat exchanger is a plug flow reactor, where no reaction, only a heat exchange, occurs. Its model captures a distributed parameter, distributed time delay behavior and can be used for simplified modeling of, for example, water-to-air heat exchangers.

A water-to-air heat exchanger is considered to be a tubular single-pass cross-flow heat exchanger with the following assumptions: $\left(a_1\right)$ Heat transfer rates are constant in space, time and temperature. They solely depend on flow velocities. $\left(a_2\right)$ There is only one phase in both fluids, i.e., no condensation nor evaporation occurs.$\left(a_3\right)$ Fluids are incompressible and of constant density and specific heat capacity. $\left(a_4\right)$ Heat conduction (dispersion, diffusion) in the direction of flow is negligible. $\left(a_5\right)$ Water flow is radially coherent. The heat exchanger is thin in the air flow direction. $\left(a_6\right)$ Air flow has no dynamics and is ideally mixed at the inlet and outlet. The air temperature acting an the heat exchanger body as assumed constatnt and equal to the air inlet temperature. $\left(a_7\right)$ Water velocity (mass flow) is known.

The model is derived using the Fourier’s and conservation laws. The resulting system of partial differential equations is the first order hyperbolic

$$\begin{array}{ccc}\hfill {\partial }_t{T}_{\mathrm{w}}+v{\partial }_x{T}_{\mathrm{w}}& =& \alpha \left(T_{\mathrm{b}}-T_{\mathrm{w}}\right)\hfill \\ \hfill {\partial }_t{T}_{\mathrm{b}}& =& {\beta }_1\left(T_{\mathrm{w}}-T_{\mathrm{b}}\right)+{\beta }_2\left(T_{\mathrm{a}}-T_{\mathrm{b}}\right),\hfill \end{array}$$

where the state variables $T_{\mathrm{w}}=T_{\mathrm{w}}\left(x,t\right)\left[{}^{°}\mathrm{C}\right]$ and $T_{\mathrm{b}}=T_{\mathrm{b}}\left(x,t\right)\left[{}^{°}\mathrm{C}\right]$ denote the water and metal body temperature, respectively, and $v\left(t\right)=\dot{m}\left(t\right)/m_w\left[{\mathrm{s}}^{-1}\right]$ is a normalized transport speed, where $\dot{m}\left[\mathrm{kg}/\mathrm{s}\right]$ is a water mass flow and $m_{\mathrm{w}}\left[\mathrm{kg}\right]$ a weight of water in the heat exchanger. The space variable was normalized using $x=\xi /L\left[-\right]$, where $\xi \left[\mathrm{m}\right]$ is the space variable along the length of the heat exchanger, $L\left[\mathrm{m}\right]$ denotes the length of the heat exchanger.

The boundary conditions are

$$\begin{array}{ccc}\hfill T_{\mathrm{w}}\left(t,0\right)& =& T_{\mathrm{w},\mathrm{in}}\left(t\right)\hfill \\ \hfill T_{\mathrm{a}}\left(t,x\right)& =& T_{\mathrm{a},\mathrm{in}}\left(t\right),x\in 〈0,1〉\hfill \end{array}$$

moreover, initial conditions are

$$\begin{array}{ccc}\hfill T_{\mathrm{w}}\left(0,x\right)& =& T_{w,0}\left(x\right)\hfill \\ \hfill T_{\mathrm{b}}\left(0,x\right)& =& T_{\mathrm{b},0}\left(x\right),x\in 〈0,1〉\hfill \end{array}$$

where $T_{\mathrm{w},0}\left(x\right)$ and $T_{\mathrm{b},0}\left(x\right)$ are initial temperature profiles in the water and body, respectively, and $T_{\mathrm{w},\mathrm{in}}\left(t\right)$, $T_{\mathrm{a},\mathrm{in}}\left(t\right)\left[{}^{°}\mathrm{C}\right]$ are water and air inlet temperatures, respectively. The coefficients are

$$\alpha =\frac{H_{\mathrm{wb}}\left(\dot{m}\right)}{C_{\mathrm{w}}},{\beta }_1=\frac{H_{\mathrm{wb}}\left(\dot{m}\right)}{C_{\mathrm{b}}},{\beta }_2=\frac{H_{\mathrm{ba}}\left(\dot{V}\right)}{C_{\mathrm{b}}},$$

where $H_{\mathrm{wb}}\left(\dot{m}\right)\left[\mathrm{W}/\mathrm{K}\right]$, $H_{\mathrm{ba}}\left(\dot{V}\right)\left[\mathrm{W}/\mathrm{K}\right]$ are flow dependent heat exchange coefficients and $C_{\mathrm{w}}\left[\mathrm{J}/\mathrm{K}\right]$, $C_{\mathrm{b}}\left[\mathrm{J}/\mathrm{K}\right]$ are total heat capacities of water and metal body, respectively. Air volumetric flow is denoted $\dot{V}\left[{\mathrm{m}}^3/\mathrm{s}\right]$. The outputs of the model are a heat flow from body to air $Q\left(t\right)={\int }_0^{1}U{A}_{\mathrm{ba}}\left(T_{\mathrm{b}}\left(x,t\right)-T_{\mathrm{a},\mathrm{in}}\left(t\right)\right)\mathrm{d}x\left[\mathrm{W}\right]$ and an outlet water temperature $T_{\mathrm{w},\mathrm{out}}\left(t\right)=T_{\mathrm{w}}\left(1,t\right).$

FVM is applied in a standard way; the CFL condition is $\dot{m}\left(t\right){\tau }_{\mathrm{s}}\left(t\right)N\le m_{\mathrm{w}}$, where ${\tau }_{\mathrm{s}}\left(t\right)$ is a sampling period.

For the application of MMFVM, first define advected and static state variables $q_{\mathrm{a}}=T_{\mathrm{w}}$, $q_{\mathrm{s}}=T_{\mathrm{b}}$ and boundary conditions $b\left(t\right)=T_{\mathrm{w},\mathrm{in}}\left(t\right)$, $u\left(t\right)=T_{\mathrm{a},\mathrm{in}}\left(t\right)$. Quick observation of the source terms

$$\begin{array}{ccc}\hfill s_{\mathrm{a}}\left(q_{\mathrm{a}},q_{\mathrm{s}}\right)& =& -\alpha q_{\mathrm{a}}+\alpha q_{\mathrm{s}}\hfill \\ \hfill s_{\mathrm{s}}\left(q_{\mathrm{a}},q_{\mathrm{s}}\right)& =& {\beta }_1{q}_{\mathrm{a}}-\left({\beta }_1+{\beta }_2\right)q_{\mathrm{s}}+{\beta }_{2}u,\hfill \end{array}$$

reveals their linearity. By application of the MMFVM step (26) we get the following evolution of volume averages

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{q}_{\mathrm{a}}^i\\ q_{\mathrm{s}}^i\end{array}\right]& =& \left[\begin{array}{cc}-\alpha & \alpha \\ {\beta }_1& -\left({\beta }_1+{\beta }_2\right)\end{array}\right]\left[\begin{array}{c}{q}_{\mathrm{a}}^i\\ q_{\mathrm{s}}^i\end{array}\right]+p\left[\begin{array}{cc}0& \alpha \\ {\beta }_1& 0\end{array}\right]\left[\begin{array}{c}{q}_{\mathrm{a}}^{i-1}-q_{\mathrm{a}}^i\\ q_{\mathrm{s}}^{i+1}-q_{\mathrm{s}}^i\end{array}\right]+\left[\begin{array}{c}0\\ {\beta }_2\end{array}\right]u\hfill \end{array}$$

and a relative shift equation

$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}p& =& \frac{v}{\Delta x}.\hfill \end{array}$$

Solving now this system according to Theorem 2 we obtain the solution. Time integration is performed by BS23 [29] method with event location [30] at $p=1$.

Setting of the experiment was: $H_{\mathrm{wb}}\left(\dot{m}\right)=43.1{\left(\dot{m}/3600\right)}^{0.292}\mathrm{W}/\mathrm{K},$ $H_{\mathrm{ba}}\left(\dot{V}\right)=1.68\times {10}^{-3}{\dot{V}}^2-0.87\dot{V}+260\mathrm{W}/\mathrm{K}$, $C_{\mathrm{w}}=2.41·{10}^3\mathrm{J}/\mathrm{K}$, $C_{\mathrm{b}}=2.26\times {10}^3\mathrm{J}/\mathrm{K}$, $m_{\mathrm{w}}=0.577\mathrm{kg}$, $T_{\mathrm{w},0}=T_{\mathrm{b},0}=20{}^{°}\mathrm{C}$. The time sequence of all the four inputs is depicted in Figure 4. Figure 5 shows body and water volume temperatures at four equidistant spatial positions $x=\left[0.125,0.375,0.625,0.875\right]$ (from top to bottom). Fine FVM solution ($N=1004$) with HCUS limiter [31] mimics the role of the true solution. To have a comparison with the two consecutive solutions with $N=4$, the states of the fine solution are averaged in space over $\Delta x=1/4$. Figure 6 depicts the outputs. An interesting comparison begins with the FVM simulation with only four volumes ($N=4$) and HCUS flux limiter. Both state and output figures show how numerical diffusion alters the solution. The inlet temperature step at $t=100s$ is smeared out, whereas the true solution is a sharp moving front. MMFVM with only four volumes ($N=4$) can advect the sharp front. However, the time sampling is tightly defined by the problem setting and transport velocity, as also is the solution delay. Considerable quality enhancement is achieved by having more volumes, a low number of volumes has been chosen to intensify the distinction in properties of the compared methods. Nevertheless, the MMFVM still performs reasonably well in this highly convective scenario.

4. Conclusions

A numerical solution method for systems with dominant advective components has been presented. The method utilized the Eulerian grid for stationary states and the Lagrangian grid for states with a characteristic velocity. Formulation of the interaction of the two grids is the main contribution of the paper. The employment of the Lagrangian grid removes advection flux from the cell-to-cell flux for its volumes as derived in the proof to Theorem 1. The solution to this problem formulation can contain discontinuities and a numerical diffusion is significantly decreased. As no finite numerical method can give an absolutely accurate numerical solution, there also is a price to pay for MMFVM. The time sampling intervals are advection velocity dependent, and there is a velocity-dependent delay in the solution. It is reasonable to resort to classic FVM for low-to-zero velocities, where diffusion/dispersion/conduction begin to prevail anyway. The MMFVM method should be used as complementary to FVM in highly-advective conditions $\left(\mathrm{Pe}\gg 10\right)$. A low-order convection-accurate heat exchanger simulation model has been a primary motivation for the development of the method. The model has been implemented and its performance analyzed.