# Logarithmic lattice models for flows with boundaries

Ciro S. Campolina  Université Côte d’Azur, CNRS, Inria, Institut de Physique
de Nice, France††thanks<EMAIL_ADDRESS>Alexei A.
Mailybaev  Instituto de Matemática Pura e Aplicada – IMPA, Rio de Janeiro,
Brazil

###### Abstract

Many fundamental problems in fluid dynamics are related to the effects of
solid boundaries. In general, they install sharp gradients and contribute to
the developement of small-scale structures, which are computationally
expensive to resolve with numerical simulations. A way to access extremely
fine scales with a reduced number of degrees of freedom is to consider the
equations on logarithmic lattices in Fourier space. Here we introduce new toy
models for flows with walls, by showing how to add boundaries to the
logarithmic lattice framework. The resulting equations retain many important
properties of the original systems, such as the conserved quantities, the
symmetries and the boundary effects. We apply this technique to many flows,
with emphasis on the inviscid limit of the Navier-Stokes equations. For this
setup, simulations reach impressively large Reynolds numbers and disclose
interesting insights about the original problem.

##  1 Introduction

The presence of solid walls installs many subtleties in the dynamics of
fluids. Boundaries naturally break homogeneity and act on the flow as sources
of vorticity. Those aspects add extra complexity to fundamental mathematical
and physical questions that are also formulated for boundaryless domains.

An important problem, for instance, is whether Navier-Stokes smooth solutions
converge to Euler’s in the vanishing viscosity limit. While the convergence is
well-established in full space [22], the case with boundaries is still not
well-understood [14]. At high Reynolds numbers, viscous effects and sharp
gradients are confined to a region called boundary layer [41], which is
supposed to remain in a vicinity of the wall, but might detach from it. The
propagation of vortices generated on the boundary to the bulk flow has been
linked to a possible anomalous dissipation in the inviscid limit. This would
prevent the convergence of Navier-Stokes solutions to Euler’s [23] and has
been object of numerical investigations [34, 33].

An archetype that has been considered in computational simulations is the case
of dipole-wall collisions [26]. Their dynamics can be visualized in the
vorticity snapshots in Fig. 1. The flow is confined between two solid
boundaries from left and right. We impose periodic boundary conditions in the
top and the bottom. We initialize the flow with a dipole of opposite
vorticities in the center of the domain—Fig. 1(a). After the initialization,
the cores of the dipole drift towards right. A sharp vorticity strip is
developed very close to the bondaries—Fig. 1(b). This is Prandtl boundary
layer, which becomes closer and closer to the boundary as viscosity decreases
[41]. When the dipole collides with the wall, small but intense vorticity
structures arise from the boundary—Fig. 1(c). They have the opposite sign of
the vortex. This sharp adverse vorticity field arising from the walls is
usually known as the boundary layer detachment. As a consequence, the large
scale flow of the dipole rebound away from the wall—Fig. 1(d)—, while small
scale structures continue to be generated at the boundaries. All this
phenomenology of boundary layers is strongly tight to the no-slip boundary
condition combined with a very small viscosity. The inviscid case, however, is
completely different. Euler’s flow cannot penetrate boundaries, but it can
slip along them, and this is what happens to inviscid solutions—consult [33,
pp. 712-713] for snapshots of this same dipole-wall collision but evolved by
Euler equations. These two contrasting states after collision—rebounded
vortices for Navier-Stokes and its slippery Euler counterpart—strongly
indicate a possible lack of convergence of the inviscid limit.

Despite the relevant numerical investigations, the mechanism of unsteady
detachment is still under debate [10, 14], while the Prandtl boundary layer
equations also manifest a rich variety of features, such as instabilities
[20], finite-time singularities [15, 27], and ill-posedness [18]. For a review
of these and other interesting mathematical problems surrounding the effects
of boundaries we refer the reader to [1, 30].

To overcome computational limitations, toy models are traditionally employed
in order to investigate physical phenomena. Good example are the shell models
[3], which have been successfully used to study important aspects of
turbulence of homogeneous (boundaryless) flows. The subtle question lies on
which features of the original model should be kept and which should be
dropped. A framework recently proposed for the formulation of toy models is
the logarithmic lattice [9]. In this configuration, we consider a reduced
number of degrees of freedom by taking the governing equations on logarithmic
lattices in Fourier space. This allows to access etremely fine scales of the
flow. Such technique was successfully applied to the study of the chaotic
blowup in the 3D incompressible Euler equations [8] and to Navier-Stokes
turbulence [9], among other applications [2, 13, 39, 38]. Nevertheless, all
the studies so far considered only homogeneous flows, that is, flows without
boundaries.

(a)

(b)

(c)

(d)

Figure 1: Direct numerical simulations of a dipole vortex colliding with a
solid wall. Snapshots are for vorticity field. The flow is confined by two
solid boundaries on left and right. Periodic boundary conditions are imposed
on top and bottom. The Navier-Stokes equations are integrated using a
pseudospectral method [45] and the solid boundary is modeled by a volume
penalization method [34]. We time-step with a classical 4th order Runge-Kutta
scheme using an integration factor. Time evolves from (a) to (d).

In this paper, we present some toy models for flows with boundaries, by
showing how to introduce boundaries on the logarithmic lattice framework. The
key challenge resides on how to model walls and restricted domains with
Fourier variables. We do it in a systematic way by extending the flow to the
whole Euclidean space through symmetrization. This process introduces some
jump singularities accross the boundaries, that need to be carefully treated.
The resulting model retains the same properties of the original equations,
such as the conserved quantities, the symmetry groups and the boundary
effects. With this model, we investigate the inviscid limit problem of the
Navier-Stokes equations. Reynolds number can be raised to impressivevily large
values and we provide some insights on this long-standing problem.

The paper is organized as follows. In Section 2 we explain how to add
boundaries on logarithmic lattices. The full 3D incompressible Navier-Stokes
equations with no-slip walls are systematically derived on this framework. In
Section 3 we consider classical shear flows, which provide great intuition and
rigor on the unusual understanding of boundaries on Fourier variables. In
Section 4 we address the inviscid limit problem of the Navier-Stokes equations
in two dimensions. In Section 5 we draw some conclusions.

##  2 Boundaries on logarithmic lattices

Logarithmic lattice fields are inherently defined on Fourier space, which is
not readily available for physical domains other than the full Euclidean
space. For this reason, it is chalenging to consider solid boundaries and
restricted domains on logarithmic lattices. Our strategy is to extend the flow
to the whole space and then to model the boundaries as surfaces immersed in
the fluid. In this process, discontinuities might incur accross the boundary,
and so we need to consider a discontinuous formulation of the governing
equations [43, 19]. This approach is inspired by the well-known immersed
boundary method [36, 37].

###  2.1 Fluid dynamics equations with immersed boundary

We consider a three-dimensional velocity field $\mathbf{u}(x,y,z)=(u,v,w)$ in
the upper-half space $y>0$ with a steady solid boundary on the plane $y=0$.
Viscous incompressible flow in this domain is governed by the classical
Navier-Stokes equations together with no-slip boundary condition

$\begin{cases}\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p+\nu\Delta\mathbf{u}\quad&\text{in}\ y>0,\\\
\nabla\cdot\mathbf{u}=0\quad&\text{in}\ y>0,\\\
\mathbf{u}=\mathbf{0}\quad&\text{on}\ y=0.\end{cases}$ (1)

System (1) also approximates the governing set of equations for smooth
boundaries of more general geometries, when considering the flow in a small
vicinity of a boundary point taken as the origin in local Cartesian
coordinates.

Our goal is to deduce a system equivalent to (1), but with the new field
variables $u$, $v$, $w$ and $p$ defined everywhere. A simple way to achieve
this is by extending the flow to the lower-half space $y<0$ through the
symmetries

$\displaystyle u(x,y,z,t)$ $\displaystyle=$ $\displaystyle u(x,-y,z,t),$ (2)
$\displaystyle v(x,y,z,t)$ $\displaystyle=-$ $\displaystyle v(x,-y,z,t),$
$\displaystyle w(x,y,z,t)$ $\displaystyle=$ $\displaystyle w(x,-y,z,t),$
$\displaystyle p(x,y,z,t)$ $\displaystyle=$ $\displaystyle p(x,-y,z,t),$

for all $x$, $y$, $z$ and $t$. Fig. 2 illustrates this symmetrization. Observe
that, because of the parity of $v$, any symmetric field (2) immediately
satisfies the no-penetration condition on the boundary

$v=0\ \text{at}\ y=0.$ (3)

However, in performing such extension, we might have introduced singularities
in the system. More precisely, the resulting reflected field variables may
present jump discontinuities in their derivatives accross the solid boundary.
Therefore, we should consider the axis $y=0$ as a steady discontinuity
surface, and the balance laws must take the jump singularities into account.
The governing equations we should consider are the discontinuous formulation
of the Navier-Stokes equations [43]

$\begin{cases}\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p+\nu\Delta\mathbf{u}-\nu\mathbf{J}(x,z,t)\delta(y)\quad&\text{in}\
\mathbb{R}^{3},\\\ \nabla\cdot\mathbf{u}=0\quad&\text{in}\ \mathbb{R}^{3},\\\
\mathbf{u}=\mathbf{0}\quad&\text{on}\ y=0,\\\ \text{+
symmetries~{}\eqref{EQ:flow_symmetry}}.\end{cases}$ (4)

where $\mathbf{J}(x,z,t)$ is the jump of the shear velocities accross the
boundary

$\mathbf{J}=\begin{pmatrix}[\partial_{y}u]\\\ 0\\\
[\partial_{y}w]\end{pmatrix}.$ (5)

Here, the square brackets $[\ \cdot\ ]$ represent the jump of a field accross
the boundary $[f]=f(y=0^{+})-f(y=0^{-})$, and $\delta(y)$ is the Dirac delta
distribution. Note that, due to the parity (2) of $v$, we have zero $y$-shear
component $[\partial_{y}v]=0$. We present details of the deduction of this
model in Appendix A.

Figure 2: Extension of (a) a flow restricted to the upper-half space $y>0$
with a no-slip boundary at $y=0$ to (b) a flow defined on the whole Euclidean
space $\mathbb{R}^{3}$ through symmetries (2) with a no-slip immersed surface
at its original place $y=0$.

The concentrated force $-\nu\mathbf{J}(x,z,t)\delta(y)$ is the shear exerted
by the solid plate on the fluid and is responsible for relaxing the flow
velocity $\mathbf{u}$ (more precisely, components $u$ and $w$) to zero on the
disontinuity surface $y=0$.

Just as pressure is obtained from incompressibility, the jump variable
$\mathbf{J}$ is determined from the no-slip condition, as we are going to see
for the logarithmic lattice formulation in the next subsection. A way to
compute it in the continuous framework is suggested in [7].

As we can see from the direct deduction of the models, the original boundary
value problem (1) and its immersed boundary formulation (4) are equivalent. A
more precise argument for such statement can be found in [7].

We remark the necessity of imposing the no-slip condition on the discontinuity
surface for the uniqueness of solutions. Indeed, the no-slip condition cannot
be lifted up from the set of equations (4), otherwise the resulting system
would be ill-posed—see Appendix B for an example of nonunique solutions from
the same initial data when no-slip condition is not explicitly prescribed.

In the inviscid case $\nu=0$, we simplify (4) to the incompressible Euler
equations

