remarkRemark hypothesisHypothesis claimClaim Discrete Stream(line) MethodKiril
S. Shterev

# On convective terms approximation approach that corresponds to pure
convection††thanks: This work was funded by the Bulgarian NSF under Grant
DN-02/7-2016.

Kiril S. Shterev Institute of Mechanics, Bulgarian Academy of Sciences, Acad.
G. Bonchev Str., Block 4, Sofia 1113, Bulgaria (, www.imbm.bas.bg).
<EMAIL_ADDRESS>
###### Abstract

Recent decades are put lots of efforts to develop a higher-order scheme for
convective terms approximation that is stable and reliable. The idea presented
here is that approximation approach has to correspond to the physical
phenomenon described by approximated terms. Pure convection (advection) that
is described by convective terms is transporting a property along the
streamline, and the information propagation is unidirectional, i.e.,
transported property depends on previous values along the streamline but does
not depend on the next ones. The proposed approach represents streamlines on
mesh as discrete streamlines and is called Discrete Stream(line) Method
(DStreaM). A discrete streamline here is represented as a narrow triangle with
one vertex of the approximated node and two others neighbor upstream nodes.
Discrete streamlines are orientated using local flow direction as skew upwind
schemes. DStreaM corresponds to pure convection. Here are considered standard
test problems: advection of a step profile, advection of a double-step
profile, advection of a sinusoidal profile, and Smith and Hutton problem.
DStreaM solutions were compared with upwind-first order scheme and second-
order Total Variation Diminishing (TVD) schemes with limiters Min-Mod, QUICK,
and SUPERBEE solutions. DStreaM demonstrated second-order accuracy and rapid
convergence. Upwind and DStreaM need 2 or 4 iterations to reach a final
solution while TVD schemes need from 15 to 93.5 more iterations. DStreaM
approach looks promising for calculation of convective-dominated problems
because it approximates naturally first derivatives and is straightforwardly
applicable as a meshfree method or on unstructured meshes.

###### keywords:

Discrete Strea(line) Method, pure convection, second-order scheme, advection
of a step profile, Smith and Hutton problem

65N99, 35L60

##  1 Introduction

Development of higher-order convective terms approximation scheme is an
important problem of present-day science. Navier-Stokes (NS) equations are one
of the dominant equations used in Computational Fluid Dynamic (CFD) to model
fluid flows. NS system of equations describes unsteady, viscous,
incompressible/compressible, and heat-conductive fluid. They contain unsteady,
convective, and diffusion terms. As the diffusion terms could be easily
approximated with at least second-order spatial accuracy, convective terms can
be approximated unconditionally with upwind first-order scheme. Upwind first-
order scheme obtains correct results for all speeds and parameters, but it is
computationally expensive. The corresponding second-order scheme that obtains
correct physical results for all speeds and parameters as upwind first-order
scheme does not exist. Up to now are followed two routes to develop higher-
order convective terms approximation scheme: One way is to use higher-order
schemes that are basically derived as one-dimensional scheme and applied for
multidimensional cases and help in reducing numerical diffusion. Total
Variation Diminishing (TVD) schemes are a typical representative that follow
this route. They are designed to prevent the undesirable oscillatory behavior
of higher-order schemes. In TVD schemes, the tendency towards oscillation is
counteracted by adding an artificial diffusion fragment or by adding a
weighting towards upstream contribution. In the literature early schemes,
based on these ideas were called flux corrected transport (FCT) schemes: see
Boris and Book [2], [3]. Further, works by van Leer [29], [15], [16], [30],
Harten [10], [11], Sweby [27], Roe [22], Osher, and Chakravarthy [20], Zijlema
[35], Arora and Roe [1], Čada and Torrilhon [32], Ferreira et al. [9], Zhang
et al. [34], Kriel [14], and many others have contributed to the development
of present-day TVD schemes. Another route is to use a skew upwind scheme that
reduces numerical diffusion because of their multidimensional nature. They
take into account flow direction (for 2D and 3D flows) and determine neighbor
points according to it, in some sense, they follow the streamlines. First
Raithby presented skew upwind schemes SUDS and SUWDS in [24] and [23]. These
schemes tend to involve the neighboring diagonal cells as well as those
opposite to each cell face. Also, within the finite volume context, Hassan,
Rice, and Kim [12] proposed a mass-flow-weighted skew upwind scheme as an
improvement over the conventional skew-upwind scheme. The scheme ensures a
reduction of numerical instability and numerical diffusion errors. In the
finite element context, Schneider and Raw [25] proposed an upwind procedure
that accounts for the directionality of the flow field through a skewed
approach, while simultaneously precluding the possibility of negative
coefficients. Patel, Cross, and Markatos [21] retained the general objectives
of Raithby scheme and proposed CUPID scheme. CUPID scheme accounts for the
flow angles at the corners of a cell rather than at the cell faces, as does
Raithby’s SKEW schemes. Skew upwind weighted differencing (SUWD) scheme was
proposed by Busnaina, Zheng, and Sharif [4]. Here the convected quantity at a
cell face is approximated by a weighted average of two quantities obtained by
using linear interpolation at two locations upstream of the resultant velocity
vector at the cell face. Carey, Scanlon, and Fraser [6], [5] created SUCCA
scheme that concentrates attention at the corners of finite-volume cells.
Darwish and Moukalled [8] presented STOIC scheme (Second- and Third-Order
Interpolation for Convection) that is a high-resolution scheme developed and
implemented in the context of the normalized variable formulation methodology
developed by Leonard [17]. Ogedengbe and Naterer presented convective upwind
scheme called NISUS (Non-Inverted Skew Upwind Scheme) in [19] and [18]. The
main benefit of NISUS lies in avoiding the costly inversion of the upwind
coefficient matrix, without any significant loss of accuracy. Karadimou and
Markatos retain the objectives of CUPID and SUCCA schemes and reformulate the
convection terms in the momentum and scalar conservation equations in a way to
treat 3D flows. They presented SUPER (Skew Upwind and Corner Algorithm) scheme
in [13]. After all, both routes obtain more accurate results than first-order
upwind scheme, but they can give rise of oscillations or under/over-shoots,
especially in regions of strong gradients when are applied to calculate fast
compressible flows. TVD schemes particularly obtain excellent results for the
steady state, slow and moderate flows. Slow and moderate flows are unsteady
but close to steady one for lower Reynolds numbers. When unsteady fluid flows
is fast TVD schemes are not time accurate, see Chung [7]. Furthermore, if
Reynolds number is higher, even when fluid flow is steady, TVD schemes can
obtain an unphysical oscillated solution. In such cases, the solution can be
obtained using the first-order upwind scheme using significant computational
resources. Van Leer in his review article [31] writes that American
aeronautical community is skeptical about the cutting-edge algorithm, and they
prefer the advent of high-performance computing and promise of massively
parallel computing instead of the development of higher order approximation
schemes. The need for a new approach is evident.
The approximation scheme should correspond to process that approximated
term(s) describe. Diffusion terms describe exchange of concentration of
specific property from one place to another. Diffusion depends on
concentration and does not depend on flow direction, i.e., does not depend on
velocity. Information is exchanged in all directions according to a specific
node. One of the simplest approximation schemes for diffusion terms is the
central difference, it is unconditionally stable when approximate diffusion
terms, it is a second-order scheme, and it takes information from neighbor
nodes in all directions. The convection opposite to diffusion depends on flow
direction and does not depend on concentration. Convection terms describe
transportation of property ($\phi$) along the streamline. They contain
information about streamline (velocities in convective terms) and previous
values of $\phi$ (see equation (1)) along the streamline, but they do not
contain information for $\phi$ further along the streamline. Approximation
scheme that corresponds to convection has to contain upwind information along
the streamline that passes through the approximated node and no downstream
information about $\phi$. TVD schemes use downstream information about $\phi$
that corresponds to diffusion, which is not included as information in
convective terms. This is their problem and the reason for their limited
application.
In this paper is presented approximation approach for convective terms that
use only upstream points and orients approximation scheme along streamlines
that makes it of a group of skew upwind schemes. Four standard two-dimensional
pure convection (advection) steady test cases were used to compare the
relative performance of proposed Discrete Stream(line) Method (DStreaM) with
standard schemes: upwind first-order and second-order TVD schemes. TVD schemes
are well studied and widely used in Computational Fluid Dynamic (CFD) codes
that make them appropriate schemes for presented comparisons. Here are used
TVD scheme with limiters Min-Mod [22], QUICK [17] and SUPERBEE [22]. Sweby
[27] has given necessary and sufficient conditions for a scheme to be TVD
second-order. The Min-Mod limiter function exactly traces the lower limit of
the TVD region, SUPERBEE limiter follows the upper limit, and QUICK limiter is
between them. The presented test problems are advection of a step profile,
advection of a double-step profile, advection of a sinusoidal profile, and
Smith and Hutton problem. Obtained results, show the second-order spatial
accuracy of DStreaM and significant reduction of iterations from 15 to 93.5
times according to TVD schemes.

