Computing a Class of Blow-up Solutions for the Navier-Stokes Equations
C. BOLDRIGHINI1, S. FRIGIO2, P. MAPONI2, A. PELLEGRINOTTI3
1Istituto Nazionale di Alta Matematica (INdAM), Gruppo Nazionale per la Fisica Matematica (GNFM)
Università di Roma La Sapienza
Piazzale Aldo 2 Moro, 00185 Rome
ITALY
2Scuola di Scienze e Tecnologie
Università di Camerino
62032 Camerino (MC)
ITALY
3Dipartimento di Matematica e Fisica
Università di Roma Tre
Largo S. Leonardo Murialdo 1, 00146 Rome
ITALY
Abstract: The three-dimensional incompressible Navier-Stokes equations play a fundamental role in a large num-
ber of applications to fluid motions, and a large amount of theoretical and experimental studies were devoted to
it. Our work is in the context of the Global Regularity Problem, i.e., whether smooth solutions in the whole space
R3can become singular (“blow-up”) in a finite time. The problem is still open and also has practical importance,
as the singular solutions would describe new phenomena. Our work is mainly inspired by a paper of Li and Sinai,
who proved the existence of a blow-up for a class of smooth complex initial data. We present a study by computer
simulations of a larger class of complex solutions and also of a related class of real solutions, which is a natural
candidate for evidence of a blow-up. The numerical results show interesting features of the solutions near the
blow-up time. They also show some remarkable properties for the real flows, such as a sharp increase of the total
enstrophy and a concentration of high values of velocities and vorticity in small regions.
Key-Words: Navier-Stokes equations, Blow-up, Global Regularity Problem, Numerical solution, Fluid dynamics
Received: June 17, 2023. Revised: May 16, 2024. Accepted: July 11, 2024. Published: August 8, 2024.
1 Introduction
Fluid dynamics studies fluid motion and its interac-
tion with the environment, such as solid bodies and
other fluids. It has a wide variety of applications:
calculating forces and moments, [1], determining the
mass flow rate of oil through pipelines, [2], fore-
casting weather patterns, [3], designing aircrafts, [4],
and turbine machines, [5], manage renewable energy
sources, [6], [7], optimize processes in nutrition sci-
ence, [8], [9], [10], [11], only to mention a few. The
relevance of the applications stimulated an intense
study of fluid dynamics from different points of view,
from approximation and computational methods to
obtain particular solutions to theoretical frameworks
providing general properties of fluid dynamics mod-
els.
The incompressible Navier-Stokes equations de-
scribe the evolution of the flow velocity and pres-
sure and they may well be considered the main model
of fluid dynamics, [12]. We restrict our attention to
flows in the whole space R3with no boundary con-
ditions. In spite of great efforts, important properties
of the model remain to be established by the scientific
research. The Global Regularity Problem (GRP), i.e.,
the problem whether smooth solutions in absence of
forcing can become singular at a finite time, is still
open and in the list of the Clay millennium prize,
[13]. In the context of well definitness, [14], proved
a global weak existence theorem and a uniqueness
and regularity theorem only for finite times. The in-
terest of possible singularities in the solution is that
they would describe sudden concentrations of energy
in a finite region (as it happens in tornadoes or hurri-
canes); in fact, the main features of the possible finite-
time singularities (“blow-up”), are the divergence of
the total enstrophy, [12], and the divergence at some
point of the absolute value of the velocity, [15]. We
note that the study of GRP is an extremely active field
of study, where the Navier-Stokes equations is anal-
ysed theretically and numerically, see [16], [17], [18],
[19], [20], [21].
Recently, [22], introduced a new approach to
the study of incompressible Navier-Stokes equations.
They considered the integral formulation of the equa-
tions in Fourier transform space for a particular class
of initial data corresponding to complex solutions,
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
247
Volume 19, 2024
and could prove a finite-time blow-up by a renormal-
ization group method. Following this initial contri-
bution, several papers shed light on the properties of
the solutions by theoretical studies, [23], and by com-
puter simulations, [24], [25]. In this paper, we de-
scribe the integral formulation of the Navier-Stokes
equations, and the numerical approximation scheme
for their simulation. We then report the main results
obtained by a high performance computing code, both
for some complex and some related real solutions.
The organization of the paper is as follows. In the
next section, we provide a detailed derivation of the
integral formulation of the Navier-Stokes equations
in the Fourier transform space. In section 3, we de-
scribe the approximation scheme of the integral equa-
tion arising from the Navier-Stokes equations. In sec-
tion 4, we show the results obtained in numerical ex-
periments with the complex and real solutions of the
Navier-Stokes equations. In section 5, we provide
some conclusions and final remarks.
2 Integral formulation of the
Navier-Stokes equations
We consider the initial value problem for the incom-
pressible Navier-Stokes equations in the whole space
with no boundary conditions:
∂u
∂t +
3
X
j=1
uj
∂u
∂xj
= ∆u− ∇p,
x= (x1, x2, x3)t∈R3, t > 0,(1)
∇tu= 0,x∈R3, t > 0,(2)
u(x,0) = u0(x),x∈R3,(3)
where u= (u1, u2, u3)tis the velocity vector, pis the
pressure, ∇denotes the gradient operator, ∇tdenotes
the divergence operator, ∆denotes the Laplacian op-
erator, and u0is the initial velocity of the fluid. The
viscosity in (1) is ν= 1, which can be always ob-
tained by space scaling, since no boundaries are pre-
scribed.
We define the Fourier transform ^u of uas follows:
^u(k, t) = (Fu(·, t))(k)
=1
(2π)3ZR3
u(x, t)e−ıhk,xidx,k∈R3,
where h·,·i is the scalar product in R3. The integral
formulation considers the function v(k, t) = ı^u(k, t),
so from standard arguments on Fourier transform the-
ory we have:
u(x, t)=−ı(F−1v(·, t))(x)=−ıZR3
v(k, t)eıhk,xidk.
(4)
Let f, g :R→Rbe two regular functions, we remind
two well-known properties of the Fourier transform:
F(Df) = ik ˆ
f, (5)
F(fg) = ˆ
f∗ˆg, (6)
where Df denotes the derivative of fand
(f∗g)(x) = ZR
f(x0)g(x−x0)dx0(7)
denotes the convolution integral. We note that from
(5), (6), the symmetry of the convolution integral and
the definition of function vwe have
F
3
X
j=1
uj(·, t)∂u
∂xj
(·, t)
(k) =
3
X
j=1
ˆuj(·, t)∗(ıkj^u(·, t))
(k) =
ZR3
3
X
j=1
vj(k−k0, t)ık0
jv(k0, t)dk0=
−ıZR3
hv(k−k0, t),k0iv(k0, t)dk0.(8)
We apply the Fourier transform to the Navier-Stokes
equations (1), and from relations (5)-(8) we obtain
∂v(k, t)
∂t +kkk2v(k, t) =
ZR3
hv(k−k0, t),k0iv(k0, t)dk0−ıkˆp(k, t),(9)
where ˆpis the Fourier transform of p. From the in-
compressibility equation (2) we see that uis a rotation
(or solenoidal) field, so that in the Fourier space:
F∇tu(·, t)(k) = ıhk,^u(k, t)i= 0.(10)
So, in the Fourier space, the incompressibility equa-
tion (2) reduces to the orthogonality of ^u(k, t)to k,
hence, the projectior on the subspace of the rotation
fields coincides with the projector Pkon the subspace
of vectors orthogonal to k, so that for a generic vector
w∈R3it is defined as:
Pkw=w−hw,ki
kkk2k.(11)
We can easily see that Pkhas no effect on the func-
tion ^u(k, t), that is Pk^u(k, t) = ^u(k, t); on the other
hand, it cancels the term ıkˆp(k, t)in formula (9), in
fact Pkk=0. Thus, when we apply Pkto equation
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
248
Volume 19, 2024
(9) we get
∂v(k, t)
∂t +kkk2v(k, t) =
ZR3
hv(k−k0, t),k0iPkv(k0, t)dk0.(12)
We multiply both sides of (12) by etkkk2obtaining on
the right the time-derivative of etkkk2v(k, t)and, by a
time integration, the following equation:
v(k, t)=e−tkkk2v0(k) + Zt
0
e−(t−s)kkk2C(k, s;v)ds,
(13)
where
C(k, s;v) = ZR3
hv(k−k0, s),k0iPkv(k0, s)dk0,
(14)
is a convolution integral and the initial data v0are the
Fourier transform of the data u0in (3).
Equation (13) provides the integral formulation of
the Navier-Stokes equations. We note that the study of
Li and Sinai considers real solutions of this equation,
which correspond in general to complex solutions of
the original Navier-Stokes equations (1). However, if
the initial data v0, and hence the solution of equation
(13), is antisymmetric, (i.e. v(−k, t) = −v(k, t)),
then the solution in the physical space is real and de-
scribes a fluid flow.
3 The approximation scheme
We consider the following problem: given v0, com-
pute a numerical approximation V(k, t)of v(k, t)for
k∈R⊂R3, and t∈T⊂[0, t], where R,Tare
suitable discrete sets.
We note that equation (13) has a structure of a
Volterra integral equation in the time variable and a
Fredholm integral equation in the space (conjugate)
variables, when a suitable truncation of R3is con-
sidered. So, roughly speaking, the proposed approxi-
mation is based on a time-marching scheme where in
each step an iterative solution of the nonlinear Fred-
holm equation is computed, see [26], for a general
introduction to numerical approximation schemes for
integral equations.
We consider a uniform partition T={tn=δtn,
n= 0,1, . . . , N}of the interval [0, t]with step size
δt=t
N. Knowing the solution v(k, tn),k∈R3at
time tn, we can formally compute v(k, tn+1)for k∈
R3by using the rectangle quadrature formula
v(k, tn+1) = e−δtkkk2v(k, tn) +
Ztn+1
tn
e−kkk2(tn+1−s)C(k, s;v)ds
≈e−δtkkk2v(k, tn)+δte−δtkkk2C(k, tn;v),(15)
or by using the trapezoidal quadrature formula
v(k, tn+1)≈e−δtkkk2v(k, tn) +
δt
2e−δtkkk2C(k, tn;v) + C(k, tn+1;v).(16)
For the discretization of the spatial variables, we sup-
pose that, for t∈T,v(k, t)is negligible if kis out-
side the parallelepiped [a1, b1]×[a2, b2]×[a3, b3],
and consider a uniform mesh on it: R={ki,j,l =
(a1+δ1i, a2+δ2j, a3+δ3l)t,i= 0,1, . . . , I,
j= 0,1, . . . , J,l= 0,1, . . . , L}with step sizes
δ1=b1−a1
I,δ2=b2−a2
J,δ3=b3−a3
L, respectively.
Let Ibe the set of indices defined in R, the approx-
imation Vn={vi,j,l,n,(i, j, l)∈I}of v(ki,j,l, tn),
ki,j,l ∈R,tn∈Tcan be computed by formulas (15)
and (16) applying the trapezoidal quadrature formula
with nodes R, that is
vi,j,l,n+1 =e−δtkki,j,lk2vi,j,l,n +
δte−δtkki,j,lkc(i, j, l, Vn),(17)
vi,j,l,n+1 =e−δtkki,j,lk2vi,j,l,n +
δt
2e−δtkki,j,lk2c(i, j, l, Vn)+c(i, j, l, Vn+1)
(18)
where
c(i, j, l, Vm) =
X
(p,q,r)∈I
Wp,q,rhki,j,l,v(ki,j,l −kp,q,r, tm)i ·
Pki,j,l vp,q,r,m,(19)
is a discrete convolution term, and Wp,q,r,(p, q, r)∈
Iare the trapezoidal quadrature weights. We note that
the function c defined in (19) can be efficiently com-
puted by using the FFT algorithm.
Formula (17) gives an explicit Euler discretization
scheme for equation (13), in fact it explicitely gives
variables Vn+1 in terms of Vn. Formula (18) gives
an implicit discretization scheme for equation (13),
where variables Vn+1 can be computed from Vnby
solving a nonlinear equation; this last scheme is usu-
ally called Crank-Nicolson scheme. By a usual strat-
egy, these schemes are jointly used to profitably ex-
ploit the stability features of the implicit methods and,
in the particular case of the Crank-Nicolson scheme,
also to obtain a second order accuracy with respect
to the discretization steps in the time and spatial vari-
ables.
Algorithm 1 Let Tand Rbe the aforementioned dis-
cretization sets for the time variable and space vari-
ables, respectively; let n= 0,1, . . . , N and Ithe cor-
responding sets of indices. Let tol > 0and ν∈Nbe
given tolerances. Let V0={v(ki,j,l,0),(i, j, l)∈I}.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
249
Volume 19, 2024
Compute the approximation Vn,n= 1,2, . . . , N of
the solution of (13) at t∈Tby performing the fol-
lowing steps:
i) For n= 0,1, . . . , N −1,
ii) v(0)
i,j,l,n+1=e−δtkki,j,lk2vi,j,l,n+δte−δtkki,j,l k2·
c(i, j, l, Vn),(i, j, l)∈I
iii) ν= 0,
iv) Repeat
v) ν=ν+ 1,
vi) v(ν)
i,j,l,n+1 =e−δtkki,j,lk2vi,j,l,n +
δt
2e−δtkki,j,lk2c(i, j, l, Vn) +
c(i, j, l, V(ν−1)
n+1 ),(i, j, l)∈I,
vii) Until
V(ν)
n+1 −V(ν−1)
n+1
< tol or ν > ν do
viii) Vn+1 =V(ν)
n+1.
We note that Algorithm 1 computes an initial approx-
imation V(0)
nof the solution at time tnby using the
Euler method (17) and recursively refines such ini-
tial guess by using the Crank-Nicolson method (18).
This iterative process terminates when two consecu-
tive solution vectors are sufficiently close or a maxi-
mum number of iterations is reached. Moreover, from
standard arguments on Fourier transform theory, we
have that the computational cost of term Cin steps
ii) and iv) can be considered linear in the number of
variables in the arrays Vn.
4 Numerical experiment
We describe some results of a numerical experiment
performed by Algorithm 1, implemented in FOR-
TRAN 90. The MPI software library, [27], is used to
take advantage of parallel computing architectures to
enhance the performance of the computational code.
The 2Decomp&FFT library, [28], is used to compute
the discrete convolution term c in (19).
The results reported in this section have two main
aims: i) showing the accuracy of Algorithm 1, ii) un-
derlying some useful properties of the solutions of
equation (13); they are reported in two sections. Sec-
tion 4.1 analyses the accuracy of Algorithm 1; in par-
ticular, we consider one set of initial data and com-
pare the numerical solutions obtained by using differ-
ent discretization parameters. In section 4.2 we con-
sider different initial data by highlighting their blow-
up properties.
4.1 Accuracy of the discretization scheme
Let Abe a positive constant, and B⊂R3be the
sphere of center k0= (0,0, c)and radius 0< r < |c|,
let
Y=k−k0
p|k0|,
and
v(k,0)=(−AY1, Y2,−k1Y1+k2Y2
k3t
e−|Y|2
2,k∈B,
0,k6∈ B.
(20)
In particular, we choose c= 20,r= 17, and Asuch
that the energy of the initial data is
1
2ZR3
|v(k,0)|2dk= 200.
The discretization of (13) is obtained by using the
same discretization step δk>0along the three co-
ordinate directions, i.e. δk=δ1=δ2=δ3.
Table 1: Difference in the approximated solutions
computed with different discretization steps δk: first
row δk= 1,0.5, second row δk= 0.5,0.25;
the results are computed by using the space interval
[−127,127] ×[−127,127] ×[−19,1004],δt= 10−6
and they refer to the first time-iterate. The notation
x(y)stays for x·10y.
δk123
1-0.5 9.4(−7) 5.6(−9) 2.5(−7)
0.5-0.25 3.1(−8) 5.6(−9) 2.6(−8)
Table 2: Difference in the approximated solutions
computed with different discretization steps δt: first
row δt= 10−6,0.5 10−6, second row δt=
0.5 10−6,0.25 10−6; the results are computed by us-
ing the space interval [−127,127] ×[−127,127] ×
[−19,1004],δk= 1 and they refer to the first time-
iterate. The notation x(y)stays for x·10y.
δt123
1(−6) −0.5(−6) 2.3(−7) 4.7(−17) 1.2(−11)
0.5(−6) −0.25(−6) 1.2(−14) 1.2(−17) 5.7(−13)
The numerical results show the difference between
the approximated solutions computed by using dif-
ferent discretization parameters; more precisely, three
points are considered: k1= (−2,5,22)t,k2=
(10,−20,18)t,k3= (−5,10,40t), and the corre-
sponding errors l,l= 1,2,3are reported. In par-
ticular, Table 1 shows the differences in the numerical
solutions obtained with different spatial discretization
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
250
Volume 19, 2024
Table 3: Difference in the approximated solutions
computed with different space intervals [−63,63] ×
[−63,63] ×[−19,501],[−127,127] ×[−127,127] ×
[−19,1004]; the results are computed by using δk=
1,δt= 10−6and they refer to the time-iterate n=
1,195. The notation x(y)stays for x·10y.
Iterate 123
12.3(−7) 7.9(−22) 8.0(−18)
195 3.3(−15) 5.6(−15) 3.2(−15)
steps δk. Table 2 shows the differences in the numer-
ical solutions obtained with different time discretiza-
tion steps δt. Table 3 shows the differences in the nu-
merical solutions obtained with different space inter-
vals and at different time-iterates, i.e., n= 1,195.
4.2 Blow-up properties of Li Sinai solutions
We analyse four cases. The first two are a complex
flow with initial data given in section 4.1 (see formula
(20)) and the corresponding case with antisymmetric
initial data, which is a real solution of the problem
(1)-(3). The other two cases are similar, except that
in the initial data we set c= 30, so that the case with
antisymmetric initial data is again a real solution of
the problem (1)-(3).
In all cases, the results of the numerical simula-
tions are described in terms of the following quanti-
ties: the total energy:
E(t) = 1
2ZR3
|v(k, t)|2dk,(21)
the total enstrophy:
S(t) = ZR3
kkk2|v(k, t)|2dk,(22)
and the marginal density for the enstrophy along the
k3axis in k-space;
S3(k3, t) = ZR×R
kkk2|v(k, t)|2dk1dk2.(23)
Observe that for the complex solutions there is a
divergence of both energy and enstrophy, as predicted
by the Li-Sinai theory. We note that with the initial
data v0, centered along the k3-axis, the support of the
solution vof (13) extends along the k3axis, and the
features of the blow-up are well described by the func-
tion S3.
As already mentioned, the first example consid-
ers the initial data (20) in section 4.1. The nu-
merical results are computed in the space interval
500 1000 1500 t107
1000
104
105
106
Log(E)
Fig1: First case: the total energy E(t),t∈
[0,1550 10−7]obtained with the initial data (20), i.e.,
a complex solution of (1)-(3) in the class of Li and
Sinai solutions; the scale on the ordinate axis is loga-
rithmic.
500 1000 1500 t107
105
108
1011
Log(S)
Fig2: First case: the total enstrophy S(t),t∈
[0,1550 10−7]obtained with the initial data (20), i.e.,
a complex solution of (1)-(3) in the class of Li and
Sinai solutions; the scale on the ordinate axis is loga-
rithmic.
[−127,127] ×[−127,127] ×[−19,2528], time inter-
val [0,1550 10−7], and discretization steps δ1=δ2=
δ3= 1,δt= 10−7. The results are reported in Figure
1, Figure 2, Figure 3. In particular, Figure 1 shows
the plot of the total energy E(t),t∈[0,1550 10−7];
Figure 2 shows the plot of the total enstrophy S(t),
t∈[0,1550 10−7]; Figure 3 shows the plot of the
marginal density of the enstrophy along the k3axis in
k-space S3(k3, t),k3∈[−19,2528] for three times
near the blow-up, i.e. t= 1450 10−7,1500 10−7,
1530 10−7.
The results are quite satisfactory, in that they
clearly show the blow-up properties of the com-
plex solution of problem (1)-(3) under consideration.
In particular, Figure 3 gives interesting information
about the structure of the solution along the k3-axis.
The second example considers an antisymmetric
initial data of (13) providing a real solution of prob-
lem (1)-(3). More precisely, with the same notation
introduced in section 4.1, we define:
˜
v(k,0)=(−AY1, Y2,−k1Y1+k2Y2
k3te−|Y|2
2,k∈B,
0,k6∈ B,
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
251
Volume 19, 2024
t107=1530,
t107=1500
t107=1450
500 1000 1500 2000 2500
k3
10 000
20 000
30 000
S3(k3,t)
Fig3: First case: the marginal density for the en-
strophy along the k3axis in k-space S3(k3, t),k3∈
[−19,2528] for three times near the blow-up, i.e. t=
1450 10−7,1500 10−7,1530 10−7obtained with the
initial data (20), i.e., for a complex solution of (1)-(3)
in the class of Li and Sinai solutions.
Fig4: Second case: the total energy E(t),t∈
[0,2.19 10−5]obtained with the initial data (24), i.e.,
for a real solution of (1)-(3).
v(k,0) = ˜
v(k,0) −˜
v(−k,0),(24)
where again c= 20,r= 17, and Ais such that the
total initial energy E(0) = 250000.
The numerical results are computed by using
the space interval [−177,177] ×[−177,177] ×
[−1200,1200], with discretization steps δ1=δ2=
δ3= 1. The discretization step in time is δt=
1.5625 10−8for a maximal time 1400, so that the
resulting time interval is about [0,2.19 10−5]. The
results are reported in Figure 4, Figure 5, Figure 6,
Figure 13. In particular, Figure 4 shows the plot of the
200 400 600 800 1000 1200 1400 t1.562108
5.0 ×107
1.0 ×108
1.5 ×108
2.0 ×108
2.5 ×108
S
Fig5: Second case: the total enstrophy S(t),t∈
[0,2.19 10−5]obtained with the initial data (24), i.e.,
for a real solution of (1)-(3).
t*1.562 108=500
t*1.562 108=1100
t*1.562 1081400
-1000 -500 500 1000
k3
2×106
4×106
6×106
8×106
1×107
S3(k3,t)
Fig 6: Second case: the marginal density for the
enstrophy along the k3axis in k-space S3(k3, tn),
k3∈[−1200,1200] for three time-iterates near
the maximum of the total enstrophy, i.e. n=
500,1100,1400 obtained with the initial data (24),
i.e., for a real solution of (1)-(3).
Fig7: Third case: the total energy E(t),t∈
[0,1550 10−7]obtained with the initial data (20) with
c= 30, i.e., a complex solution of (1)-(3) in the class
of Li and Sinai solutions; the scale on the ordinate axis
is logarithmic.
Fig8: Third case: the total enstrophy S(t),t∈
[0,1550 10−7]obtained with the initial data (20) with
c= 30, i.e. a complex solution of (1)-(3) in the class
of Li and Sinai solutions; the scale on the ordinate axis
is logarithmic.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
252
Volume 19, 2024
t107=1075,
t107=1060
t107=1040
500 1000 1500 2000 2500
k3
20 000
40 000
60 000
80 000
S3(k3,t)
Fig9: Third case: the marginal density for the en-
strophy along the k3axis in k-space S3(k3, t),k3∈
[−19,2528] for three times near the blow-up, i.e. t=
1040 10−7,1060 10−7,1075 10−7obtained with the
initial data (20) with c= 30, i.e., for a complex solu-
tion of (1)-(3) in the class of Li and Sinai solutions.
Fig10: Fourth case: the total energy E(t),t∈
[0,2.19 10−5]obtained with the initial data (24) with
c= 30, i.e., for a real solution of (1)-(3).
Fig11: Fourth case: the total enstrophy S(t),t∈
[0,2.19 10−5]obtained with the initial data (24) with
c= 30, i.e., for a real solution of (1)-(3).
t*1.562 108=300
t*1.562 108=700
t*1.562 108=900
-1000 -500 500 1000
k3
5.0 ×106
1.0 ×107
1.5 ×107
2.0 ×107
S3(k3,t)
Fig12: Fourth case: the marginal density for the
enstrophy along the k3axis in k-space S3(k3, tn),
k3∈[−1200,1200] for three time-iterates near
the maximum of the total enstrophy, i.e. n=
300,700,900 obtained with the initial data (24) with
c= 30, i.e., for a real solution of (1)-(3).
total energy E(t),t∈[0,2.19 10−5], Figure 5 shows
the plot of the total enstrophy S(t),t∈[0,2.19 10−5],
and Figure 6 shows the plot of the marginal density of
the enstrophy along the k3axis in k-space S3(k3, tn),
k3∈[−1200,1200] for three time-iterates, i.e. n=
500,1100,1400.
The third example considers the initial data (20)
in section 4.1, but with c= 30. The numerical re-
sults are computed by using the same discretization
parameters of the first example. The results are re-
ported in Figure 7, Figure 8, Figure 9; these figures
are similar to the ones already described in the first
example. However, the third example shows a more
rapid blowup effect than the first one.
The fourth example considers the antisymmetric
initial data (24), but with c= 30. The numerical re-
sults are computed by using the same discretization
parameters of the second example. The results are re-
ported in Figure 10, Figure 11, Figure 12; these fig-
ures are similar to the ones already described in the
second example. Also for the real solution, the case
c= 30 shows a higher increase in the total enstrophy
than the case c= 20.
For the second example, Figure 13 shows a vol-
ume plot of the local energy kuk2/2 obtained by the
inverse Fourier transform of the solution vof (13)
at the final time iterate n= 1400. We note that,
with the above grid in the k-space, uis defined for
x∈[−π, π]×[−π, π]×[−π, π]. We reported only
the central part [−50,50]×[−50,50]×[−100,100] of
the grid, which is the most significant, to highlight the
structure of the solution: Figure 13 shows ku(x, t)k2,
x∈[−0.89,0.89] ×[−0.89,0.89] ×[−0.26,0.26],
t= 2.19 10−5.
A statistical analysis of the properties of the solu-
tions in the Li-SInai class would be very interesting,
unfortunately this study has a prohibitive computa-
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
253
Volume 19, 2024
Figure 13: Second example: the local energy kuk2/2,
in the subgrid [−50,50] ×[−50,50] ×[−100,100], at
time iterate n= 1400, obtained with the initial data
(24), i.e., for a real solution of (1)-(3).
tional cost since the numerical simulation of each case
needs several tens of thousands of computation hours
in tier0 supercomputers.
Also for the real case, in spite of the fact that there
is no blow-up, the results are reliable and quite inter-
esting. They provide evidence of a large increase in
the total enstrophy that is related to a concentration of
energy at some points along the k3-axis.
Regarding the proposed numerical scheme, we can
observe by a close analysis of the plot in Figure 6
some spurious effects due to the boundary of the com-
putational domain. They are however quite small if
we consider that the values at the boundary are ob-
tained by multiplying the values of vby kkk2, and
summing over k1and k2. As the contribution of the
boundary grows with time, the problem deserves fur-
ther analysis to reduce as much as possible the effects.
The numerical results shown in sections 4.1, 4.2
are obtained by running the implementation code for
Algorithm 1 under the system Joliot Curie – KNL, ar-
chitecture BULL Sequana X1000, made available by
the GENCI, [29], within the PRACE Project Access -
Call 23, [30].
5 Conclusions
Our paper considers the incompressible Navier-
Stokes equations in an equivalent formulation of an
integral equation in Fourier space. An approxima-
tion scheme for this integral equation is proposed and
some corresponding numerical results are presented.
The numerical simulations show that the approxi-
mation scheme is able to compute accurately the solu-
tion of the integral equation arising from the Navier-
Stokes equations, even up to the blow-up time for
some complex solutions, although, as expected, the
accuracy is reduced when the solution support reaches
the boundary of the computational domain. The re-
sults of the simulations also show interesting prop-
erties of the singular complex solutions, and, for the
related real solutions provide evidence on the concen-
tration of high velocity regions and on the increase
of the total enstrophy. This study has considered the
Li-Sinai solutions and analogous real solutions ob-
tained by using antisymmetric initial data. As the
initial data for the real solutions are axial symmetric
with no swirl, it is well-known that the (real) solu-
tions are regular for all times, [31]. The present com-
putational framework can however be a useful tool to
extend the Li and Sinai approach to a wider class of
Navier-Stokes solutions. So, future studies have to
extend the class of Li-Sinai solutions to initial data
with non-zero swirl and generalize the previous re-
sults to this new class of solutions. Another interest-
ing study is the organization of the presented compu-
tation tool in a computational fluid dynamc software,
where the eventual boundary conditions can be eas-
ily treaded through an immersed boundary approach,
[32].
Acknowledgment:
C. Boldrighini and A. Pellegrinotti are members of
the Gruppo Nazionale Fisica Matematica-Istituto
Nazionale di Alta Matematica (GNFM-INdAM).
Pierluigi Maponi is a member of Gruppo Nazionale
Calcolo Scientifico-Istituto Nazionale di Alta
Matematica (GNCS-INdAM).
The computer simulations were performed at the
Joliot Curie – KNL, architecture BULL Sequana
X1000, made available by the GENCI, [29], within
the PRACE Project Access - Call 23, [30].
References:
[1] Auton, T.R., The lift force on a spherical body in
a rotational flow, Journal of fluid mechanics Vol.
183, 1987, pp. 199–218.
[2] Vegad, G.D. and Jana, A.K., Experimental and
Computational Fluid Dynamics-Based
Simulation of Oil-in-Water Emulsion Flow
through a Pipeline. Chemical Engineering &
Technolog, Vol. 46, 2023, pp. 1476-1484.
[3] Krishnamurti, T.N., ”Numerical Weather
Prediction”. Annual Review of Fluid Mechanics,
Vol. 27, 1995, pp. 95–225.
[4] Wang, Z. J., High-order computational fluid
dynamics tools for aircraft design, Philosophical
transactions of the Royal Society of London.
Series A: Mathematical, physical, and
engineering sciences, Vol.372, 2022, p.1-18.
[5] Denton, J.D., Dawes, W.N., Computational fluid
dynamics for turbomachinery design,
Proceedings of the Institution of Mechanical
Engineers, Part C: Journal of Mechanical
Engineering Science, Vol.213, 1998, pp.107-124.
[6] Giacomini, J., Invernizzi, M.C., Maponi, P.,
Verdoya, M., Testing a model of flow and heat
transfer for u-shaped geothermal exchangers.
Advances in Modelling and Analysis A, Vol. 55,
2018, pp. 151–157.
[7] Egidi, N., Giacomini, J., Maponi, P., Inverse heat
conduction to model and optimise a geothermal
field. Journal of Computational and Applied
Mathematics, Vol. 423, 2023, 114957.
[8] Giacomini, J., Khamitova, G., Maponi, P.,
Vittori, S., Fioretti, L., Water flow and transport
in porous media for in-silico espresso coffee.
International Journal of Multiphase Flow, Vol.
126, 2020, 103252.
[9] Egidi, N., Giacomini, J., Maponi, P., Perticarini,
A., Cognigni, L., Fioretti, L., An
advection–diffusion–reaction model for coffee
percolation. Computational and Applied
Mathematics, Vol. 41, 2022, 229.
[10] Angeloni, S., Giacomini, J., Maponi, P.,
Perticarini, A., Vittori, S., Cognigni, L., Fioretti,
L., Computer Percolation Models for Espresso
Coffee: State of the Art, Results and Future
Perspectives. Applied Sciences, Vol. 13, 2023,
2688.
[11] Giacomini, J., Maponi, P., Perticarini, A.,
CMMSE: a reduced percolation model for
espresso coffee. Journal of Mathematical
Chemistry, Vol. 61, 2023, pp. 520–538.
[12] Temam, R., Navier-Stokes Equations, North
Holland, 1979.
[13] Clay Mathematics Institute, The Millennium
Prize Problems, accessed 22 July 2024,
https://www.claymath.org/millennium-problems/
[14] Leray, J., Sur le mouvement d’un liquide
visqueux emplissant l’espace. Acta Math, Vol.
63, 1934, pp. 193-248.
[15] Seregin, G., A Certain Necessary Condition of
Potential Blow up for Navier-Stokes Equations.
Commun. Math. Phys., Vol. 312, 2012, pp.
833-845.
[16] T.Y. Hou, Y. Wang. 2024, Blowup Analysis for
a Quasi-exact 1D Model of 3D Euler and
Navier–Stokes, Nonlinearity, Vol. 37, 035001.
[17] T.Y. Hou, Potentially Singular Behavior of the
3D Navier–Stokes Equations, Found Comput
Math Vol. 23, 2023, pp. 2251–2299.
[18] Okita, M. , On the blow-up criterion for the
Navier–Stokes equations with critical time order,
Journal of Differential Equations, Vol. 349,
2023, pp. 269-283,
[19] Tao, T., Quantitative bounds for critically
bounded solutions to the Navier–Stokes
equations, in Nine mathematical challenges: an
elucidation, edited by A. Kechris et al., Proc.
Sympos. Pure Math., Vol. 104, Amer. Math.
Soc., Providence, RI, 2021, pp. 149–193.
[20] Wang, N., Hu, Y., Blowup of solutions for
compressible Navier–Stokes equations with
revised Maxwell’s law, Applied Mathematics
Letters, Vol. 103, 2020, 106221.
[21] T. Barker, 2023, ’Localized quantitative
estimates and potential blow-up rates for
Navier-Stokes equations’, Siam Journal on
Mathematical Analysis, vol. 55, no. 5, pp.
5221-5259. https://doi.org/10.1137/22M1527179
[22] Li, D., Sinai YA. G., Blowups of complex
solutions of the 3D Navier-Stokes system and
renormalization group method. J. Eur. Math.
Soc., Vol. 10, 2008, pp. 267-313.
[23] Boldrighini, C., Li, D., Sinai Ya., G., Complex
singular solutions of the 3-d Navier-Stokes
equations and related real solutions. Journal of
Statistical Physics, Vol. 167, 2017, pp. 1-13.
[24] Boldrighini, C., Frigio, S., Maponi, P., On the
blow-up of some complex solutions of the 3D
Navier-Stokes Equations: theoretical predictions
and computer simulation. IMA Journal of
Mathematical Physics, Vol. 83, 2017, pp.
697-714.
[25] Boldrighini, C., Frigio, S., Maponi, P.,
Pellegrinotti, A., Sinai, Y.G., An Antisymmetric
Solution of the 3D Incompressible
Navier–Stokes Equations with “Tornado-Like”
Behavior. Journal of Experimental and
Theoretical Physics, Vol. 131, 2020, pp.
356–360.
[26] C.T.H. Baker, The Numerical Treatment of
Integral Equations. Oxford Univ. Press, Oxford,
1977.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
255
Volume 19, 2024
[27] MPI Forum, MPI Forum, accessed 22 July
2024, https://www.mpi-forum.org/
[28] Li, N., Laizet, S., 2DECOMP&FFT – A highly
scalable 2D decomposition library and FFT
interface, Cray User Group 2010 conference,
Edinburgh, 2010.
[29] GENCI, GENCI, HPC at the zervice of
knowledge, accessed 22 July 2024,
https://www.genci.fr/en
[30] PRACE, PRACE, Partnership for Advanced
Computing in Europe, accessed 22 July 2024,
https://prace-ri.eu/
[31] Lei, Z., Zang, Q. Criticality of the axially
symmetric Navier-Stokes equations, Pacific
Journal of Mathematics, Vol. 289, 2017, pp.
169-187.
[32] Verzicco, R., Immersed Boundary Methods:
Historical Perspective and Future Outlook,
Annual Review of Fluid Mechanics, Vol. 55,
2023, pp. 129-155.
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The authors equally contributed in the present re-
search, at all stages from the formulation of the prob-
lem to the final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
No funding was received for conducting this study.
Conflicts of Interest
The authors have no conflicts of interest to
declare that are relevant to the content of this
article.
Creative Commons Attribution License 4.0
(Attribution 4.0 International , CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2024.19.23
C. Boldrighini, S. Frigio, P. Maponi, A. Pellegrinotti
E-ISSN: 2224-347X
256
Volume 19, 2024