$\begin{cases}\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p\quad&\text{in}\ \mathbb{R}^{3},\\\ \nabla\cdot\mathbf{u}=0\quad&\text{in}\
\mathbb{R}^{3},\\\ \text{+ symmetries~{}\eqref{EQ:flow_symmetry}}.\end{cases}$
(6)

In this system, no penetration

$v=0\quad\text{in}\ y=0$ (7)

is satisfied as a consequence of the flow symmetry, while the flow may slip
through the boundary, i.e. $u\neq 0$ and $w\neq 0$ on $y=0$.

###  2.2 Logarithmic lattice formulation

In the immmersed boundary formulation, the flow is defined in the whole
Euclidean space. This allows us to readily work with Fourier variables, and so
on logarithmic lattices.

We shall consider a three-dimensional logarithmic lattice
$\mathbb{\Lambda}^{3}$, where

$\mathbb{\Lambda}=\\{\pm 1,\pm\lambda,\pm\lambda^{2},\dots\\},$ (8)

for some $\lambda>1$. The usual choices are $\lambda=2$,
$\lambda=(1+\sqrt{5})/2\approx 1.618$ (the golden ratio), and
$\lambda=(\sqrt[3]{9+\sqrt{69}}+\sqrt[3]{9-\sqrt{69}})/\sqrt[3]{18}\approx
1.325$ (the plastic number), which result in nontrivial nonlinear interactions
among modes [9]. This lattice mimics Fourier space with largest integral scale
$L\sim 2\pi$ corresponding to $|\mathbf{k}|\sim 1$.

We represent the velocity field
$\mathbf{u}(\mathbf{k},t)=(u,v,w)\in\mathbb{C}^{3}$ as a function of the wave
vector $\mathbf{k}=(k_{x},k_{y},k_{z})\in\mathbb{\Lambda}^{3}$ on the
logarithmic lattice and the time variable $t\in\mathbb{R}$. Similarly, we have
the scalar pressure $p(\mathbf{k},t)\in\mathbb{C}$. These and all field
variables are supposed to satisfy the reality condition
$\mathbf{u}(-\mathbf{k},t)=\overline{\mathbf{u}(\mathbf{k},t)}$, where the
overbar indicates complex conjugation.

In physical space, the flow is reflected with respect to the plane $y=0$ as in
(2). Accordingly, we demand the lattice fields to satisfy the corresponding
symmetries

$\displaystyle u(k_{x},k_{y},k_{z},t)$ $\displaystyle=$ $\displaystyle
u(k_{x},-k_{y},k_{z},t),$ (9) $\displaystyle v(k_{x},k_{y},k_{z},t)$
$\displaystyle=-$ $\displaystyle v(k_{x},-k_{y},k_{z},t),$ $\displaystyle
w(k_{x},k_{y},k_{z},t)$ $\displaystyle=$ $\displaystyle
w(k_{x},-k_{y},k_{z},t),$ $\displaystyle p(k_{x},k_{y},k_{z},t)$
$\displaystyle=$ $\displaystyle p(k_{x},-k_{y},k_{z},t),$

for all $k_{x}$, $k_{y}$, $k_{z}$ and $t$.

We have also to consider the action of the shear force at the boundary
$\mathbf{F}=(F_{x},F_{y},F_{z})$ on the flow, which is related to the jump
discontinuities $\mathbf{J}=(J_{x},J_{y},J_{z})$ in the form

$\mathbf{F}(\mathbf{k},t)=-\nu\mathbf{J}(k_{x},k_{z},t)\delta(k_{y}).$ (10)

Since the jump discontinuities in physical space occur at the plane $y=0$, the
corresponding lattice variable is a function independent of $k_{y}$, and thus
depends on $k_{x}$, $k_{z}$ and $t$ only. In analogy with the Fourier
transform of Dirac delta distribution, we take the lattice Dirac delta
function $\delta(k_{y})$ as unity

$\delta(k_{y})=1\quad\text{for all}\ k_{y}\in\mathbb{\Lambda}.$ (11)

Such natural definition preserves some important properties of Dirac delta,
like scaling invariance

$\delta(\lambda k_{y})=\delta(k_{y})$ (12)

and parity

$\delta(-k_{y})=\delta(k_{y}).$ (13)

Moreover, Dirac delta (11) on the lattice keeps similarity with a classical
property of the original distribution: if $f(k_{y})$ is a lattice function
representing a function $F(y)$ in physical space, then its $k_{y}$-inner
product against delta mimics the localization of $F$ on $y=0$, since

$(f,\delta)_{k_{y}}=\sum_{k_{y}\in\mathbb{\Lambda}}f(k_{y})\overline{\delta(k_{y})}=\sum_{k_{y}\in\mathbb{\Lambda}}f(k_{y})\simeq\int\hat{F}(k_{y})dk_{y}=F(y)\Big{|}_{y=0}.$
(14)

Here the left-hand side of $\simeq$ corresponds to a logarithmic lattice
representation $k_{y}\in\mathbb{\Lambda}$, and the right-hand side corresponds
to a usual representation in continuous space $y\in\mathbb{R}$. This
interpretation allows us to impose the no-slip boundary condition in the form

$(\mathbf{u},\delta)_{k_{y}}=\mathbf{0}\quad\text{for all}\
k_{x},k_{z}\in\mathbb{\Lambda},t\in\mathbb{R},$ (15)

since we interpret

$(\mathbf{u},\delta)_{k_{y}}\simeq\mathbf{u}\Big{|}_{y=0}.$ (16)

Observe that the left expression in Eq. (15) is a function of $k_{x}$, $k_{z}$
and $t$.

With all the above definitions, we can now establish the incompressible
Navier-Stokes equations with a solid boundary on a logarithmic lattice. We
simply write the immersed boundary formulation of the Navier-Stokes equations
(4), but considering the fields and operations on logarithmic lattices

$\begin{cases}\partial_{t}\mathbf{u}+\mathbf{u}\ast\nabla\mathbf{u}=-\nabla
p+\nu\Delta\mathbf{u}-\nu\mathbf{J}(k_{x},k_{z},t)\delta(k_{y})\quad&\text{in}\
\mathbb{\Lambda}^{3},\\\ \nabla\cdot\mathbf{u}=0\quad&\text{in}\
\mathbb{\Lambda}^{3},\\\ (\mathbf{u},\delta)_{k_{y}}=\mathbf{0}\quad&\text{for
all}\ k_{x},k_{z},t,\\\ \text{+
symmetries~{}\eqref{EQ:flow_symmetry_lattice}}.\end{cases}$ (17)

The corresponding incompressible Euler equations with solid boundary on a
logarithmic lattice reads

$\begin{cases}\partial_{t}\mathbf{u}+\mathbf{u}\ast\nabla\mathbf{u}=-\nabla
p\quad&\text{in}\ \mathbb{\Lambda}^{3},\\\
\nabla\cdot\mathbf{u}=0\quad&\text{in}\ \mathbb{\Lambda}^{3},\\\ \text{+
symmetries~{}\eqref{EQ:flow_symmetry_lattice}}.\end{cases}$ (18)

In this case, the no-penetration boundary condition $(v,\delta)_{k_{y}}=0$ is
a consequence of the symmetry $v(-k_{y})=-v(-k_{y})$ on $v$ from (9). We
recall that the derivatives on logarithmic lattices are given by the Fourier
factors

$\partial_{j}f(\mathbf{k})=ik_{j}f(\mathbf{k}),$ (19)

and that the star $\ast$ product is the convolution on the logarithmic lattice

$(f\ast g)(\mathbf{k})=\sum_{\mathbf{p},\mathbf{q}\in\mathbb{\Lambda}^{3},\\\
\mathbf{p}+\mathbf{q}=\mathbf{k}}f(\mathbf{p})g(\mathbf{q}).$ (20)

We refer the reader to [9] for a comprehensive description of the logarithmic
lattice framework.

Taking the divergence of the momentum equation in (17) and invoking
incompressibility, we obtain a Poisson equation for the pressure

$\Delta
p=-\nabla\cdot(\mathbf{u}\ast\nabla\mathbf{u})-\nu\nabla\cdot(\mathbf{J}(k_{x},k_{z},t)\delta(k_{y})).$
(21)

This allows us to eliminate the pressure from system (17). Observe, however,
that the pressure is written as a function of not only the velocities, but of
the jumps as well. The additional term
$-\nu\nabla\cdot(\mathbf{J}(k_{x},k_{z},t)\delta(k_{y}))$ in pressure’s
equation stands for the contribution of the boundary.

###  2.3 Computation of the jump lattice variable

Just as pressure is obtained from the incompressibility constraint, the jump
variables are computed from the no-slip condition. We show here how such
calculation can be evaluated. The idea relies on first truncating the lattice
to later pass the limit on the truncation. The truncated system is a finite
dimensional system, whose jump contribution can be readily computed.

Let us truncate the lattice $\mathbb{\Lambda}$ by consering only frequencies
smaller than the wave number $k_{N}=\lambda^{N-1}$ as

$\mathbb{\Lambda}_{N}=\\{\pm
1,\pm\lambda,\pm\lambda^{2},\dots,\pm\lambda^{N-1}\\}.$ (22)

Then, the Dirac delta on the lattice (11) is automatically “regularized”,
since truncation turns it into a summable function.

Then we consider the approximated system

$\begin{cases}\partial_{t}\mathbf{u}^{N}+\mathbf{u}^{N}\ast\nabla\mathbf{u}^{N}=-\nabla
p^{N}+\nu\Delta\mathbf{u}^{N}-\nu\mathbf{J}^{N}(k_{x},k_{z},t)\delta^{N}(k_{y})\quad&\text{in}\
\mathbb{\Lambda}^{3}_{N},\\\ \nabla\cdot\mathbf{u}^{N}=0\quad&\text{in}\
\mathbb{\Lambda}^{3}_{N},\\\
(\mathbf{u}^{N},\delta^{N})_{k_{y}}=\mathbf{0}\quad&\text{for all}\
k_{x},k_{z},t,\\\ \text{+
symmetries~{}\eqref{EQ:flow_symmetry_lattice}}.\end{cases}$ (23)

We wrote the superscript $N$ to indicate all truncated variables, i.e. those
defined on the truncated lattice (22).

To compute the truncated jumps $\mathbf{J}^{N}$, we take the $k_{y}$-inner
product of the momentum equation with $\delta^{N}$. As a consequence of no-
slip condition, the contribution of time variation vanishes, and we are lead
to

$(\mathbf{u}^{N}\ast\nabla\mathbf{u}^{N},\delta^{N})_{k_{y}}=(-\nabla
p^{N}+\nu\Delta\mathbf{u}^{N},\delta^{N})_{k_{y}}+(-\nu\mathbf{J}^{N}(k_{x},k_{z},t)\delta^{N}(k_{y}),\delta^{N})_{k_{y}}.$
(24)

Using the fact that the jumps $\mathbf{J}^{N}$ do not depend on $k_{y}$, one
may write

$(-\nu\mathbf{J}^{N}(k_{x},k_{z},t)\delta^{N}(k_{y}),\delta^{N})_{k_{y}}=-\nu\mathbf{J}^{N}(k_{x},k_{z},t)(\delta^{N},\delta^{N})_{k_{y}}.$
(25)

Because the approximated truncated $\delta^{N}$ is summable, the product
$(\delta^{N},\delta^{N})_{k_{y}}$ appearing in Eq. (25) is a well-defined
positive number. Substitution of expression (25) into Eq. (24) yields, after
some manipulations,