##  2 Considerations on convective and diffusion terms information propagation

Here are presented more detailed considerations about the pure convection and
diffusion processes that are important for further explanations of the
proposed approximation approach.
The steady convection-diffusion equation for a general property $\phi$ is:

(1) $\underbrace{\rho\frac{\partial(u\phi)}{\partial
x}+\rho\frac{\partial(v\phi)}{\partial y}}_{\text{Convective
terms}}=\underbrace{\Gamma\left(\frac{\partial^{2}\phi}{\partial
x^{2}}+\frac{\partial^{2}\phi}{\partial y^{2}}\right)}_{\text{Diffusion
terms}},$

where $u$ and $v$ are velocities along x- and y-axis, $\rho$ is density, and
$\Gamma$ is diffusion coefficient.
The diffusion terms (see equation (1)) describe pure diffusion. Diffusion is
the movement of molecules and/or atoms from a region of high concentration (or
high chemical potential) to a region of low concentration (or low chemical
potential). The information propagates in all directions independently of the
velocity field. Fig. 1 (a) shows a diagram of directions that information
propagates. Green continuous arrows denote the information that propagates to
$\phi$ while red dashed arrows denote the information that propagates outward
$\phi$. The information propagates in both directions. Therefore, the
diffusion terms approximation scheme should include neighbor nodes in all
directions. The central difference numerical scheme is second-order, matches
perfectly to the diffusion terms properties, and approximates diffusion terms
without need of any additional criteria to ensure convergence or prevent a
non-physical solution.
Fig. 1 (b) shows the information propagation of convective-dominated fluid
flow schematically. The information about $\phi$ is available inside
information propagation cone, but it is not available outside the cone
including opposite or normal to velocity directions.

(a)

(b)

Figure 1: Schematic information propagation of diffusion dominated problem (a)
and convection dominated problem (b).

The convective terms (see equation (1)) describe pure convection (advection).
Advection can be considered as a motion of fluid particles through a
streamtubes without exchange of properties between neighbor streamtubes. The
transported fluid property is $\rho\phi$, and velocities in convective terms
can be interpreted as fluid particle’s streamline/path geometry, see Navier-
Stokes system of equations derivation and considerations [28], [33]. Fig. 2
shows the motion of a fluid particle from time $t^{0}$ to time $t^{n}$ through
a streamtube. The velocity field is known in all directions according to
position of the fluid particle while the information relates transported
property at time $t^{n}$ exist only in previous times of its travel up to
$t^{0}$ and does not exist for next times as $t^{n+1}$ because they do not
happen yet, and this kind of information is not included in convective terms.
The unidirectional information propagation of transported property is the main
characteristic of convective terms. When it is considered pure convection
case, transported property of fluid particle at time $t^{n}$ is the same as at
time $t^{n-1}$. Therefore, they are the same as transported property at time
$t^{0}$: $\rho^{n}\phi^{n}=\rho^{n-1}\phi^{n-1}=...=\rho^{0}\phi^{0}$. When
the density is constant, it follows $\phi^{n}=\phi^{n-1}=...=\phi^{0}$. That
makes possible to determine an exact solution of a steady pure convection
problem at first glance. The basic idea of the proposed approach uses
convective terms unidirectional information transportation characteristic.

Figure 2: Transport of fluid particle through streamtube.

The propagation of information of convective-dominated problems in mesh is
important to develop approximation that corresponds to modeled flow. Fig. 3
shows propagation information cones of convective terms applied in nodes of a
Cartesian uniform mesh. In ideal pure convection case each cone’s angle
$\alpha$ (see Fig. 1, (b)) has to approach 0 ($\alpha\rightarrow 0$). Thus,
will turn each cone into a curve and the node’s value will be equal to value
at the other side of the curve. That can be another node if the curve goes
through it or the computational domain boundary. $\alpha\rightarrow 0$
corresponds to numerical scheme without false diffusion. False diffusion is
information propagation normal to the streamline. Here cones represent
information propagation of numerical solution of pure convection problems
where false diffusion exists ($\alpha>0$). The cones, in Fig. 3, can be
separated into three groups according to the information that reaches
$\phi_{i,j}$. The information from the first group reaches $\phi_{i,j}$. In
considered case, these are $\phi_{i-1,j}$ and $\phi_{i-2,j}$. The information
from the second group could reach $\phi_{i,j}$. It depends on the cone’s angle
($\alpha$). The nodes of this group are $\phi_{i-2,j-1}$, $\phi_{i-2,j+1}$,
$\phi_{i-1,j-1}$, and $\phi_{i-1,j+1}$. In present case information from node
$\phi_{i-2,j-1}$ reaches $\phi_{i,j}$ while information from nodes
$\phi_{i-2,j+1}$, $\phi_{i-1,j-1}$, and $\phi_{i-1,j+1}$ does not. And the
information from the last group where are the nodes $\phi_{i,j-1}$,
$\phi_{i,j+1}$, $\phi_{i+1,j-1}$, $\phi_{i+1,j}$, and $\phi_{i+1,j+1}$ does
not reach $\phi_{i,j}$. The information propagation cone of $\phi_{i,j}$
rotated on 180 degrees is reversed information propagation cone and includes
the nodes that information propagates to $\phi_{i,j}$, see Fig. 4. According
to these considerations, the approximation scheme of convective terms for the
node $\phi_{i,j}$ can include only nodes that information reaches it and could
not include nodes that information is impossible to achieve it.