$\mathbf{J}^{N}(k_{x},k_{z},t)=\frac{1}{\nu(\delta^{N},\delta^{N})_{k_{y}}}(-\mathbf{u}^{N}\ast\nabla\mathbf{u}^{N}-\nabla
p^{N}+\nu\Delta\mathbf{u}^{N},\delta^{N})_{k_{y}}.$ (26)

Because of the flow symmetries (9), we have $J_{y}^{N}=0$, as expected. We
remark, however, that this is not a closed formula for the regularized jumps
$\mathbf{J}^{N}$, but an implicit equation. Indeed, the pressure must be
solved from the Poisson equation (21) in terms of both the velocities
$\mathbf{u}^{N}$ and the jumps $\mathbf{J}^{N}$. The resulting equations can
be solved explicitly for $\mathbf{J}^{N}$, but we omit here the laborious
computations. In Sections 3 and 4 we present the resulting explict formula for
the simpler cases of one and two dimensions, respectively.

To complete the computations, we just need to take the limit $N\to\infty$.

###  2.4 Basic symmetries and balance laws

The Euler equations on logarithmic lattices are known for having the same
symmetries as the original equations [4]. When adding the boundary, the new
restricted symmetry group should preserve the imposed flow symmetries (9). To
enumerate them, if $\mathbf{u}(\mathbf{k},t)$ is a solution of system (18),
then the following transformations also yield solutions:

  1. (E.1)

(Time translations) $\\!\begin{aligned}
\mathbf{u}^{\tau}(\mathbf{k},t)=\mathbf{u}(\mathbf{k},t+\tau),\end{aligned}$
for any $\tau\in\mathbb{R}$;

  2. (E.2)

(Space translations in $x$ and $z$) $\\!\begin{aligned}
\mathbf{u}^{\boldsymbol{\xi}}(\mathbf{k},t)=e^{-i\mathbf{k}\cdot\boldsymbol{\xi}}\mathbf{u}(\mathbf{k},t),\end{aligned}$
for any $\boldsymbol{\xi}=(\xi_{x},0,\xi_{z})\in\mathbb{R}^{3}$;

  3. (E.3)

(Isotropy in $x$ and $z$ and parity)
$\mathbf{u}^{R}(\mathbf{k},t)=R^{-1}\mathbf{u}(R\mathbf{k},t),$ where $R$ is
any transformation $(k_{1},k_{2},k_{3})\mapsto(\pm k_{\alpha},\pm k_{2},\pm
k_{\beta})$ with $(\alpha,\beta)$ permutations of $(1,3)$;

  4. (E.4)

(Scale invariance) $\\!\begin{aligned}
\mathbf{u}^{n,h}(\mathbf{k},t)=\lambda^{h}\mathbf{u}\left(\lambda^{n}\mathbf{k},\lambda^{h-n}t\right),\end{aligned}$
for any $h\in\mathbb{R}$ and $n\in\mathbb{Z}$, where $\lambda$ is the lattice
spacing;

  5. (E.5)

(Time reversibility) $\\!\begin{aligned}
\mathbf{u}^{r}(\mathbf{k},t)=-\mathbf{u}\left(\mathbf{k},-t\right)\end{aligned}$;

  6. (E.6)

(Galilean invariance in $x$ and $z$) $\\!\begin{aligned}
\mathbf{u}^{\mathbf{v}}(\mathbf{k},t)=e^{-i\mathbf{k}\cdot\mathbf{v}t}\mathbf{u}(\mathbf{k},t)-\widehat{\mathbf{v}}(\mathbf{k}),\end{aligned}$
for any $\mathbf{v}=(v_{x},0,v_{z})\in\mathbb{R}^{3}$, where
$\widehat{\mathbf{v}}(\mathbf{k})$ is the constant velocity field on the
lattice defined as $\widehat{\mathbf{v}}(\mathbf{0})=\mathbf{v}$ and zero for
$\mathbf{k}\neq\mathbf{0}$.

Naturally, Galilean invariance (E.6) is well-defined only if we add the zero
component to the lattice $0\in\mathbb{\Lambda}$.

Since the equations are not modified by introducing a boundary, the conserved
quantities are also the same, say: the energy

$E(t)=\frac{1}{2}\sum_{\mathbf{k}\in\mathbb{\Lambda}^{d}}|\mathbf{u}(\mathbf{k},t)|^{2}$
(27)

and helicity

$H(t)=\sum_{\mathbf{k}\in\mathbb{\Lambda}^{d}}\mathbf{u}(\mathbf{k},t)\overline{\boldsymbol{\omega}(\mathbf{k},t)}$
(28)

in the three-dimensional case $d=3$, and the energy (27) and enstrophy

$\Omega(t)=\frac{1}{2}\sum_{\mathbf{k}\in\mathbb{\Lambda}^{d}}|\boldsymbol{\omega}(\mathbf{k},t)|^{2}$
(29)

in the two-dimensional case $d=2$. Here,
$\boldsymbol{\omega}=\nabla\times\mathbf{u}$ are the vorticities. Kelvin’s
Theorem for the conservation of circulation also holds in the presence of a
boundary. The precise statements and the proofs for those conservations laws
can be found in [4, 9].

In the case of positive viscosity, the Navier-Stokes equations with and
without boundary share the same scaling symmetry

  1. (NS)

(Scale invariance) $\\!\begin{aligned}
\mathbf{u}^{n}(\mathbf{k},t)=\lambda^{-n}\mathbf{u}\left(\lambda^{n}\mathbf{k},\lambda^{-2n}t\right),\end{aligned}$
for any $n\in\mathbb{Z}$, where $\lambda$ is the lattice spacing.

Thus, the introduction of a boundary through our modelling technique does not
disrupt the self-similarity properties of the Navier-Stokes equations.

Additionally, the shear force
$\mathbf{F}(k_{x},k_{y},k_{z},t)=-\nu\mathbf{J}(k_{x},k_{z},t)\delta(k_{y})$
on the boundary exerts no work in the flow, as a consequence of no-slip
boundary condition

$(\mathbf{F},\mathbf{u})=-\nu\sum_{k_{x},k_{z}\in\mathbb{\Lambda}}\mathbf{J}(k_{x},k_{z},t)\cdot(\mathbf{u},\delta)_{k_{y}}=0.$
(30)

This proportionates the classical energy balance law

$\frac{dE}{dt}=-2\nu\Omega(t),$ (31)

where $\Omega$ is the enstrophy (29).

The conservation laws and the energy balance are also satisfied by the
truncated flows $\mathbf{u}^{N}$.

##  3 Classical shear flows

In this section, we take one step back by considering some classical shear
flows. Due to their simplicity, we can study the logarithmic lattice solutions
in the light of exact expected results, or even compare them with direct
numerical simulations. Naturally, the comparisons are always in terms of
qualitative behavior.

Consider the governing system of equations (1) for a three-dimensional flow on
the upper-half space $y>0$ with a solid boundary on the plane $y=0$. Let us
assume such flow has no variation with respect to $x$ and $z$, and that
$v=w=0$. Under those hypotheses, incompressibility is trivially satisfied,
while pressure is constant. Then, the resulting flow simplifies to a one-
dimensional velocity field $u=u(y,t)$ governed by

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u+f\quad&\text{in}\ y>0,\\\
u=0\quad&\text{on}\ y=0,\end{cases}$ (32)

where $f=f(y,t)$ is a possible external force. System (32) is supplemented by
proper initial conditions

$u\Big{|}_{t=0}=u^{0}(y).$ (33)

This is the Dirichlet problem for the one-dimensional linear heat equation.
Such system can be solved exactly using the heat kernel and reflections—see
Appendix B for a closed formula. Here, we understand its solutions as simple
shear flows over a solid plate at $y=0$.

Following the steps from Section 2, we can formulate problem (32) on the whole
domain, considering a discontinuity point at the origin $y=0$. The resulting
immersed boundary formulation reads

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u+f-\nu
J(t)\delta(y)\quad&\text{in}\ \mathbb{R},\\\ u(y)=u(-y)\quad&\text{in}\
\mathbb{R},\\\ u=0\quad&\text{on}\ y=0.\end{cases}$ (34)

The jump discontinuity is derived from system (34) as

$J(t)=[\partial_{y}u].$ (35)

On a logarithmic lattice

$\mathbb{\Lambda}=\\{0,\pm 1,\pm\lambda,\pm\lambda^{2},\dots\\},$ (36)

model (34) reads

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u+f-\nu
J(t)\delta(k)\quad&\text{in}\ \mathbb{\Lambda},\\\ u(k)=u(-k)\quad&\text{in}\
\mathbb{\Lambda},\\\ (u,\delta)=0\quad&\text{for all}\
t\in\mathbb{R}.\end{cases}$ (37)

As we saw earlier, the jump $J(t)$ can be computed explicitly. If the lattice
(36) is truncated up to the wave number $k_{N}=\lambda^{N-1}$ as

$\mathbb{\Lambda}_{N}=\\{0,\pm
1,\pm\lambda,\pm\lambda^{2},\dots\lambda^{N-1}\\},$ (38)

we have the truncated model

$\begin{cases}\partial_{t}u^{N}=\nu\partial_{y}^{2}u^{N}+f^{N}-\nu
J^{N}(t)\delta^{N}(k)\quad&\text{in}\ \mathbb{\Lambda}_{N},\\\
u^{N}(k)=u^{N}(-k)\quad&\text{in}\ \mathbb{\Lambda}_{N},\\\
(u^{N},\delta^{N})=0\quad&\text{for all}\ t\in\mathbb{R}.\end{cases}$ (39)

Taking the inner product of the main equation in (39) against $\delta^{N}$,
and using the no-slip condition $(u^{N},\delta^{N})=0$, we obtain

$(\nu\partial_{y}^{2}u^{N}+f^{N},\delta^{N})+(-\nu
J^{N}(t)\delta^{N},\delta^{N})=0.$ (40)

Next, using the fact that $J^{N}(t)$ does not depend on $k$, we can evaluate
the term

$(-\nu J^{N}(t)\delta^{N},\delta^{N})=-\nu J^{N}(t)(\delta^{N},\delta^{N}).$
(41)

Observe that, since the lattice is now truncated, the expression
$(\delta^{N},\delta^{N})$ is the well-defined positive number

$(\delta^{N},\delta^{N})=\sum_{k\in\mathbb{\Lambda}^{N}}1=2N+1.$ (42)

Finally, we substitute (41) into (40) and isolate $J^{N}(t)$ to obtain a
closed formula for the approximated jump

$J^{N}(t)=\frac{(\partial_{y}^{2}u^{N}+f^{N},\delta^{N})}{(\delta^{N},\delta^{N})}.$
(43)

The original jump is then recovered from the sequence of its approximations

$J(t)=\lim\limits_{N\to\infty}J^{N}(t)=\lim\limits_{N\to\infty}\frac{(\partial_{y}^{2}u^{N}+f^{N},\delta^{N})}{(\delta^{N},\delta^{N})}.$
(44)

###  3.1 Couette flow

Let us consider a flow between two parallel plates separated by a unit
distance. One of the plates is at rest, and the other moves with a constant
horizontal speed $V$. This classical problem has well-known stationary
solution given by the linear velocity profile

$u(y)=Vy\quad\text{for}\ 0\leq y\leq 1.$ (45)

Such configuration is called Couette flow—consult e.g. [28, §17]. Then, the
shear force $f_{1}$ exerted on the fluid by the moving plate at $y=1$ is