Figure 3: Information propagation cones of convection dominated problem at
mesh nodes.

##  3 Discrete Stream(line) Method (DStreaM)

Two-dimensional problem is considered here. DStreaM can be extended to a
three-dimensional problem straightforwardly.
Advection equation (pure convection) can be derived from convection-diffusion
equation (1) when $\Gamma=0$. For simplification, density is constant
($\rho=const$):

(2) $\rho\frac{\partial(u\phi)}{\partial
x}+\rho\frac{\partial(v\phi)}{\partial y}=0,$

DStreaM is applied on Cartesian uniform mesh for simplicity, but it can be
applied as a meshless method or on an unstructured mesh. Fig. 4 shows that
information from $\phi_{i-1,j}$, $\phi_{i-2,j}$ and $\phi_{i-2,j-1}$
propagates to $\phi_{i,j}$. The nodes $\phi_{i-2,j}$, $\phi_{i-1,j}$ and
$\phi_{i,j}$ lie on a straight line and the information that propagates to
$\phi_{i,j}$ is information from $\phi_{i-1,j}$. The information from
$\phi_{i-2,j}$ does not reach $\phi_{i,j}$, because it is behind
$\phi_{i-1,j}$. Therefore, the shape function interpolates nodes $\phi_{i,j}$,
$\phi_{i-1,j}$ and $\phi_{i-2,j-1}$, see Fig. 4. Convective terms shape
function is a triangle with vertexes $\phi_{i,j}$, $\phi_{i-1,j}$ and
$\phi_{i-2,j-1}$. It can be considered as a discrete representation of
streamline through $\phi_{i,j}$, and as it is the basic element that
represents the idea of the proposed approach it is included in the approach‘s
name. Shape function based on three nodes is:

(3) $f^{shape}(x,y)=c_{0}+c_{1}x+c_{2}y$

Figure 4: Reversed information propagation cone and shape function that
corresponds to convective terms, and it is a discrete representation of a
streamline through point $\phi_{i,j}$.

Seeking more general form shape function nodes are denoted as $\phi_{0}$,
$\phi_{1}$, and $\phi_{2}$ (see Fig 5 (a)) and correspond to $\phi_{i,j}$,
$\phi_{i-1,j}$, and $\phi_{i-2,j-1}$ in Fig. 4, respectively. Coefficients
$c_{0}$, $c_{1}$, and $c_{2}$ are determined using condition that shape
function interpolates nodes $\phi_{0}$, $\phi_{1}$, and $\phi_{2}$ with
coordinates $(x_{0},y_{0})$, $(x_{1},y_{1})$, and $(x_{2},y_{2})$,
respectively. The coefficients are equal to:

(4)
$\begin{split}c_{0}=&\dfrac{\phi_{0}(x_{2}y_{1}-x_{1}y_{2})+\phi_{1}(x_{0}y_{2}-x_{2}y_{0})+\phi_{2}(x_{1}y_{0}-x_{0}y_{1})}{x_{0}(y_{2}-y_{1})+x_{1}(y_{0}-y_{2})+x_{2}(y_{1}-y_{0})}\\\
c_{1}=&\dfrac{\phi_{0}(y_{2}-y_{1})+\phi_{1}(y_{0}-y_{2})+\phi_{2}(y_{1}-y_{0})}{x_{0}(y_{2}-y_{1})+x_{1}(y_{0}-y_{2})+x_{2}(y_{1}-y_{0})}\\\
c_{2}=&-\dfrac{\phi_{0}(x_{2}-x_{1})+\phi_{1}(x_{0}-x_{2})+\phi_{2}(x_{1}-x_{0})}{x_{0}(y_{2}-y_{1})+x_{1}(y_{0}-y_{2})+x_{2}(y_{1}-y_{0})}\end{split}$

(a) Shape function with average velocities.

(b) Horizontal average velocity approximation ($u^{average}$).

(c) Vectors $\overrightarrow{V}^{average}$ and
$\overrightarrow{\phi_{0}\phi_{M}}$ that are used to determine sector for
approximation scheme.

Figure 5: Shape function with average velocities (a), average velocities idea
(b), and vectors that determine sector for approximation scheme (c).

The velocities are approximated as an average value determined as follows:

(5) $\begin{split}u^{average}=(2u_{i,j}+u_{i-1,j}+u_{i-2,j-1})/4\\\
v^{average}=(2v_{i,j}+v_{i-1,j}+v_{i-2,j-1})/4\end{split}$

Fig. 5 (b) shows $u^{average}$ geometrical representation. Control volume and
shape function coincide. $u^{average}$ and $v^{average}$ are considered as a
constant when integrating over control volume. The numerical equation is
derived substituting shape function (3) and average velocities (5) in (2) and
integrating over control volume:

(6) $\begin{split}&\int_{CV}\rho\frac{\partial(u\phi)}{\partial
x}+\rho\frac{\partial(v\phi)}{\partial y}\,dV\approx\\\ &\rho
u^{average}\int_{CV}\frac{\partial(f^{shape}(x,y))}{\partial x}\,dV+\rho
v^{average}\int_{CV}\frac{\partial(f^{shape}(x,y))}{\partial y}\,dV=\\\ &\rho
u^{average}c_{1}V+\rho v^{average}c_{2}V=0\end{split}$

After substitution and simplification the numerical equation for $\phi_{0}$
can be expressed as:

(7) $\begin{split}\phi_{0}=(a_{1}\phi_{1}+a_{1}\phi_{2})/a_{0},\end{split}$

where:

(8) $\begin{split}a_{0}=&a_{1}+a_{2}\\\
a_{1}=&u^{average}(y_{0}-y_{2})-v^{average}(x_{0}-x_{2})\\\
a_{2}=&v^{average}(x_{0}-x_{1})-u^{average}(y_{0}-y_{1})\end{split}$