$f_{1}=\nu V\quad\text{at}\ y=1,$ (46)

while the force $f_{0}$ from the plate $y=0$ at rest is the symmetric
counterpart

$f_{0}=-\nu V\quad\text{at}\ y=0,$ (47)

To model a similar phenomenon on a logarithmic lattice, we take the one-
dimensional shear flow equation (37) and consider the action of the moving
plate as a constant-in-time force $f$ applied at $k=0$ in the form

$f(k,t)=2\nu V\delta_{0}(k),$ (48)

with

$\delta_{0}(k)=\begin{cases}1,\quad\text{if}\ k=0,\\\
0,\quad\text{otherwise}.\end{cases}$ (49)

We concentrate the force at $k=0$ in order to model the momentum input due to
the moving boundary. The choice of the force (48) is motivated by the known
shear action (46) at the moving plate from the original problem. Such force is
proportional to the relative velocitity between the plates and to the fluid
viscosity. The factor $2$ appears in (48) because our reflected flow doubles
the forces on the discontinuity surface.

####  3.1.1 Stationary solution

Let us look for stationary solutions of system (37) under the action of the
force (48). Evaluating the governing equation at $k=0$ gives us the value of
the jump

$J=2V.$ (50)

The solution at $k\neq 0$ can now be obtained, using that $f(k)=0$ at $k\neq
0$, as

$0=\nu\partial^{2}_{y}u(k)-J\delta(k)=-\nu k^{2}u(k)-\nu J,$ (51)

which yields

$u(k)=-2Vk^{-2}\quad\text{for}\ k\neq 0.$ (52)

Finally, the mean flow $u(k=0)$ can be computed from the no-slip boundary
condition, as follows

$0=(u,\delta)=\sum_{k\in\mathbb{\Lambda}}u(k)=u(0)+\sum_{k\in\mathbb{\Lambda}\backslash\\{0\\}}u(k)=u(0)-2V\sum_{k\in\mathbb{\Lambda}\backslash\\{0\\}}k^{-2},$
(53)

whence

$u(0)=2V\sum_{k\in\mathbb{\Lambda}\backslash\\{0\\}}k^{-2}.$ (54)

The final solution

$J=2V\quad\text{and}\quad
u(k)=\begin{cases}2V\sum_{\tilde{k}\in\mathbb{\Lambda}\backslash\\{0\\}}\tilde{k}^{-2}\quad&\text{for}\
k=0,\\\ -2Vk^{-2}\quad&\text{for}\ k\neq 0,\end{cases}$ (55)

shares many similitudes with the original shear flow. First, the jump agrees
exactly with the original velocity profile (45), if we consider its symmetric
reflection $u(-y)=u(y)$ around the origin $y=0$. Consequently, the shear force
at the plate is $-\nu J=-2\nu V$, which is in agreement with the force (47).
Remember that the reflections of the flow around the origin double the forces,
which explains the extra factor $2$. Second, we observe a solution tail
proportional to $k^{-2}$, which is the expected Fourier spectrum for a
function whose first derivative is discontinuous [17].

####  3.1.2 Unsteady solutions

We can also simulate unsteady solutions of system (37) under the action of the
force (48). For this, we set $V=1$, $\nu=1$ and consider identically zero
initial conditions $u^{0}\equiv 0$. We take the truncated logarithmic lattice
(38) with $\lambda=2$ and $N=50$. Together, they provide the finest scale
$\ell_{N}=1/k_{N}\approx 10^{-15}$. We integrate the equations using Matlab’s
ode15s solver [42], with the tolerances $RelTol=10^{-8}$ and
$AbsTol=10^{-11}$.

Fig. 3 shows the time evolution of our Couette-like flow. As time advances,
the solution converges to the stationary state (55), which is a fixed point
attractor of the system. At each instant, the solution presents two different
regimes along scales: a constant plateau $|u|\approx const.$ state at larger
scales, followed by a $|u|\propto k^{-2}$ tail at smaller scales. Such
behavior is more noticeable in earlier instants, since the plateau range
shrinks as time advances. These two regimes can be understood by establishing
the corresponding asymptotic solutions, as we do now. For this analysis, we
consider that the lattice $\mathbb{\Lambda}$ has infinitely large and
infinitely small $k$.

Figure 3: Dynamical evolution of the Couette flow on a logarithmic lattice.
(a) Spectra of solution $|u(k)|$ as a function of the wave number $k$ at
several instants. Time advances from blue to red in logarithmic scale, from
$t=10^{-12}$ up to $t=10^{2}$ by factors of $10$. The solution converges to
the stationary power-law $k^{-2}$, represented by the black dashed line. (b)
Time evolution of the jump variable $J(t)$, converging to the constant steady
value $J=2V$, with $V=1$, as time advances.

The asymptotic solutions are obtained by analysing how the terms in the
governing equation scale with respect to $k$. First, we notice that the
boundary force $-\nu J(t)\delta(k)$ acts with the same magnitude at all
scales, since $\delta\equiv 1$. Therefore, it contributes both to large and to
small scales. Next, we analyze the contribution of the dissipative term
$\nu\partial_{y}^{2}u(k,t)=-\nu k^{2}u(k,t)$, which determines the two
different regimes. At high $k$, the dissipative term becomes important and
dominates over the time derivative $\partial_{t}u$. In this regime, we have
the balance

$\nu\partial_{y}^{2}u\sim\nu J\delta\quad\text{as}\ k\to\infty,$ (56)

which results in the tail asymptotic solution

$u(k,t)\sim-J(t)k^{-2}\quad\text{as}\ k\to\infty.$ (57)

Figure 4: Time evolution of $u(k,t)$ at small $k\in\mathbb{\Lambda}$ on the
logarithmic lattice. Axes are in logarithmic scales. We plot in green the
predicted asymptotic value $v(t)$. Black curves show actual computed solutions
$u(k,t)$. Different styles of lines stand for different scales $k$.

On the other hand, if $k$ is small, the dissipative term becomes negligible,
and it is dominated by the time derivative $\partial_{t}u$. The balance in
this case is

$\partial_{t}u\sim-\nu J\delta\quad\text{as}\ k\to 0.$ (58)

Eq. (58) is independent from $k$, which justifies the constant plateau
behavior

$u(k,t)\approx v(t)\quad\text{for small}\ k,$ (59)

in Fig. 3(a). To confirm such statement, let us consider the ordinary
differential equation

$\frac{dv}{dt}=-\nu J(t),$ (60)

which is the exact form for the asymptotic balances (58) and (59). We
integrate this equation numerically using the numerical result of $J(t)$ from
our model. Fig. 4 compares this predicted value for the plateau with the
actual solution $u(k,t)$ for small $k$. We verify that as $k$ decreases, the
agreement between the prediction and the solution improves.

We observe the transition between the two regimes when $|u(k,t)|$ is of order
of $k^{2}$, which eventually occurs at any finite scale $k$ after sufficiently
large time. We can see the transition of regimes in Fig. 3(a) and the expected
deviation from the asymptotic solution at finite $k$ after some time in Fig.
4.

###  3.2 Decaying shear flow

Next we consider a decaying shear flow. The setup consists of system (32) with
zero external force $f\equiv 0$ and non-zero initial condition
$u^{0}\not\equiv 0$.

First, we shall solve the continuous-space immersed boundary formulation (34)
with direct numerical simulations. For this, we employ simple finite
difference schemes with a regularized Dirac delta parametrized by the number
of grid points. Such strategy will illustrate the explicit computation of jump
singularities. Then, we shall consider the logarithmic lattice formulation.
The results on logarithmic lattices will be compared with the precedent direct
numerical simulations.

####  3.2.1 Direct numerical simulations in continuous space

Consider the one-dimensional shear flow equation with immersed boundary (34)
in continuous space $y\in\mathbb{R}$, with $f\equiv 0$. For the numerical
model, we shall restrict the domain to the bounded interval $y\in[-L,L]$ and
simulate the dynamics on a finite interval of time $t\in[0,T]$. Such framework
is a good approximation of the unbounded system $y\in\mathbb{R}$ when the
initial condition has fast decrease at infinity and $T$ is sufficiently small.
In what follows, we fix $L=10$.

We discretize both space and time and consider second-order centered finite
differences for the laplacian and first-order forward finite difference for
the time derivative [29]. The jump can be explicitly computed with the
discrete variables. This yields a recurrence relation that can be used to have
numerical approximations of the dynamics.

Figure 5: Direct numerical simulations of the continuous-space immersed
boundary formulation of a decaying shear flow. (a) Solution $u$ along space
$y$. Different colors indicate solutions at different instants $t$. Time
advances from blue to red. (b) Fourier spectra of solution in logarithmic
scales. Time advances from from blue to red in logarithmic scale. (c) Time
evolution of the jump $J(t)$.

For the numerical experiments, we set $\nu=1$, and the initial condition to

$u^{0}(y)=\left[1-\cos\left(\frac{\pi
y}{5}\right)\right]e^{-y^{2}/2}\quad\text{for}\ y\in\mathbb{R}.$ (61)

This function satisfies the odd symmetry $u(-y)=u(y)$ and decays exponentially
fast at infinity $|y|\to\infty$. Moreover, the initial condition satisfies the
no-slip condition $u^{0}(0)=0$ and has zero initial shear force as
$[\partial_{y}u^{0}]=0$.

Our first numerical integrations concern the general picture of the solution
in physical space after some amount of time. We use $T=2$, $N=1024$ points in
space and $M=20~{}000$ points in time. The results are depicted in Figs. 5(a)
and 5(c). To reach better resolution in Fourier space at earlier instants, we
change the numerical parameters to $T=0.01$ and $N=4096$, while $M$ is kept
the same. The Fourier spectrum is drawn in Fig. 5(b). We now discuss the
results.

We show the time evolution of the velocity along physical space in Fig. 5(a).
The initial data (61) is essentially two symmetric bumps close to the origin
quickly decaying at infinity. They represent an initial motion restricted to
large scales—as we verify from the spectrum concentrated on small $k$ in Fig.
5(b). As time advances, the solution dissipates and tends towards the
identically zero steady state. Initially, the velocity profile develops a
large increase in the jump, it achieves a maximum and then decays, as in Fig.
5(c).

In Fourier space, we see some features that we already presented for the
logarithmic model of a Couette-like flow. The initial state is confined to
large scales, but the boundary shear force affects all scales at every instant
of time $t>0$. This is noticed by the instantaneous development of a constant
plateau regime along intermediate scales. Such plateau would be followed by a
$k^{-2}$ tail at large scales even for very small $t$, if we simulated for
larger $k$. Under limited resolution, we can see the power-law development at
later times in Fig. 5(b), when the tail already reaches intermediate and large
scales. We observe some slight deviation from the $k^{-2}$ at the higher
scales, as an effect of truncation. Increasing resolution extends the power-
law and sends this truncation effect towards higher $k$.

####  3.2.2 Logarithmic lattice simulations

Figure 6: Dynamical evolution of a decaying shear flow on a logarithmic
lattice. (a) Spectra of solution $|u(k)|$ as a function of the wave number $k$
at several instants. Time advances from blue to red in logarithmic scale, from
$t=10^{-10}$ up to $t=10^{0}$ by factors of $10$. The dashed lines show the
corresponding asymptotic solutions $u(k,t)\sim-J(t)k^{-2}$ for large $k$. (b)
Time evolution of the jump variable $J(t)$.