The scheme approximates values between nodes in a set range according to the
approximated node ($\phi_{i,j}\equiv\phi_{0}$). The range applied on uniform a
Cartesian mesh is a positive integer number. The nodes taken into account are
on maximal distance from central node $\pm r\Delta$, where $r$ is a range,
$\Delta$ is a spatial step equal to
$\Delta=x_{i}-x_{i-1}=x_{i+1}-x_{i}=y_{j}-y_{j-1}=y_{j+1}-y_{j}$. Fig. 6 shows
shape functions for ranges 1, 2, and 3 according to node $\phi_{i,j}$. One can
see that shape functions correspond to possible discretized streamlines
through node $\phi_{i,j}$ for the considered ranges. Neighbor nodes taken into
account are with indexes from i-range to i+range along x-axis and from j-range
to j+range along y-axis. As example when range is equal to 1 (Fig. 6 (a)),
approximation scheme consider all nodes within a range $\pm\Delta$ that
indexes along the x- and the y-axis are i-1,i,i+1 and j-1,j,j+1, respectively.
As a result, shape functions separate space into 8 sectors. When range
increases with 1, every shape function is separated into two, and a node is
added. Thus, node numbers for ranges 1, 2, and 3 are 8, 16, and 32,
respectively (see Fig. 6). All angles with vertex node $\phi_{i,j}$ for range
1 are equal. For other ranges, not all angles are equal. The sector and
corresponding nodes $\phi_{1}$ and $\phi_{2}$ used in approximation are
selected using criteria. The selected sector has a maximum angle between
average velocity vector ($\overrightarrow{V}^{average}$) and median vector
$\overrightarrow{\phi_{0}\phi_{M}}$ (see Fig. 5 (c)). In this way, the average
velocity of selected shape function is maximum closer to median line
$\phi_{0}\phi_{M}$, and the flow direction is from nodes $\phi_{1}$ and
$\phi_{2}$ toward $\phi_{0}$ that corresponds to uniform information
propagation of convective terms. The shape function represents the mapping of
local streamline into the mesh as discrete streamline. When range increases,
discrete streamline become thinner and thinner. Thinner discrete streamline
tends to a line, and it is closer to a local streamline when local velocity
field is uniform. When local velocity field is not uniform, thinner discrete
streamlines deviate from local streamlines. Note that average velocity could
be different for different sectors when velocity field is not uniform and has
to be taken into account in algorithm when the sector is selected.

(a) Shape functions for Range=1.

(b) Shape functions for Range=2.

(c) Shape functions for Range=3.

Figure 6: Shape functions according to node $\phi_{i,j}$ for ranges equal to 1
a), 2 b), and 3 c).

Proposed approximation approach corresponds to straight local streamline that
makes it simple, but some limitations occur when the velocity field is not
uniform, and the range is large. DStreaM obtains solution without under/over-
shoots when velocity field is uniform, and the range is large. These test
cases are advection of a step profile, advection of a double-step profile, and
advection of a sinusoidal profile presented further in the paper. When the
local streamlines are curved, the considered approximation shows its
limitations. Such approximation is acceptable for short range near the
approximated node. The last considered test case is Smith and Hutton problem
where local streamlines are curved, see Fig. 10. DStreaM obtains excellent
results without under/over-shoots for ranges up to 3, see Fig. 11. When the
range is equal or higher than 4, obtained results have under- and over-shoots.
A standard approach to prevent solution under/over-shoots and to enhance
convergence is to keep coefficients $a_{1}$ and $a_{2}$ always positive, see
Patel, Cross, and Markatos [21] recommendations. Here a limiter can be used,
as follows:

(9)
$\begin{split}a_{1}=&max(0,u^{average}(y_{0}-y_{2})-v^{average}(x_{0}-x_{2}))\\\
a_{2}=&max(0,v^{average}(x_{0}-x_{1})-u^{average}(y_{0}-y_{1}))\end{split}$

Using limiting coefficients when $a_{0}=a_{1}+a_{2}$ means that $\phi_{0}$
(see equation (7)) will be between values of $\phi_{1}$ and $\phi_{2}$,
because their coefficients $a_{1}/a_{0}$ and $a_{2}/a_{0}$ are always
positive, they are between 0 and 1, and their sum is equal to 1
($a_{1}/a_{0}+a_{2}/a_{0}=1$).
Limiting coefficients prevent approximation scheme to obtain results with
under/over-shoots and increases scheme ranges that obtain physically realistic
results. Fig. 11 shows obtained fields with non-limited coefficients (8), Fig.
11 (a) - (e), and limited coefficients (9), Fig. 11 (g) - (k). Both obtain the
same results up-to range 3. Differences occur when the range is equal to or
greater than 4. For range 4 scheme with non-limited coefficients obtains a
result with small under/over-shoots, while scheme with limited coefficients
obtains a result without under/over-shoots, see Fig. 11 (e) and (j),
respectively. For range 5 result obtained with non-limited coefficients is
with large under- and over-shoots, Fig. 11 (e). The solution obtained using
limited coefficients, and range 5 is without under/over-shoots, but it has
unphysical jumps between some neighbor nodes, Fig. 11 (k). The maximal range
for this spatial step that obtains correct physical results is 4 when are used
limited coefficients. After all, limited coefficients make the scheme more
stable and increase the scheme’s application range.
DStreaM applies to nodes with an arbitrary position near $\phi_{0}$ that makes
is straightforwardly applicable as a meshfree method or on unstructured
meshes.
DStreaM consists of the following steps:
1\. Calculate average velocities $u^{average}$ and $v^{average}$ (5).
2\. Determine upwind nodes $\phi_{1}$ and $\phi_{2}$.
3\. Calculate $\phi_{0}$ (7).

##  4 Applications

DStreaM was tested on standard advection problems: advection of a step
profile, advection of a double-step profile, advection of a sinusoidal
profile, and Smith and Hutton problem. Numerical equations were calculated
using Gauss-Seidel iterative method with initial guess $\phi=0$. Calculation
sequence of mesh nodes influences the number of iterations that each algorithm
needs to reach the solution. Here were used two approaches to obtain a
numerical solution with a minimal number of iterations for each scheme type.
TVD schemes start every iteration from corner $i=1,j=1$. First-order upwind
scheme and DStreaM change starting corner at each iteration sequentially.
The following schemes obtained compared numerical results: first-order upwind
scheme, DStreaM with ranges from 1 to 5 (DStreaM R1 to R5), TVD scheme with
limiters Min-Mod, QUICK, and SUPERBEE. TVD schemes used appropriate under-
relaxation coefficients to reach a solution with minimum iterations.
DStreaM can obtain an exact solution for uniform velocity field when it uses a
sufficiently large range. For presented test problems with uniform velocity
field when horizontal and vertical velocities are equal ($u=v$) DStreaM
obtains an exact solution for range equal to 1. These test problems are
advection of a step profile, advection of a double-step profile, and advection
of a sinusoidal profile. To prevent obtaining exact solution were used
horizontal and vertical components of velocity equal to $u=0.8$ and $v=1$,
respectively.

###  4.1 Advection of a step profile

Horizontal and vertical components of velocity are equal to $u=0.8$ and $v=1$,
respectively. The boundary conditions are rather standard for this problem,
see Fig. 7 (a). The mesh is uniform Cartesian mesh with $30\times 30$ nodes.

(a) Physical domain

(b) $\phi$ profile at $y=0.8$

Figure 7: Advection of a step profile: (a) physical domain and (b) obtained
profiles at $y=0.8$.
Table 1: Under-relaxation coefficient, number of iterations, and obtained a maximum residual of considered schemes for calculation of advection of a step profile with convergence criteria $\epsilon=10^{-8}$. Approximation scheme | $\alpha$ | Iterations | max residual
---|---|---|---
Upwind | 1.0 | 2 | 0
DStreaM R1 | 1.0 | 2 | 0
DStreaM R2 | 1.0 | 2 | 0
DStreaM R3 | 1.0 | 2 | 0
DStreaM R4 | 1.0 | 2 | 0
DStreaM R5 | 1.0 | 2 | 0
Min-Mod | 0.95 | 36 | $8.1\times 10^{-9}$
QUICK(TVD) | 0.95 | 49 | $6.2\times 10^{-9}$
SUPERBEE | 0.95 | 151 | $9.4\times 10^{-9}$

DStreaM with range 1 (DStreaM R1) solution is between TVD scheme with Min-Mod
and QUICK limiters solutions, see Fig. 7 (b). DStreaM R2 solution is close to
QUICK while DStreaM R3 solution is close to SUPERBEE. DStreaM R4 obtains
better solution than SUPERBEE. DStreaM R5 obtains the exact solution. DStreaM
does not need under-relaxation coefficient and obtains the solution for 2
iterations as upwind first-order scheme. Upwind and DStreaM need 18, 24.7, and
75.5 times fewer iterations than Min-Mod, QUICK, and SUPERBEE, see Table 1.

###  4.2 Advection of a double-step profile

Horizontal and vertical components of velocity are equal to $u=0.8$ and $v=1$,
respectively. The boundary conditions are rather standard for this problem,
see Fig. 8 (a). The mesh is uniform Cartesian mesh with $30\times 30$ nodes.

(a) Physical domain

(b) $\phi$ profile at $y=0.8$

Figure 8: Advection of a double-step profile: (a) physical domain and (b)
obtained profiles at $y=0.8$.
Table 2: Under-relaxation coefficient, number of iterations, and obtained a maximum residual of considered schemes for calculation of advection of a double-step profile with convergence criteria $\epsilon=10^{-8}$. Approximation scheme | $\alpha$ | Iterations | max residual
---|---|---|---
Upwind | 1.0 | 2 | 0
DStreaM R1 | 1.0 | 2 | 0
DStreaM R2 | 1.0 | 2 | 0
DStreaM R3 | 1.0 | 2 | 0
DStreaM R4 | 1.0 | 2 | 0
DStreaM R5 | 1.0 | 2 | 0
Min-Mod | 1.0 | 30 | $6.7\times 10^{-9}$
QUICK(TVD) | 0.95 | 55 | $5.6\times 10^{-9}$
SUPERBEE | 0.9 | 187 | $9.9\times 10^{-9}$

DStreaM R1 solution is between solutions obtained by TVD scheme with QUICK and
SUPERBEE limiters, see Fig. 8 (b). DStreaM R3 solution is close to SUPERBEE.
DStreaM R4 obtains better solution than SUPERBEE. DStreaM R5 obtains the exact
solution. DStreaM does not need under-relaxation coefficient and obtains the
solution for 2 iterations as upwind first-order scheme. Upwind and DStreaM
need 15, 27.5, and 93.5 times fewer iterations than Min-Mod, QUICK, and
SUPERBEE, see Table 2.

###  4.3 Advection of a sinusoidal profile

Horizontal and vertical components of velocity are equal to $u=0.8$ and $v=1$,
respectively. The boundary conditions are rather standard for this problem,
see Fig. 9 (a). The mesh is uniform Cartesian mesh with $30\times 30$ nodes.

(a) Physical domain

(b) $\phi$ profile at $y=0.8$

Figure 9: Advection of a sinusoidal profile: (a) physical domain and (b)
obtained profiles at $y=0.8$.
Table 3: Under-relaxation coefficient, number of iterations, and obtained a maximum residual of considered schemes for calculation of advection of a sinusoidal profile with convergence criteria $\epsilon=10^{-8}$. Approximation scheme | $\alpha$ | Iterations | max residual
---|---|---|---
Upwind | 1.0 | 2 | 0
DStreaM R1 | 1.0 | 2 | 0
DStreaM R2 | 1.0 | 2 | 0
DStreaM R3 | 1.0 | 2 | 0
DStreaM R4 | 1.0 | 2 | 0
DStreaM R5 | 1.0 | 2 | 0
Min-Mod | 0.95 | 30 | $6.6\times 10^{-9}$
QUICK(TVD) | 0.9 | 51 | $7.9\times 10^{-9}$
SUPERBEE | 0.85 | 65 | $8.3\times 10^{-9}$

DStreaM R1 solution corresponds to QUICK’s solution, see Fig. 8 (b). DStreaM
R2, R3, and R4 solutions are closer to the exact solution than QUICK’s and
SUPERBEE’s solutions. DStreaM R5 obtains the exact solution. DStreaM scheme
does not need under-relaxation coefficient and obtains a solution for 2
iterations as upwind first-order scheme. Upwind and DStreaM need 15, 25.5, and
32.5 times fewer iterations than Min-Mod, QUICK, and SUPERBEE, see Table 3.

###  4.4 Smith and Hutton problem

The Smith and Hutton problem [26] will be investigated at $u=2y(1-x^{2})$ and
$v=-2x(1-y^{2})$. Along the inlet $\phi$ is prescribed by the distribution of
$\phi(-1\leq x\leq 0,y=0)=1+\text{tanh}[\alpha(2x+1)]$, see Fig. 10. Along the
lines $x=-1$, $y=1$, and $x=1$, $\phi$ is prescribed as
$1-\text{tanh}(\alpha)$ while along the outlet $(0\leq x\leq 1,y=0)$ a zero
gradient is specified for $\phi$. The considered case is for $\alpha=1000$ and
$\Gamma=0$. Detailed comparisons are presented for mesh with $40\times 20$
nodes. After that, approximation schemes mesh convergence is presented.