Now, let us consider the decaying shear flow on a logarithmic lattice. We take
the governing system of equations (37) with unit viscosity $\nu=1$ and zero
external force $f\equiv 0$. Representing the initial condition $u(k,0)$ by the
vector $u^{0}=(u^{0}(0),u^{0}(\pm
1),u^{0}(\pm\lambda),u^{0}(\pm\lambda^{2}),\dots)$, we fix

$u^{0}=\left(8,\frac{14}{3},-\frac{67}{6},\frac{5}{2},0,0,\dots\right).$ (62)

Such initial flow satisfies the no-slip condition $(u^{0},\delta)=0$ and has
no initial shear flow at the boundary, since
$(\partial_{y}^{2}u^{0},\delta)=0$, and thus $J|_{t=0}=0$.

The numerical simulations are undertaken on the truncated lattice (38) with
$\lambda=2$ and $N=50$. We solve the equations with Matlab’s ode15s solver
[42], with tolerances $RelTol=10^{-8}$ and $AbsTol=10^{-11}$.

Figure 7: Convergence of solutions with respect to truncation. Error in
$\ell^{2}$ norm in time for the jump $J^{N}$ with respect to a reference
solution with $75$ as a function of the number of nodes $N$. Different curves
stand for the three lattice spacings: $\lambda=2$ in blue, $\lambda=\varphi$
(the golden mean) in red, and $\lambda=\sigma$ (the plastic number) in green.
Dashed lines give reference to exponential convergence $\propto\lambda^{-2N}$.

In Fig. 6 we see the dynamics of the solution on the logarithmic lattice. Fig.
6(a) shows the spectra of the velocites for different instants of time. The
initial spectrum is confined to large scales, but the flow instantly develops
a $k^{-2}$ tail at high $k$. The gap between large and small scales is filled
by a constant plateau state, which shrinks as time advances. In dashed lines,
we plot the asymptotic solutions (57) and verify that they match the computed
results at high $k$.

Here we verify that the solution on the logarithmic lattice is quite similar
to that obtained from DNS, with the advantage that in the former we can reach
way finer scales. The observations of a $k^{-2}$ tail, a plateau in
intermediate scales and the dynamical shrink of the gap between them is also
present in Fig. 5(b) for the continuous space model, but in a more restricted
resolution. Moreover, the time evolution of the jump lattice variable shown in
Fig. 6(b) is qualitatively indistinguishable from the corresponding jump on
the continuous model from DNS in Fig. 5(c).

###  3.3 Convergence with respect to truncation

Heuristically, truncation of the logarithmic lattice $\mathbb{\Lambda}_{N}$ at
a sufficiently large wave number $k_{N}=\lambda^{N-1}$ represents a cutoff of
the $k^{-2}$ tail whose sum is proportional to $k_{N}^{-2}$. This establishes
the exponential convergence rate $\lambda^{-2N}$.

To verify such convergence rate with respect to truncation, we perform
simulations of the decaying shear flow from Section 3.2.2 for several number
of nodes. We range $N$ from $5$ up to $70$ by unit increments. To compute
errors, we consider the simulation $N=75$ as reference. We integrate equations
using Matlab’s ode15s solver with tolerances $RelTol=10^{-13}$ and
$AbsTol=10^{-16}$. We fix the time window $T=3$ and sample the solution at
points $t_{j}=j\Delta t$, for $j=0,1,\dots,M$ with $M=3000$, so $\Delta
t=T/M=0.001$. For each $N$, we compute the $\ell^{2}$ norm in time of the
difference between the jump $J^{N}$ with respect to the jump of the reference
solution with $75$ nodes.

Fig. 7 shows the results of this set of runs. We readily attest the
exponential rate of convergence $(error)_{N}\propto\lambda^{-2N}$ for the main
lattice spacings $\lambda=2$, $\varphi$ (the golden ratio $\varphi\approx
1.618$) and $\sigma$ (the plastic number $\sigma\approx 1.325$). Since the
convergence depends on the lattice spacing, errors reach minimum value at a
different number of node points, which is approximately $N=24$ for
$\lambda=2$, $N=35$ for $\lambda=\varphi$ and $N=60$ for $\lambda=\sigma$.

##  4 Two-dimensional boundary layers

In this section, we study the development of boundary layers over solid
boundaries and the convergence of Navier-Stokes solutions to Euler’s at
infinite Reynolds limit in this framework. Just like in reported DNS [34, 33],
we study this problem in two-dimensions. In the light of Kato’s Equivalence
Theorem [23], we track total dissipation of Navier-Stokes flows with
increasing Reynolds number. Here, however, we do not give a rigorous final
answer to this problem on logarithmic lattices. Many questions remain open and
further work is needed. The main goal of this chapter is to present promising
results for this long standing problem. Particularly, we show that logarithmic
models can reach extremely larger Reynolds numbers than in DNS.

We consider two-dimensional flow in the presence of a flat boundary at $y=0$.
Using the two-dimensional logarithmic lattice
$\mathbf{k}=(k_{x},k_{y})\in\mathbb{\Lambda}^{2}$, the velocity field
$\mathbf{u}(\mathbf{k},t)=(u,v)$ and the scalar pressure $p(\mathbf{k},t)$ are
governed by the incompressible Navier-Stokes equations with a solid boundary
(17), and can be explicitly written in the form

$\begin{cases}\partial_{t}u+u\ast\partial_{x}u+v\ast\partial_{y}u=-\partial_{x}p+\nu\Delta
u-\nu J(k_{x},t)\delta(k_{y})\quad&\text{in}\ \mathbb{\Lambda}^{2}\\\
\partial_{t}v+u\ast\partial_{x}v+v\ast\partial_{y}v=-\partial_{y}p+\nu\Delta
v\quad&\text{in}\ \mathbb{\Lambda}^{2}\\\
\partial_{x}u+\partial_{y}v=0\quad&\text{in}\ \mathbb{\Lambda}^{2}\\\
u(-k_{y})=u(k_{y}),\quad v(-k_{y})=-v(k_{y}),\quad
p(-k_{y})=p(k_{y})\quad&\text{in}\ \mathbb{\Lambda}^{2}\\\
(u,\delta)_{k_{y}}=0\quad&\text{for all}\
k_{x}\in\mathbb{\Lambda}.\end{cases}$ (63)

In two dimensions, the boundary shear force acts only in the $x$ direction.
Therefore we represent the jump variable $J$ as a scalar.

As usual, the pressure can be eliminated from the equation due to
incompressibility, by solving the Poisson equation (21). In the case of non-
zero viscosity, the pressure will be solved in terms of velocities and the
jump $J$.

#### Computation of the jump

In the numerical integrations, we consider system (63) on the truncated
logarithmic lattice $\mathbb{\Lambda}_{N}^{2}$, given by

$\mathbb{\Lambda}_{N}=\\{\pm
1,\pm\lambda,\pm\lambda^{2},\dots\pm\lambda^{N-1}\\}.$ (64)

As we did in the previous sections, we add superscript $N$ to the variables to
indicate the solutions of the truncated system.

To compute the jump $J^{N}$, we take the $k_{y}$-inner product of the $u^{N}$
equation against $\delta^{N}$ and use the no-slip condition
$(u^{N},\delta^{N})_{k_{y}}=0$ to eliminate the time derivative contribution
$(\partial_{t}u^{N},\delta)_{k_{y}}=0$. After solving for the pressure
$p^{N}$, the approximated jump $J^{N}$ can be isolated and it has the
following explicit formula

$J^{N}(k_{x},t)=\frac{1}{\nu(\xi\delta^{N},\delta^{N})}_{k_{y}}(\partial_{x}\Delta^{-1}\nabla\cdot(\mathbf{u}^{N}\ast\nabla\mathbf{u}^{N})-u^{N}\ast\partial_{x}u^{N}-v^{N}\ast\partial_{y}u^{N}+\nu\partial_{y}^{2}u^{N},\delta^{N})_{k_{y}},$
(65)

where $\Delta^{-1}$ and $\xi$ are the Fourier multipliers

$\Delta^{-1}(k_{x},k_{y})=\frac{1}{k_{x}^{2}+k_{y}^{2}},\quad\xi(k_{x},k_{y})=\frac{k_{y}^{2}}{k_{x}^{2}+k_{y}^{2}}.$
(66)

The jump $J$ is recovered from the approximated jumps $J^{N}$ throught the
limit

$J(k_{x},t)=\lim\limits_{N\to\infty}J^{N}(k_{x},t).$ (67)

#### Convergence of jump variable

Let us analyse expression (65) for the truncated jump $J^{N}$. We can write it
in terms of two contributions

$J^{N}=NL^{N}+Diss^{N},$ (68)

where $NL^{N}$ represents the contribution of nonlinear terms

$NL^{N}=\frac{1}{\nu(\xi\delta^{N},\delta^{N})}_{k_{y}}(\partial_{x}\Delta^{-1}\nabla\cdot(\mathbf{u}^{N}\ast\nabla\mathbf{u}^{N})-u^{N}\ast\partial_{x}u^{N}-v^{N}\ast\partial_{y}u^{N},\delta^{N})_{k_{y}}$
(69)

and $Diss^{N}$, the contribution of the dissipative term

$Diss^{N}=\frac{1}{(\xi\delta^{N},\delta^{N})}_{k_{y}}(\partial_{y}^{2}u^{N},\delta^{N})_{k_{y}}.$
(70)

If the velocities remain in the same regularity class (with the expected
$k^{-2}$ tail in the spectrum), the numerator of $NL^{N}$ is going to converge
to a finite value in the limit $N\to\infty$. The denominator, however, grows
with respect to $N$, for a fixed viscosity. Therefore, we must have

$\lim\limits_{N\to\infty}NL^{N}=0,$ (71)

and only $Diss^{N}$ effectively contributes to the jump in the limit
$N\to\infty$. The term $NL^{N}$ cannot be dropped, however, otherwise the no-
slip boundary condition is not satisfied.

Run | I | II | III | IV | V | VI | VII | VIII | IX | X
---|---|---|---|---|---|---|---|---|---|---
$\mathbf{Re}$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$ | $10^{5}$ | $10^{6}$ | $10^{7}$ | $10^{8}$ | $10^{9}$ | $10^{10}$
$\mathbf{N}$ | $60$ | $68$ | $75$ | $83$ | $90$ | $98$ | $105$ | $113$ | $120$ | $128$
Table 1: Reynolds number $Re$ and number of nodes $N$ for each run.

A more subtle convergence problem arises in the study of vanishing viscosity.
A small viscosity parameter appearing in the denominator of $NL^{N}$ may
amplify the spurious contribution of the nonlinearities to the computation of
jumps. To guarantee convergence, we must first take the truncation limit
$N\to\infty$, and later the vanishing viscosity limit $\nu\to 0$. Still, for
proper numerical simulations, we must have the compensation of the two terms
$\nu$ and $(\xi\delta^{N},\delta^{N})$, in such a way that their product is
big enough to provide a small $NL^{N}$. In order to guarantee an adequate
convergence rate to $NL^{N}\to 0$, we redefine the inner product as

$(f,g)=\sum_{\mathbf{k}\in\mathbb{\Lambda}^{d}}|k_{1}\dots
k_{d}|^{\alpha}f(\mathbf{k})\overline{g(\mathbf{k})}$ (72)

and the star product accordingly