Figure 10: Smith and Hutton problem: boundary conditions, velocity equations,
and streamlines.

As was considered earlier DStreaM with limited coefficients (9) obtains
results without under- or over-shoots (Fig. 11 (g)-(k)) while DStreaM with
non-limited coefficients (8) obtains results with under/over-shoots when the
range is equal to or greater than 4, (Fig. 11 (a)-(e)). The reason is that
local streamline through $\phi_{i,j}$ is curved and deviate from discrete
streamline when the range is larger because when the range is larger discrete
streamline approaches a straight line. In general, a curved line can be
approximated within sufficient accuracy using small straight lines
successfully. When we use the smaller spatial step as a case with $320\times
160$ nodes, obtained results are without under/over-shoots for greater ranges
than 4, see Fig. 13 (a)-(f). The reason is that the range 5 for mesh with
$320\times 160$ nodes corresponds to distance $5\Delta_{320\times
160}=5*1/160=0.03125$ ($\Delta_{320\times 160}$ is the spatial step for mesh
with $320\times 160$ nodes) that is smaller than distance for mesh with
$40\times 20$ nodes that is $5\Delta_{40\times 20}=5*1/20=0.25$
($\Delta_{40\times 20}$ is the spatial step for mesh with $40\times 20$
nodes). Finer mesh makes discrete streamline length shorter and increases
their number. That lead to a better approximation of curved streamline (non-
uniform velocity field). That property of DStreaM makes it suitable for
applications on unstructured meshes where smaller spatial steps are used in
places of high gradients.
The comparison shows that DStreaM spatial accuracy is second order as TVD
schemes and requires 4 iterations as upwind scheme. Fig. 11 shows comparisons
of obtained fields by DStreaM, upwind and TVD schemes with limiters Min-Mod,
QUICK, and SUPERBEE. DStreaM results with ranges from 1 to 3 (Fig. 11 (g),(h),
and (i)) correspond to results of TVD schemes with limiters Min-Mod, QUICK,
and SUPERBEE (Fig. 11 (l), (m), and (n)), respectively. Profiles of obtained
solutions along the line $(0\leq x\leq 1,y=0)$ are plotted in Fig. 12. The
comparison shows second order accuracy of the proposed approach DStreaM.
Upwind and DStreaM obtain a final solution for 4 iterations while Min-Mod and
QUICK need 81 and 62 iterations, see Table 4. TVD schemes require
approximately 20 and 15 times more iterations to obtain the same solution as
DStreaM. SUPERBEE scheme cannot obtain a solution within considered residual
despite small under-relaxation coefficient and many iterations while DStreaM
R3 can obtain an as sharp solution as SUPERBEE for four iterations.

(a) DStreaM R1,

non-limited coeff. (8)

(b) DStreaM R2,

non-limited coeff. (8)

(c) DStreaM R3,

non-limited coeff. (8)

(d) DStreaM R4,

non-limited coeff. (8)

(e) DStreaM R5,

non-limited coeff. (8)

(f) Upwind

(g) DStreaM R1,

limited coeff. (9)

(h) DStreaM R2,

limited coeff. (9)

(i) DStreaM R3,

limited coeff. (9)

(j) DStreaM R4,

limited coeff. (9)

(k) DStreaM R5,

limited coeff. (9)

(l) Min-Mod

(m) QUICK(TVD)

(n) SUPERBEE

Figure 11: Smith and Hutton problem contour plots of $\phi$ on a mesh with
$40\times 20$ nodes for $\alpha=1000$ and $\Gamma=0$. Results were obtained
using (a) - (e) DStreaM with non-limited coefficients (8) with ranges from 1
to 5, respectively, upwind first-order scheme (f), (g) - (k) DStreaM with
limited coefficients (9) with ranges from 1 to 5, respectively, and TVD
schemes with limiters Min-Mode (l), QUICK (m), and SUPERBEE (n).
Figure 12: Obtained profiles at the outlet ($0\leq x\leq 1$ and $y=0$) for
Smith and Hutton problem, $\alpha=1000$, $\Gamma=0$.

(a) Upwind

(b) DStreaM R1

(c) DStreaM R2

(d) DStreaM R3

(e) DStreaM R4

(f) DStreaM R5

(g) Min-Mod

(h) QUICK(TVD)

(i) SUPERBEE

Figure 13: Smith and Hutton problem contour plots of $\phi$ on a mesh with $320\times 160$ nodes for $\alpha=1000$ and $\Gamma=0$. Results were obtained using upwind first-order scheme (a), (b) - (f) DStreaM with limited coefficients (9) with ranges from 1 to 5, respectively, and TVD schemes with limiters Min-Mode (g), QUICK (h), and SUPERBEE (i). Table 4: Under-relaxation coefficient, number of iterations, and obtained a maximum residual of considered schemes for calculation of Smith and Hutton problem with convergence criteria $\epsilon=10^{-8}$ on a mesh with $40\times 20$ nodes . Approximation scheme | $\alpha$ | Iterations | max residual
---|---|---|---
Upwind | 1.0 | 4 | 0
DStreaM R1 | 1.0 | 4 | 0
DStreaM R2 | 1.0 | 4 | 0
DStreaM R3 | 1.0 | 4 | 0
DStreaM R4 | 1.0 | 4 | 0
DStreaM R5 | 1.0 | 4 | 0
Min-Mod | 0.95 | 81 | $5.5\times 10^{-9}$
QUICK(TVD) | 0.85 | 62 | $9.0\times 10^{-9}$
SUPERBEE | 0.1 | 50000 | $10^{-3}$

Mesh convergence study shows DStreaM advantages. Fig. 14 (a) shows a maximal
difference between the numerical and exact solution of obtained fields as a
function of the number of mesh nodes. Presented solutions obtained by DStreaM
are without unphysical oscillations (when limited coefficients are not
needed). All schemes reduced obtained maximal difference for finer mesh as TVD
schemes demonstrate better mesh convergence compared to DStreaM with a fixed
range. Figures 11 and 13 also show DStreaM accuracy "shifting". DStreaM can
use larger ranges on a finer mesh to reduce maximal difference and to keep its
spacial accuracy closer to SUPERBEE. Fig. 14 (b) presents a number of
iterations needed to obtain a solution with convergence criteria
$\epsilon=10^{-8}$. TVD schemes required under-relaxation coefficient to
converge iterative solver. As a number of iterations vary according to under-
relaxation coefficients, results presented in Fig. 14 (b) for TVD schemes is
minimally obtained for a range of under-relaxation coefficients. SUPERBEE
obtains the solution within required convergence criteria only for mesh with
$20\times 10$ nodes for 68 iterations while for the others meshes it reaches
maximum residual $\sim 10^{-3}$, and the number of iterations for SUPERBEE are
not presented. On the other hand, upwind and DStreaM do not require under-
relaxation coefficients to ensure convergence of iterative process. TVD
schemes require more iterations on finer mesh and number of iterations rapidly
increases from 40 to 1700. Contrary, upwind and DStreaM need 4 iterations that
are independent from mesh refinement. The results show that TVD schemes
demonstrate better mesh convergent than DStreaM with a fixed range, but
DStreaM can increase ranges for finer mesh and to obtain results with
corresponding spatial accuracy to SUPERBEE, see Fig. 13. After all, DStreaM
needs 4 iterations to obtain a final solution that is from 10 to 425 times
fewer iterations compared to TVD schemes, and it can obtain an as sharp
solution as SUPERBEE within a given residual when range varies.

(a) Mesh convergence.

(b) Number of iterations as function of number of mesh nodes.

Figure 14: Smith and Hutton problem: (a) mesh convergence and (b) number of
iterations as function of number of mesh nodes.

##  5 Conclusions

In this paper was proposed Discrete Stream(line) Method (DStreaM) for
convective terms approximation. The main idea is the approximation approach to
corresponds to a physical phenomenon described by approximated terms. In
considered case, the physical phenomenon is pure convection that is described
by convective terms. Convective terms describe transporting a property along
the streamline, and the main characteristics are unidirectional information
propagation without exchange information with transported properties over
neighbor streamlines. DStreaM approach maps streamlines to a mesh. As a
result, streamlines are represented as Discrete Streamlines. Discrete
Streamlines represent unidirectional information propagation in mesh
orientated according to velocity field. Four standard advection (pure
convection) test problems were presented in the paper: advection of a step
profile, advection of a double-step profile, advection of a sinusoidal
profile, and Smith and Hutton problem. DStreaM’s results were compared with
results obtained by upwind first-order scheme and TVD schemes with limiters
Min-Mod, QUICK, and SUPERBEE. TVD schemes use under-relaxation coefficients to
ensure convergence and require from 15 to 93.5 times more iterations to obtain
its final solution than upwind and DStreaM. For Smith and Hutton problem
SUPERBEE cannot obtain a solution within given convergence criteria, DStreaM
can obtain an as sharp solution as SUPERBEE within 4 iterations. DStreaM
combines the positives of upwind and TVD schemes. It obtains a solution within
2 or 4 iterations without under-relaxation coefficients as upwind, and it is
second order scheme as TVD schemes.
DStreaM approach looks promising for calculation of convective-dominated
problems because:

  * •

Naturally approximates first derivatives in partial differential equations;

  * •

All independent variables can be defined in the same nodes;

  * •

It can be straightforward and easy applied on arbitrary meshes as a meshfree
method or unstructured meshes with real-time mesh refinement;

## Acknowledgments

We would like to acknowledge the financial support provided by the Bulgarian
NSF under Grant DN-02/7-2016.