$(f\ast g)(\mathbf{k})=|k_{1}\dots
k_{d}|^{\beta}\sum_{\begin{subarray}{c}\mathbf{p}+\mathbf{q}=\mathbf{k}\\\\[2.0pt]
\mathbf{p},\mathbf{q}\in\mathbb{\Lambda}^{d}\end{subarray}}\left|\frac{p_{1}\dots
p_{d}}{k_{1}\dots k_{d}}\frac{q_{1}\dots q_{d}}{k_{1}\dots
k_{d}}\right|^{\frac{\alpha+\beta}{3}}f(\mathbf{p})g(\mathbf{q}).$ (73)

The factors $|k_{1}\dots k_{d}|$ are interpreted as the volume of lattice
cells. Parameters $\alpha$ and $\beta$ can be manipulated to change
dimensionality and the scaling of terms. With such modifications, we do not
change the structure of the equations, nor their properties (symmetry group,
conservation laws, etc). Consult [7, §4] for more details.

For the generalized inner product (72) and the generalized convolution product
(73) with parameters $\alpha>0$ and $\beta=0$, one may estimate the lower
bound

$(\xi\delta^{N},\delta^{N})\geq\lambda^{\alpha(N-1)}\quad\text{for all}\
k_{x}\in\mathbb{\Lambda}_{N}.$ (74)

Since $(\xi\delta^{N},\delta^{N})$ grows exponentially fast, the decrease in
one order of magnitude in $\nu$ is roughly compensated by a linear increase in
$N$. In analogy to the Boundary Layer Theory [41], we expect the jumps $J$ to
increase as the viscosity vanishes. For this reason, in practice we do not try
to keep $NL^{N}$ small in absolute value, but in relative error with respect
to the total jump, i.e. we keep $NL^{N}/J^{N}$ small.

###  4.1 Numerical setup

We chose the logarithmic lattice with golden ratio spacing $\lambda\approx
1.618$. Model (63) is integrated with double precision by Matlab’s ode15s
solver [42]. We set the tolerances $RelTol=10^{-10}$ and $AbsTol=10^{-13}$.
For the operations on logarithmic lattices, we employ the computational
library LogLatt [5, 6].

The initial condition is fixed for all simulations and taken as random
components in a box of three by three nodes, at the large scales
$(k_{x},k_{y})$ with $1\leq k_{x},k_{y}\leq\lambda^{2}$, and zero elsewhere.
The components are adjusted to match no-slip condition and zero initial jump.

We simulate Reynolds numbers $Re$ from $10^{1}$ up to $10^{10}$. Parameter
$\alpha$ and the number of nodes $N$ in each direction of the lattice are
chosen so the relative error $\varepsilon^{N}=NL^{N}/J^{N}$ is kept small. We
set $\alpha=0.2$. The values of $Re$ and $N$ for each run are in Tab. 1. The
error $\varepsilon^{N}$ oscillates along scales due to nonlinearities, but
they decrease in average from $10^{-2}$ to $10^{-4}$ as $Re$ increases. To
ensure that the qualitative behavior do not change because of the error on the
jump computation, we performed simulations with twenty more and twenty less
nodes. All simulations presented the same qualitative behavior, with also good
quantitative agreement before the transition to chaos. Moreover, the error in
$\ell^{\infty}$ norm for the incompressibility was kept below $8\times
10^{-15}$ for all simulations at all instants, and the error in
$\ell^{\infty}_{t}\ell^{1}_{k_{x}}$ for the slip at the boundary was kept
below $1.2\times 10^{-9}$.

We now turn to the detailed presentation of the results.

###  4.2 Laminar to turbulent transition

Direct numerical simulations of dipole-wall collision [34, 33] indicate that
possible singularities at the boundary are expected to develop from sharp
gradients in the direction tangential to the wall. Such sharp gradients have
been related to the boundary-layer detachment. For this reason, we start the
analysis of results by considering the simple solution spectrum

$\mathcal{S}(k)=\sum_{k\leq|\mathbf{k}^{\prime}|<\lambda
k}|\mathbf{u}(\mathbf{k}^{\prime})|$ (75)

and the spectrum in the $x$ direction

$\mathcal{S}_{x}(k_{x})=\sum_{k_{y}\in\mathbb{\Lambda}}|\mathbf{u}(k_{x},k_{y})|.$
(76)

Fig. 8 plots the two above spectra for Run VI, whose $Re=10^{6}$. The behavior
of early and late times are distinct. The initial instants in Fig. 8(a) depict
an organized and ordered state, followed at late times by a disorganized and
chaotic state in Fig. 8(b). For this reason, we call the first regime laminar
and the second turbulent. We describe their features in details now.

The solutions are initialized at large scales only, but they instantly develop
a $k^{-2}$ tail (which actually occurs in $k_{y}$ direction only, but
dominates the whole spectrum) and a constant plateau $k^{0}$ at intermediate
scales $k$. The explanation here of such scalings is similar to what we
elaborated for one-dimensional shear flows—see the asymptotic solutions
constructed in Section 3.1. For large $k_{y}$, we have the asymptotic balance

$\nu\partial_{y}^{2}u\sim\nu J\delta,$ (77)

which gives us the $k_{y}^{-2}$ tail

$u(k_{x},k_{y},t)\sim-J(k_{x},t)k_{y}^{-2}\quad\text{for large}\ k_{y}.$ (78)

In intermediate scales and initial instants, we expect

$\partial_{t}u\sim-\nu J\delta\quad\text{for small}\ k_{y},$ (79)

whose solution provides the plateau. All these phenomena concern the boundary
effect in $y$ direction. On the other hand, the spectrum in $x$ direction
remains confined to large scales, slowly propagating to intermediate values of
$k_{x}$.

Figure 8: Velocity spectra in logarithmic scales. In cold colors, we plot the
total spectrum (75). In hot colors, we plot the spectrum (76) in $x$
direction. (a) Spectra of solution for early times $t=0,\ 0.01,\ 0.02,\
\dots,\ 0.15$. (b) Spectra of solution for late times $t=0.2,\ 0.7,\ 1.2,\
\dots,\ 10.0$.

As time advances, the plateau level increases monotonically until it reaches a
certain saturated value. Then, other terms of the governing equation (like the
nonlinear term) start to influence the solution. In Fig. 8(b) we track the
abrupt transition from the laminar to turbulent states. The smooth and clear
slopes evolve to sinuous and strong oscillations around a possible average
power-law. The $k^{-2}$ tail is still present due to the boundary effects in
$y$ direction. The plateau, however, changes from $k^{0}$ to a $k^{-0.27}$
decay. The $x$-spectrum $\mathcal{S}_{x}(k_{x})$ quickly propagates towards
high $k_{x}$ with the development of the power-law $S_{x}\propto
k_{x}^{-0.55}$. The abrupt transition is also verified in the time evolution
of the jump variable $J(k_{x},t)$, depicted in Fig. 9. As time advances, this
quantity develops the $k_{x}^{-0.27}$ scaling, the same observed in the total
spectrum $\mathcal{S}(k)$. Here, the exponents are rough approximations, and
they are plot in the figure just for reference.

Figure 9: Jump variable spectra $J(k_{x},t)$ with respect to $k_{x}$ in
logarithmic scales. Time advances linearly from red to yellow as $t=0.2,\
0.7,\ \dots,\ 10.0$, the same instants as those in Fig. 8(b).

###  4.3 Total dissipation at high Reynolds limit

Figure 10: Time evolution of the total dissipation $\varepsilon=2\nu\Omega$
computed from the flow enstrophy $\Omega(t)$. Reynolds numbers
$Re=10^{4},10^{5},10^{6},\dots,10^{10}$ increase from blue to red. The diffuse
gray line roughly separates the laminar from the turbulent regimes.

Fig. 10 shows the time evolution of total dissipation for the different runs.
From this picture, we can distinguish the two different regimes: at early
times, when the flow is laminar, the total dissipation decreases
monotonically; at later times, when it gets turbulent, humps and bumps in
dissipation break monotonicity. The transition from an organized to a
disorganized behavior is due to the appearence of chaotic oscillations in the
field variables, already described in the previous section. The exact instant
at which the transition occurs is, however, not clear from this picture. We
shall address again this question in the next subsection.

In the laminar regime, the total dissipation seems to collapse to zero in the
infinite Reynolds limit. After the flow is turbulized, a possible convergence
is not so clear, since monotonicity is broken. This is also the scenario in
DNS [26, 12]. Some authors conjecture the lack of convergence from Navier-
Stokes to Euler in view of those bursts in dissipation after the flow
turbulence is triggered [34, 33]. But simulations on logarithmic lattices have
the advantage of pushing the Reynolds number to incredibly large values.
Despite the strong oscillations observed in the turbulent regime, the total
dissipation appears to consistently decrease in average.

To verify how the solutions scale with Reynolds, we computed the maximum
dissipation and the maximum force intensity $F=\nu J$ within two time windows:
for the laminar regime, we consider the time window $t\in[0,0.06]$; for the
turbulent flow, $t\in[3.05,3.50]$. The results are plotted in Fig.11. To
explain the scaling of dissipation, we split the total dissipation into small
and large scales contributions in the following way. Let us consider the
dissipation spectrum

$\mathcal{E}(k)=\frac{\nu}{\Delta k}\sum_{k\leq|\mathbf{k}^{\prime}|<\lambda
k}|k^{\prime}_{x}k^{\prime}_{y}|^{\alpha}|\boldsymbol{\omega}(\mathbf{k}^{\prime})|^{2},\quad\text{with
}\Delta k=\lambda k-k.$ (80)

With this definition, the total dissipation rate $\varepsilon$ is obtained by
the sum

$\varepsilon=\sum_{1\leq k\leq k_{N}}\mathcal{E}(k)\Delta k.$ (81)

The large $\mathcal{E}_{l}$ and small $\mathcal{E}_{s}$ scale contributions
are expressed as

$\mathcal{E}_{l}=\sum_{1\leq k<k_{\text{cutoff}}}\mathcal{E}(k)\Delta
k\quad\text{and}\quad\mathcal{E}_{s}=\sum_{k_{\text{cutoff}\leq
k<k_{N-1}}}\mathcal{E}(k)\Delta k,$ (82)

where we set $k_{\text{cutoff}}=\lambda^{5}$. Together, they sum up the total
dissipation of the flow, i.e. $\varepsilon=\mathcal{E}_{l}+\mathcal{E}_{s}$.

Figure 11: Scaling of some global variables with respect to Reynolds numbers.
The first row shows maximum dissipation $\varepsilon$, and the second row
maximum force intensity $F=\nu J$. In the left, we consider early instants of
time $t\in[0,0.06]$, while in the right we study later times
$t\in[3.05,3.50]$.

When the flow is still laminar, Fig. 11(a) shows a clear convergence of
dissipation towards zero. In the light of Kato’s Theorem, the Navier-Stokes
solutions are converging to the Euler’s conservative solution. We can also
study what is happening at different scales of motion. For smaller Reynolds
numbers, the total dissipation is dominated by large scales, which decay to
zero as $\propto Re^{-1}$. Such scaling is expected by the following reason.
The total dissipation is given by the product $\varepsilon=2\nu\Omega$ between
the viscosity $\nu$ and the enstrophy $\Omega$. At large scales, turn-over
time is large, so there is not much variation in the enstrophy. If $\Omega$ is
approximately constant in this regime, the dissipation is going to decay as
$\varepsilon\propto\nu=Re^{-1}$. We recall that this is the rate at which the
Navier-Stokes solutions converge to Euler’s in the absence of boundaries— see
e.g. [32, Proposition 3.2]. As $Re$ increases, we have a transition around
$Re=10^{3}-10^{4}$, after which the small scales start to dominate the flow.
In this case, the dissipation is still converging to zero, but with the slower
rate $\propto Re^{-0.123}$.

After the flow is turbulized, some aspects of the latter scenario change. We
can see in Fig. 11(b) the development of a constant plateau dissipation in
intermediate values $10^{3}\lessapprox Re\lessapprox 10^{5}$ of Reynolds
number. The appearance of such plateau at this range of Reynolds numbers is
consistent with DNS [34] and it might suggest the lack of convergence to the
Euler conservative solution. State-of-art simulations, however, cannot achieve
Reynolds numbers higher than $10^{5}$, as it is easily reached by our
logarithmic lattice model. By increasing $Re$, the plateau is followed by a
consistent decrease in dissipation. This description agrees with the visual
aspect of Fig. 10, in which the total dissipation is decreasing in average.
This suggests that dissipation must vanish in the infinite Reynolds number,
but in a slow rate.

Figs. 11(c) and (d) show the corresponding scalings of the force intensity
$F=\nu J$ with respect to Reynolds. They both decrease to zero, indicating no
residual boundary shear force in the vanishing viscosity limit. However, the
aspect of the curves differ. In the laminar regime, the force intensity
converges to a clear power law $\propto Re^{-0.115}$, while in the turbulent
regime such scaling is obeyed in average only.

###  4.4 Convergence in the inviscid limit

Figure 12: Heating times that each wave number $k$ requires for the spectrum
$\mathcal{S}_{x}$ at $k$ to reach magnitude $10^{-10}$. (a) Heating times of
each $k$ for different Reynolds numbers, ranging from $10^{1}$ (in dark blue)
to $10^{10}$ (in ligth green). (b) Heating times against Reynolds number for
the two specific wave numbers $k=\lambda^{11}$ (in red) and $k=\lambda^{15}$
(in yellow). These two scales are also depicted in dashed lines with their
respective colors in (a) for reference.

The above results suggest that dissipation vanishes as Reynolds goes to
infinity. In view of Kato’s Theorem [23], this means that the Navier-Stokes
solutions converge to Euler’s in the vanishing viscosity limit. Here we
address some explanations for this result.

The flow is initially regular and confined to low $k_{x}$. A transition
occurs, and a strong front quickly fills smaller $k_{x}$ scales. Such
transition also correlates with the appearance of strong dissipation events.
This could be interpreted as the development of sharp gradients along
boundaries as the flow detaches, which is also followed by extreme events in
dissipation [34, 33, 26, 12]. Nevertheless, our computations show that the
transition from laminar to turbulence and so the developement of small
structures along the boundary—takes more time to take place as the Reynolds
number increases. Indeed, Fig. 12(a) shows the heating times that each scale
$k$ needs for the spectrum $\mathcal{S}_{x}(k)$ at $k$ to reach magnitude
$10^{-10}$. As we can see, the heating times increase with the Reynolds
number, a fact which is confirmed in Fig. 12(b) for two fixed scales
$k=\lambda^{11}$ and $k=\lambda^{15}$. This explains why, in a fixed interval
of time, the dissipation looks like is converging to a constant value, but
start vanishing again with increasing Reynolds: the transition instants drift
away from the time window for larger Reynolds numbers.

In the case of dipole-wall collision, we don’t expect the instant of
transition to advance as the Reynolds number increases. Indeed, turbulence is
triggered by the collision of the vortex against the wall, and the dipole
travel time is supposedly governed by the Reynolds-independent bulk flow [35].
We remark, however, that the results on logarithmic lattices do not disagree
with those from DNS [34, 33], if we compare total dissipation in their common
Reynolds number range. For instance, if we look at the Reynolds interval
$10^{3}\leq Re\leq 10^{5}$, the dissipation in the logarithmic lattice model
from Fig. 11(b) is similar to DNS results from [34].

As we mentioned in the description of the numerical setup, we performed some
of the simulations with twenty less and twenty more nodes, in order to check
the influence of lattice truncation on the results. The same qualitative
picture was observed for all runs. We believe, however, that more rigorous
numerical study should be undertaken.

##  5 Conclusions

We showed that logarithmic lattices are able to model systems with solid
boundaries. We explored several examples of shear flows to create intuition on
the unusual description of flows with walls in Fourier space. The vanishing
viscosity problem was approached for two-dimensional flows in the presence of
a flat solid boundary. We tracked total dissipation, in the light of Kato’s
Theorem. As time evolves, a laminar organized regime is followed by a chaotic
turbulent state. This resembles the development of small-scale structures in
the dipole-wall collision. Nevertheless, the Navier-Stokes solutions seem to
converge to Euler’s in the infinite Reynolds limit. We explained that the
propagation of energy from large to small scales in $k_{x}$ takes more and
more time with increasing Reynolds, which justifies the vanishing dissipation
on compact time intervals. We argued that this fact does not contradict the
reported DNS from literature, since the results are in agreement in the
resolution accessible to the state-of-the-art DNS. Further investigation,
however, is needed to obtain a better understanding of this challenging
problem. Still, the logarithmic lattice model looks like a promising tool for
this and other related problems, since it reaches the large Reynolds number
$Re=10^{10}$ with moderate computational effort compared to the expensive DNS
[34, 12, 33] limited to $Re=10^{5}$.

We expect that the techniques presented here will be useful for the study of
many relevant problems surrounding the effects of boundaries. This includes
the boundary layer instabilities [20], the finite-time blowup in the Prandtl
equations [15, 27] and further studies about the inviscid limit of the Navier-
Stokes equations [14]. Another problem that fits our framework is the
developement of finite-time singularities in the Euler equations with
boundaries [24], which, despite the recent successful mathematical advances
[16, 11] and the strong numerical evidences [31], is still under active
investigation [25].

The approach presented in this paper is restricted to a steady flat solid
boundary. While this is a good setup for the study of local singularities,
more general problems may involve complex geometries, unsteady walls or even
multiple boundaries. We leave these interesting but not straightforward
generalizations to future works.

## Acknowledgments

The authors express their gratitude to Bérengère Dubrulle and Simon Thalabard
for their continuous support and debate on logarithmic lattices. They are also
grateful to Marie Farge for the insights about her DNS results on the
vanishing viscosity limit. C.C. was supported by the French National Research
Agency (ANR project TILT; ANR-20-CE30-0035). A.M. was supported by CNPq grant
308721/2021-7, FAPERJ grant E-26/201.054/2022, and FAPERJ Pensa Rio grant
E-26/210.874/2014. C.C. acknowledges the finantial support from IMPA for his
visit, where the writing of this manuscript has been completed. C.C. dedicates
this work to Our Lady of Mount Carmel.

## Appendix A: Discontinuous Navier-Stokes equations

In this appendix, we provide some details and references about the Navier-
Stokes equations for flows with discontinuities.

Discontinuous formulations of the Fluid Dynamics equations were first derived
in [43] for compressible ideal flows and subsequently applied to gas dynamics
[44], compressible magnetohydrodynamics [40], and in the modelling of no-slip
boundary for incompressible viscous flow [19]. It also inspired many immersed
boundary techniques that revealed useful for the numerical simulation of
complex boundary geometries [36, 37].

Let us allow the fields to have jump discontinuities on a surface
$\mathcal{S}_{t}$ immersed on the flow. Assuming the flow to be incompressible
and that there is no penetration accross $\mathcal{S}_{t}$, the governing
equations are

$\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p+\nu\Delta\mathbf{u}+[\boldsymbol{\sigma}\cdot\mathbf{n}]\delta_{\mathcal{S}_{t}},\quad\nabla\cdot\mathbf{u}=0,$
(83)

where $\boldsymbol{\sigma}$ is the stress tensor, $\mathbf{n}$ is the unit
normal vector to $\mathcal{S}_{t}$ and the brackets $[\cdot]$ indicate the
jump accross $\mathcal{S}_{t}$. The convention for the jumps is the following:
if $\mathcal{S}_{t}$ splits the domain into two subdomains—say a positive and
a negative one—then $\mathbf{n}$ points from the positive towards the
negative, and the jump is computed as $[f]=f^{+}-f^{-}$, where $f^{+}$ is the
value at the positive side, and $f^{-}$ at the negative side. Finally, we have
also introduced the Dirac delta $\delta_{\mathcal{S}_{t}}$ locating
$\mathcal{S}_{t}$ given by

$\delta_{\mathcal{S}_{t}}(\mathbf{x})=\int_{\mathcal{S}_{t}}\prod_{i=1}^{d}\delta(x_{i}-y_{i})dS(\mathbf{y}),$
(84)

where $\delta$ is the usual Dirac delta, and $d$ is the spatial dimension.

In the framework under consideration in this paper, we can compute the term
$[\boldsymbol{\sigma}\cdot\mathbf{n}]\delta_{\mathcal{S}_{t}}$ explicitly. Our
boundary is the steady plane $y=0$, so we have simply
$\delta_{\mathcal{S}_{t}}=\delta(y)$. It remains to establish the jump
$[\boldsymbol{\sigma}\cdot\mathbf{n}]$ of tractions accross the boundary.
First, we set the upper-half space $y>0$ as being the positive domain split by
the discontinuity surface $y=0$, while the lower-half space $y<0$ is its
negative counterpart. Following the convention established above, the unit
normal vector at $y=0$ pointing from the positive towards the negative domains
is $\mathbf{n}=(0,-1,0)^{T}$. Next, the stress tensor on $y=0$ is given by

$\boldsymbol{\sigma}=\begin{pmatrix}-p&\nu\partial_{y}u&0\\\
\nu\partial_{y}u&-p+2\nu\partial_{y}v&\nu\partial_{y}w\\\
0&\nu\partial_{y}w&-p\end{pmatrix}\quad\text{on}\ y=0,$ (85)

where we have used the no-slip boundary condition $\mathbf{u}=\mathbf{0}$ on
the $x$-$z$ plane to neglect all terms involving spatial derivatives in $x$
and $z$ directions, that is
$\partial_{x}u=\partial_{x}v=\partial_{x}w=\partial_{z}u=\partial_{z}v=\partial_{z}w=0$
on $y=0$. Hence,

$\boldsymbol{\sigma}\cdot\mathbf{n}=\begin{pmatrix}-\nu\partial_{y}u\\\
p-2\nu\partial_{y}v\\\ -\nu\partial_{y}w\end{pmatrix}\quad\text{on}\ y=0.$
(86)

Finally, because of the symmetries (2) on $p$ and $v$, we have
$[p]=[\partial_{y}v]=0$, so we obtain

$[\boldsymbol{\sigma}\cdot\mathbf{n}]=-\nu\begin{pmatrix}[\partial_{y}u]\\\
0\\\ [\partial_{y}w]\end{pmatrix}.$ (87)

This leads to system (4)-(5).

## Appendix B: Boundary conditions and nonuniqueness of immersed boundary
flows

The purpose of this appendix is to show that we cannot suppress the boundary
condition in the immersed boundary models, otherwise the problem becomes ill-
posed. More precisely, we show that if we do not impose the no-slip condition,
solutions are nonunique. For sake of simplicity, we do it for one-dimensional
shear flows.

Let us consider the initial value problem

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u-\nu
J(t)\delta(y)\quad&\text{in}\ y\in\mathbb{R},\ t>0\\\
u(y,t)=u(-y,t)\quad&\text{in}\ y\in\mathbb{R},\ t>0\\\
u(y,0)=g(y)\quad&\text{for}\ y\in\mathbb{R}.\end{cases}$ (88)