## References

  * [1] M. Arora and P. L. Roe, A Well-Behaved TVD Limiter for High-Resolution Calculations of Unsteady Flow, Journal of Computational Physics, 132 (1997), pp. 3–11, https://doi.org/10.1006/jcph.1996.5514, http://www.sciencedirect.com/science/article/pii/S002199919695514X.
  * [2] J. P. Boris and D. L. Book, Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works, Journal of Computational Physics, 11 (1973), pp. 38–69, https://doi.org/10.1016/0021-9991(73)90147-2, http://www.sciencedirect.com/science/article/pii/0021999173901472.
  * [3] J. P. Boris and D. L. Book, Solution of continuity equations by the method of flux-corrected transport, 1976, pp. 85–129.
  * [4] A. Busnaina, X. Zheng, and M. Sharif, A modified skew upwind scheme for fluid flow and heat transfer computations, Applied Mathematical Modelling, 15 (1991), pp. 425–432, https://doi.org/10.1016/0307-904X(91)90030-S, http://www.sciencedirect.com/science/article/pii/0307904X9190030S.
  * [5] C. Carey, T. Scanlon, and S. Fraser, SUCCA-an alternative scheme to reduce the effects of multidimensional false diffusion, Applied Mathematical Modelling, 17 (1993), pp. 263–270, https://doi.org/10.1016/0307-904X(93)90048-L, http://www.sciencedirect.com/science/article/pii/0307904X9390048L.
  * [6] C. Carey, T. J. Scanlon, and S. M. Fraser, The Reduction of False Diffusion in the Simulation of Vortex Shedding, Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 206 (1992), pp. 399–406, https://doi.org/10.1243/PIME_PROC_1992_206_147_02, https://doi.org/10.1243/PIME_PROC_1992_206_147_02.
  * [7] T. J. Chung, Computational Fluid Dynamics, Cambridge University Press, 2002, https://doi.org/10.1017/CBO9780511606205.
  * [8] M. Darwish and F. Moukalled, A new approach for building bounded skew-upwind schemes, Computer Methods in Applied Mechanics and Engineering, 129 (1996), pp. 221–233, https://doi.org/10.1016/0045-7825(95)00881-0, http://www.sciencedirect.com/science/article/pii/0045782595008810.
  * [9] V. Ferreira, R. de Queiroz, G. Lima, R. Cuenca, C. Oishi, J. Azevedo, and S. McKee, A bounded upwinding scheme for computing convection-dominated transport problems, Computers & Fluids, 57 (2012), pp. 208–224, https://doi.org/10.1016/j.compfluid.2011.12.021, http://www.sciencedirect.com/science/article/pii/S0045793011003938.
  * [10] A. Harten, High resolution schemes for hyperbolic conservation laws, Journal of Computational Physics, 49 (1983), pp. 357–393, https://doi.org/10.1016/0021-9991(83)90136-5, http://www.sciencedirect.com/science/article/pii/0021999183901365.
  * [11] A. Harten, On a Class of High Resolution Total-Variation-Stable Finite-Difference Schemes, SIAM Journal on Numerical Analysis, 21 (1984), pp. 1–23, https://doi.org/10.1137/0721001, http://epubs.siam.org/doi/abs/10.1137/0721001, https://arxiv.org/abs/http://epubs.siam.org/doi/pdf/10.1137/0721001.
  * [12] Y. A. Hassan, J. G. Rice, and J. H. Kim, A STABLE MASS-FLOW-WEIGHTED TWO-DIMENSIONAL SKEW UPWIND SCHEME, Numerical Heat Transfer, 6 (1983), pp. 395–408, https://doi.org/10.1080/01495728308963096, https://doi.org/10.1080/01495728308963096.
  * [13] D. Karadimou and N. Markatos, A novel flow-oriented discretization scheme for reducing false diffusion in three-dimensional (3D) flows: An application in the indoor environment, Atmospheric Environment, 61 (2012), pp. 327–339, https://doi.org/10.1016/j.atmosenv.2012.07.033, http://www.sciencedirect.com/science/article/pii/S1352231012007054.
  * [14] A. Kriel, Error analysis of flux limiter schemes at extrema, Journal of Computational Physics, 328 (2017), pp. 371–386, https://doi.org/10.1016/j.jcp.2016.10.024, http://www.sciencedirect.com/science/article/pii/S0021999116305228.
  * [15] B. V. Leer, Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow, Journal of Computational Physics, 23 (1977), pp. 263–275, https://doi.org/10.1016/0021-9991(77)90094-8, http://www.sciencedirect.com/science/article/pii/0021999177900948.
  * [16] B. V. Leer, Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection, Journal of Computational Physics, 23 (1977), pp. 276–299, https://doi.org/10.1016/0021-9991(77)90095-X, http://www.sciencedirect.com/science/article/pii/002199917790095X.
  * [17] B. P. Leonard, Simple high-accuracy resolution program for convective modelling of discontinuities, International Journal for Numerical Methods in Fluids, 8 (1988), pp. 1291–1318, https://doi.org/10.1002/fld.1650081013, http://dx.doi.org/10.1109/84.585795.
  * [18] E. O. B. OGEDENGBE, Non-Inverted Skew Upwind Scheme for Numerical Heat Transfer and Fluid Flow Simulations, PhD thesis, 2006.
  * [19] E. O. B. Ogedengbe and G. F. Naterer, Non-inverted skew upwind scheme for three-dimensional convective transport, Numerical Heat Transfer, Part B: Fundamentals, 46 (2004), pp. 141–164, https://doi.org/10.1080/10407790490449408, https://doi.org/10.1080/10407790490449408.
  * [20] S. Osher and S. Chakravarthy, High Resolution Schemes and the Entropy Condition, SIAM Journal on Numerical Analysis, 21 (1984), pp. 955–984, https://doi.org/10.1137/0721060, http://epubs.siam.org/doi/abs/10.1137/0721060, https://arxiv.org/abs/http://epubs.siam.org/doi/pdf/10.1137/0721060.
  * [21] M. K. Patel, M. Cross, and N. C. Markatos, An assessment of flow oriented schemes for reducing ‘false diffusion’, International Journal for Numerical Methods in Engineering, 26 (1988), pp. 2279–2304, https://doi.org/10.1002/nme.1620261011, https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620261011, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620261011.
  * [22] P.L.Roe, Some contributions to the modelling of discontinuous flows, Lectures in Applied Mathematics, 22 (1985), pp. 163–193.
  * [23] G. Raithby, A critical evaluation of upstream differencing applied to problems involving fluid flow, Computer Methods in Applied Mechanics and Engineering, 9 (1976), pp. 75–103, https://doi.org/10.1016/0045-7825(76)90078-5, http://www.sciencedirect.com/science/article/pii/0045782576900785.
  * [24] G. Raithby, Skew upstream differencing schemes for problems involving fluid flow, Computer Methods in Applied Mechanics and Engineering, 9 (1976), pp. 153–164, https://doi.org/10.1016/0045-7825(76)90058-X, http://www.sciencedirect.com/science/article/pii/004578257690058X.
  * [25] G. E. Schneider and M. J. Raw, A SKEWED, POSITIVE INFLUENCE COEFFICIENT UPWINDING PROCEDURE FOR CONTROL-VOLUME-BASED FINITE-ELEMENT CONVECTION-DIFFUSION COMPUTATION, Numerical Heat Transfer, 9 (1986), pp. 1–26, https://doi.org/10.1080/10407788608913462, https://doi.org/10.1080/10407788608913462.
  * [26] R. M. Smith and A. G. Hutton, THE NUMERICAL TREATMENT OF ADVECTION: A PERFORMANCE COMPARISON OF CURRENT METHODS, Numerical Heat Transfer, 5 (1982), pp. 439–461, https://doi.org/10.1080/10407788208913458, https://doi.org/10.1080/10407788208913458, https://arxiv.org/abs/https://doi.org/10.1080/10407788208913458.
  * [27] P. K. Sweby, High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws, SIAM Journal on Numerical Analysis, vol. 21, no. 5, pp. 995-1011, 21 (1984), pp. 995–1011, https://doi.org/10.1137/0721062.
  * [28] Л. Г. Лойцянский [ L. G. Loitsiansky], Механика жидкости и газа, 7-е изд. [Fluid and gas mechanics, 7-th edition], Дрофа [Drofa], 2003.
  * [29] B. van Leer, Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme, Journal of Computational Physics, 14 (1974), pp. 361–370, https://doi.org/10.1016/0021-9991(74)90019-9, http://www.sciencedirect.com/science/article/pii/0021999174900199.
  * [30] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov’s method, Journal of Computational Physics, 32 (1979), pp. 101–136, https://doi.org/10.1016/0021-9991(79)90145-1, http://www.sciencedirect.com/science/article/pii/0021999179901451.
  * [31] B. van Leer, Upwind and high-resolution methods for compressible flow: From donor cell to residual-distribution schemes, Communications in Computational Physics, 1 (2006), p. 12.
  * [32] M. Čada and M. Torrilhon, Compact third-order limiter functions for finite volume methods, Journal of Computational Physics, 228 (2009), pp. 4118–4145, https://doi.org/10.1016/j.jcp.2009.02.020, http://www.sciencedirect.com/science/article/pii/S0021999109000953.
  * [33] F. M. White, Viscous Fluid Flow, MCGRAW HILL SERIES IN MECHANICAL ENGINEERING, McGraw-Hill Education; 3rd edition, 2005.
  * [34] D. Zhang, C. Jiang, D. Liang, and L. Cheng, A review on TVD schemes and a refined flux-limiter for steady-state calculations, Journal of Computational Physics, 302 (2015), pp. 114–154, https://doi.org/10.1016/j.jcp.2015.08.042, http://www.sciencedirect.com/science/article/pii/S0021999115005707.
  * [35] M. ZIJLEMA, ON THE CONSTRUCTION OF A THIRD-ORDER ACCURATE MONOTONE CONVECTION SCHEME WITH APPLICATION TO TURBULENT FLOWS IN GENERAL DOMAINS, International Journal for Numerical Methods in Fluids, 22 (1996), pp. 619–641, https://doi.org/10.1002/(SICI)1097-0363(19960415)22:7, https://onlinelibrary.wiley.com/doi/abs/10.1002/%28SICI%2910
97-0363%2819960415%2922%3A7%3C619%3A%3AAID-FLD370%3E3.0.CO%3B2-L.