with the jump $J(t)$ given by

$J(t)=[\partial_{y}u].$ (89)

This is the immersed boundary formulation of the one-dimensional unforced
shear flow (34), but here the boundary condition $u=0$ at $y=0$ is not
imposed. Using the classical theory of heat equation, we shall construct two
distinct solutions for problem (88).

Let us define

$u(y,t)=\int_{0}^{\infty}\left(K(y,x,t)-K(y,-x,t)\right)g(x)dx,\quad\text{for}\
y>0$ (90)

and $u(-y,t)=u(y,t)$ for all $y\in\mathbb{R}$, where

$K(y,x,t)=\frac{e^{-|y-x|^{2}/4\nu t}}{(4\pi\nu t)^{1/2}}$ (91)

is the heat kernel [21]. We claim that such function solves problem (88).
Indeed, (90) is the solution of the Dirichlet problem

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u\quad&\text{in}\ y>0,\ t>0\\\
u(y,0)=g(y)\quad&\text{for}\ y>0\\\ u(0,t)=0\quad&\text{for}\ t>0\
\text{(Dirichlet B.C.)}.\end{cases}$ (92)

Therefore, the constructed $u(y,t)$ is the physical solution of problem (88),
i.e. the solution which satisfies the no-slip condition $u=0$ at $y=0$.

Figure 13: Two distinct solutions for the immersed boundary model (88) from
the same intial condition (in black). Initial data satisfy no-slip $u=0$ at
$y=0$. Colors change from blue to green as time advances. (a) solution with
no-slip at the boundary, obtained from imposing Dirichlet boundary condition
$u=0$ at $y=0$; (b) solution with slip at the boundary, satisfying the Neumann
boundary condition $\partial_{y}u=0$ at $y=0$.

Next, define

$u(y,t)=\int_{0}^{\infty}\left(K(y,x,t)+K(y,-x,t)\right)g(x)dx,\quad\text{for}\
y>0$ (93)

and $u(-y,t)=u(y,t)$ for all $y\in\mathbb{R}$. We claim that this is also a
solution of problem (88). Indeed, (93) is the solution of the Neumann problem

$\begin{cases}\partial_{t}u=\nu\partial_{y}^{2}u\quad&\text{in}\ y>0,\ t>0\\\
u(y,0)=g(y)\quad&\text{for}\ y>0\\\ \partial_{y}u(0,t)=0\quad&\text{for}\ t>0\
\text{(Neumann B.C.)}.\end{cases}$ (94)

Therefore, the constructed $u(y,t)$ by the reflection trivially satisfies
problem (88), since the jump variable is identically zero $J\equiv 0$. So we
have constructed two distinct solutions to for problem (88).

In Fig. 13 we show the two distinct solutions constructed above for the
initial condition

$g(y)=C\exp\left(-\frac{1}{1-|y-y_{0}|^{2}}\right),$ (95)

with $y_{0}=1.2$ and $C$ such that $\max_{y}g(y)=1$. Such function has compact
support $\text{supp}\ g=[0.2,2.2]$; particularly, $g(0)=0$. With this example,
we also show that it is not sufficient to have no-slip condition at the
initial instant for the solution to satisy no-slip at all subsequent instants.
In Fig. 13(a) we see the physical expected solution, which satisfies no-slip
condition at all times. Differently, Fig. 13(b) shows a solution which
initially satisfies $u=0$ at $y=0$, but develops non zero slip as time
advances.

## References

  * [1] C. W. Bardos and E. S. Titi.  Mathematics and turbulence: where do we stand?  J. Turbul., 14(3):42–76, 2013.
  * [2] A. Barral and B. Dubrulle.  Asymptotic ultimate regime of homogeneous Rayleigh–Bénard convection on logarithmic lattices.  J. Fluid Mech., 962:A2, 2023.
  * [3] L. Biferale.  Shell models of energy cascade in turbulence.  Ann. Rev. Fluid Mech., 35:441–468, 2003.
  * [4] C. S. Campolina.  Fluid Dynamics on Logarithmic Lattices and Singularities of Euler Flow.  Master’s Thesis, Instituto de Matemática Pura e Aplicada, 2019.
  * [5] C. S. Campolina.  Loglatt: A computational library for the calculus and flows on logarithmic lattices.  arXiv preprint arXiv:2006.00047, 2020.
  * [6] C. S. Campolina.  LogLatt: A computational library for the calculus on logarithmic lattices, 2020.  Freely available for noncommercial use from MATLAB Central File Exchange (https://www.mathworks.com/matlabcentral/fileexchange/76295-loglatt).
  * [7] C. S. Campolina.  Fluid Flows and Boundaries on Logarithmic Lattices.  PhD Thesis, Instituto de Matemática Pura e Aplicada, 2022.
  * [8] C. S. Campolina and A. A. Mailybaev.  Chaotic blowup in the 3D incompressible Euler equations on a logarithmic lattice.  Phys. Rev. Lett., 121:064501, 2018.
  * [9] C. S. Campolina and A. A. Mailybaev.  Fluid dynamics on logarithmic lattices.  Nonlinearity, 34(7):4684, 2021.
  * [10] K. W. Cassel.  A comparison of Navier-Stokes solutions with the theoretical description of unsteady separation.  Philos. Trans. Roy. Soc. A, 358(1777):3207–3227, 2000.
  * [11] J. Chen and T. Y. Hou.  Finite time blowup of 2D Boussinesq and 3D Euler equations with $C^{1,\alpha}$ velocity and boundary.  Comm. Math. Phys., 383(3):1559–1667, 2021.
  * [12] H. J. H. Clercx and G. J. F. Van Heijst.  Dissipation of coherent structures in confined two-dimensional turbulence.  Physics of Fluids, 29(11):111103, 2017.
  * [13] G. Costa, A. Barral, and B. Dubrulle.  Reversible Navier-Stokes equation on logarithmic lattices.  Phys. Rev. E, 107(6):065106, 2023.
  * [14] W. E.  Boundary layer theory and the zero-viscosity limit of the Navier-Stokes equation.  Acta Math. Sin., 16:207–218, 2000.
  * [15] W. E and B. Engquist.  Blowup of solutions of the unsteady Prandtl’s equation.  Comm. Pure Appl. Math, 50(12):1287–1293, 1997.
  * [16] T. M. Elgindi and I.-J. Jeong.  Finite-time singularity formation for strong solutions to the axi-symmetric 3D Euler equations.  Ann. PDE, 5(2):16, 2019.
  * [17] F. G. Friedlander.  Introduction to the Theory of Distributions.  Cambridge University Press, 1998.
  * [18] D. Gérard-Varet and E. Dormy.  On the ill-posedness of the Prandtl equation.  J. Amer. Math. Soc., 23(2):591–609, 2010.
  * [19] D. Goldstein, R. Handler, and L. Sirovich.  Modeling a no-slip flow boundary with an external force field.  Journal of Computational Physics, 105(2):354–366, 1993.
  * [20] E. Grenier.  On the nonlinear instability of Euler and Prandtl equations.  Comm. Pure Appl. Math., 53(9):1067–1091, 2000.
  * [21] F. John.  Partial Differential Equations, volume 1 of Applied Mathematical Sciences.  Springer, 3rd edition, 1978.
  * [22] T. Kato.  Nonstationary flows of viscous and ideal fluids in $R^{3}$.  J. Funct. Anal., 9(3):296–305, 1972.
  * [23] T. Kato.  Remarks on zero viscosity limit for nonstationary Navier-Stokes flows with boundary.  In Seminar on nonlinear partial differential equations, pages 85–98. Springer, 1984.
  * [24] A. Kiselev.  Small scales and singularity formation in Fluid Dynamics.  In Proc. Int. Cong. of Math. 2018, volume 2, pages 2329–2356, Rio de Janeiro, 2018. World Scientific.
  * [25] S. S. V. Kolluru, P. Sharma, and R. Pandit.  Insights from a pseudospectral study of a potentially singular solution of the three-dimensional axisymmetric incompressible Euler equation.  Phys. Rev. E, 105(6):065107, 2022.
  * [26] W. Kramer, H. J. H. Clercx, and G. J. F. Van Heijst.  Vorticity dynamics of a dipole colliding with a no-slip wall.  Phys. Fluids, 19(12):126603, 2007.
  * [27] I. Kukavica, V. Vicol, and F. Wang.  The van Dommelen and Shen singularity in the Prandtl equations.  Adv. Math., 307:288–311, 2017.
  * [28] L. D. Landau and E. M. Lifshitz.  Fluid Mechanics, volume 6.  Pergamon, 1987.
  * [29] R. J. LeVeque.  Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems.  SIAM, 2007.
  * [30] H. J. N. Lopes and M. C. Lopes Filho.  Fluids, walls and vanishing viscosity.  In Proc. Int. Cong. Math. 2018, volume 3, pages 2537–2558, Rio de Janeiro, 2018. World Scientific.
  * [31] G. Luo and T. Y. Hou.  Potentially singular solutions of the 3D incompressible Euler equations.  PNAS, 111:12968–1297349, 2014.
  * [32] A. J. Majda and A. L. Bertozzi.  Vorticity and Incompressible Flow.  Cambridge University Press, 2002.
  * [33] N. Nguyen van yen, M. Waidmann, R. Klein, M. Farge, and K. Schneider.  Energy dissipation caused by boundary layer instability at vanishing viscosity.  J. Fluid Mech., 849:676–717, 2018.
  * [34] R. Nguyen van yen, M. Farge, and K. Schneider.  Energy dissipating structures produced by walls in two-dimensional flows at vanishing viscosity.  Phys. Rev. Lett., 106:184502, 2011.
  * [35] P. Orlandi.  Vortex dipole rebound from a wall.  Physics of Fluids A: Fluid Dynamics, 2(8):1429–1436, 1990.
  * [36] C. S. Peskin.  Flow patterns around heart valves: a numerical method.  Journal of computational physics, 10(2):252–271, 1972.
  * [37] C. S. Peskin.  The immersed boundary method.  Acta numerica, 11:479–517, 2002.
  * [38] Q. Pikeroen, A. Barral, G. Costa, C. Campolina, A. Mailybaev, and B. Dubrule.  Tracking complex singularities of fluids on log-lattices.  Nonlinearity (submitted), 2023.
  * [39] Q. Pikeroen, A. Barral, G. Costa, and B. Dubrulle.  Log-lattices for atmospheric flows.  Atmosphere, 14(11):1690, 2023.
  * [40] E. P. Salathé and L. Sirovich.  Boundary-value problems in compressible magnetohydrodynamics.  The Physics of Fluids, 10(7):1477–1491, 1967.
  * [41] H. Schlichting and J. Kestin.  Boundary Layer Theory, volume 121.  Springer, 1961.
  * [42] L. F. Shampine and M. W. Reichelt.  The Matlab ODE suite.  SIAM J. Sci. Comput., 18(1):1–22, 1997.
  * [43] L. Sirovich.  Initial and boundary value problems in dissipative gas dynamics.  The Physics of Fluids, 10(1):24–34, 1967.
  * [44] L. Sirovich.  Steady gasdynamic flows.  The Physics of Fluids, 11(7):1424–1439, 1968.
  * [45] L. N. Trefethen.  Spectral methods in MATLAB.  SIAM, 2000.