Original article

Is a direct numerical simulation (DNS) of Navier-Stokes equations with small enough grid spacing and time-step definitely reliable/correct?

  • Shijie Qin b ,
  • Yu Yang b ,
  • Yongxiang Huang c ,
  • Xinyu Mei d ,
  • Lipo Wang d ,
  • Shijun Liao , a, b, *
Expand
  • a State Key Laboratory of Ocean Engineering, Shanghai 200240, China
  • b School of Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
  • c State Key Laboratory of Marine Environmental Science, College of Ocean and Earth Sciences, Xiamen University, Xiamen 361102, China
  • d UM-SJTU Joint Institute, Shanghai Jiao Tong University, Shanghai 200240, China
*E-mail address: (S. Liao).

Received date: 2024-04-11

  Revised date: 2024-04-27

  Accepted date: 2024-04-29

  Online published: 2024-05-14

Abstract

Turbulence is strongly associated with the vast majority of fluid flows in nature and industry. Traditionally, results given by the direct numerical simulation (DNS) of Navier-Stokes (NS) equations that relate to a famous millennium problem are widely regarded as ‘reliable' benchmark solutions of turbulence, as long as grid spacing is fine enough (i.e. less than the minimum Kolmogorov scale) and time-step is small enough, say, satisfying the Courant-Friedrichs-Lewy condition (Courant number < 1). Is this really true? In this paper a two-dimensional sustained turbulent Kolmogorov flow driven by an external body force governed by the NS equations under an initial condition with a spatial symmetry is investigated numerically by the two numerical methods with detailed comparisons: one is the traditional DNS, the other is the ‘clean numerical simulation' (CNS). In theory, the exact solution must have a kind of spatial symmetry since its initial condition is spatially symmetric. However, it is found that numerical noises of the DNS are quickly enlarged to the same level as the ‘true' physical solution, which finally destroy the spatial symmetry of the flow field. In other words, the DNS results of the turbulent Kolmogorov flow governed by the NS equations are badly polluted mostly. On the contrary, the numerical noise of the CNS is much smaller than the ‘true' physical solution of turbulence in a long enough interval of time so that the CNS result is very close to the ‘true' physical solution and thus can remain symmetric, which can be used as a benchmark solution for comparison. Besides, it is found that numerical noise as a kind of artificial tiny disturbances can lead to huge deviations at large scale on the two-dimensional Kolmogorov turbulence governed by the NS equations, not only quantitatively (even in statistics) but also qualitatively (such as spatial symmetry of flow). This highly suggests that fine enough spatial grid spacing with small enough time-step alone could not guarantee the validity of the DNS of the NS equations: it is only a necessary condition but not sufficient. All of these findings might challenge some of our general beliefs in turbulence: for example, it might be wrong in physics to neglect the influences of small disturbances to NS equations. Our results suggest that, from physical point of view, it should be better to use the Landau-Lifshitz-Navier-Stokes (LLNS) equations, which consider the influence of unavoidable thermal fluctuations, instead of the NS equations, to model turbulent flows.

Highlights

● A two-dimensional Kolmogorov flow is numerically solved by means of the direct numerical simulation (DNS) and clean numerical simulation (CNS), respectively.

● It is found that tiny numerical noises of the DNS result are quickly enlarged to a macroscopic level so that the DNS results are quickly polluted badly.

● Detailed comparisons between the CNS and DNS reveal that artificial numerical noises lead to large deviations of the turbulent flow even in long-term statistics.

Cite this article

Shijie Qin , Yu Yang , Yongxiang Huang , Xinyu Mei , Lipo Wang , Shijun Liao . Is a direct numerical simulation (DNS) of Navier-Stokes equations with small enough grid spacing and time-step definitely reliable/correct?[J]. Journal of Ocean Engineering and Science, 2024 , 9(3) : 293 -310 . DOI: 10.1016/j.joes.2024.04.002

1. Introduction

Turbulence is strongly associated with the vast majority of fluid flows in nature and industry [1], [2], [3]. Traditionally, results given by the direct numerical simulation (DNS) of Navier-Stokes equations under various initial/boundary conditions are widely regarded as 'reliable' benchmark solutions of turbulence, as long as spatial grid spacing is fine enough (i.e. less than the minimum Kolmogorov scale) [4] and time-step is small enough, say, satisfying the Courant-Friedrichs-Lewy condition (Courant number < 1) [5]. Is this definitely true?
It is widely believed that turbulence flow, as a typical spatio-temporal chaotic system [1], [3], [4], has the famous 'butterfly effect' [6], [7]: microscopic disturbances in turbulent flows, caused by either thermal fluctuations or environmental noises that are small-scale but unavoidable in practice, might increase exponentially (and quickly) to macroscopic levels [8], [9], [10], [11]. Recently, Qin and Liao [12] provided a rigorous evidence that numerical noise as a kind of tiny artificial stochastic disturbances can indeed increase quickly to the same order of magnitude of the 'true' physical solution of a two-dimensional sustained turbulent Rayleigh-Bénard convection, and in addition can have quantitative and qualitative large-scale influences in statistics, and even might lead to different types of flow.
Note that numerical noises, i.e. truncation error and round-off error, always exist and are practically unavoidable in each result given by the 'direct numerical simulation' (DNS). Unfortunately, as mentioned above, due to the chaotic property of turbulence [13], [14], the tiny numerical noises as a kind of micro-level disturbances increase exponentially (and thus quickly) to a macroscopic level, therefore each DNS result of turbulent flow is a kind of mixture of 'true' physical solution (denoted by p ) and 'false' numerical noise (denoted by δ ) that are mostly at the same order of magnitude, i.e. pδ , as currently illustrated by Qin and Liao [12] using a two-dimensional Rayleigh-Bénard turbulent flow as an example by means of the so-called 'clean numerical simulation' (CNS) [12], [15], [16], [17], [18], [19]. For investigations on turbulence, statistic properties are more important than spatio-temporal trajectory of numerical simulation. It is a pity that, traditionally, the 'true' solution p is unknown so that one had to use the the badly polluted simulation, i.e. the above-mentioned mixture p+δ , to gain statistic results p+δ , where denotes an operator of statistics. Unfortunately, this is based on such a 'hypothesis' that p+δ=p , say, the statistics is stable under disturbances, where p denotes the 'true' physical solution. But, as illustrated by Liao [16], this hypothesis is not definitely true for many chaotic systems [18], [19] and even some turbulent flows: the DNS results of a kind of two-dimensional Rayleigh-Bénard turbulent flow are quite different from these given by CNS not only in statistics but also even in type of flow [12].
Here, let us briefly describe the basic idea of the CNS. In order to obtain a 'convergent' (or reproducible) numerical simulation (i.e. a result very close to 'true' solution p ) of chaotic dynamical systems [2], [6], [20], [21], [22], Liao [15] proposed a new numerical strategy, namely the 'clean numerical simulation' (CNS) [12], [15], [16], [17], [18], [19], [23], [24], to reduce the background numerical noise (i.e. the truncation error and round-off error) so greatly that numerical noise is negligible (say, the 'false' numerical noise δ is much smaller than the 'true' physical solution p , i.e. |δ|?|p| ) during a finite but long enough interval of time t[0,Tc] , where Tc is the so-called 'critical predictable time'. In the frame of the CNS, temporal and spatial truncation errors can be decreased to a required small level by means of the Taylor expansion with a high enough order in the temporal dimension and a fine enough discretization in the spatial dimension (such as the high-order spatial Fourier expansion), respectively. Besides, all of the physical/numerical variables and parameters must be in the multiple precision (MP) [25] with a large enough number of significant digits so that the round-off error can be decreased to a required small level. Moreover, an additional numerical simulation with the even smaller numerical noise of the same chaotic system is required to determine the value of Tc via a comparison with the previous CNS result so that the corresponding numerical noise of the simulation could be negligible and thus this computer-generated simulation is convergent/reproducible within the whole spatial domain in the temporal interval [0,Tc] . In this way, unlike other traditional numerical algorithms using double precision, the CNS can give reliable numerical simulation of a chaotic system in a finite but long enough interval of time, say, the CNS result is so close to the 'true' physical solution p that one can gain its statistics p accurately. Here, it should be emphasized that the so-called 'critical predictable time' denoted by Tc is a very important concept in the frame of the CNS: unlike DNS results whose interval of time can be as long as one would like without considering the influences of the exponential increase of numerical noises, CNS can give a reliable result only in a finite but long enough interval of time, i.e. t[0,Tc] . So, for the first time, the CNS provides us a tool to carefully check the hypothesis p+δ=p for chaotic systems and turbulent flows, where p+δ on the left hand side denotes numerical simulations given by a traditional numerical method using double precision (such as DNS), and p on the right-hand side can be accurately gained by means of the CNS, respectively.
Up to now, the CNS has been successfully applied to many chaotic dynamical systems, such as Lorenz system [15], [26], the famous three-body problem [24], [27], [28], [29], [30], Hénon-Heiles chaotic system [23], hyperchaotic R?ssler system [31], [32], a chaotic free-fall disk [33], Arnold-Beltrami-Childress (ABC) flow [34], the complex Ginzburg-Landau equation (CGLE) [18], the damped driven sine-Gordon equation (SGE) [19], the two-dimensional turbulent Rayleigh-Bénard convection [12], and so on. It was found that the hypothesis p+δ=p indeed holds for many chaotic systems and turbulent flows, say, their statistics are stable, which are called 'normal-chaos'. But, it was found that the hypothesis p+δ=p does not hold for some chaotic systems and turbulent flows, say, their statistics are unstable to small disturbances, which are called 'ultra-chaos' [31], [34], [35], [36]. The ultra-chaos as a brand-new concept might greatly deepen and enrich our understandings about chaos as well as turbulence.
In this paper, let us further consider a two-dimensional turburlent flow of viscous fluid, driven by an external body forcing that is stationary, monochromatic, and sinusoidally varying in space. This kind of flow was first introduced by Kolmogorov in 1959 [37], called 'Kolmogorov flow' and usually regarded as a mathematically and experimentally tractable flow model whose instability [38], [39], [40], [41], global stability [42], chaotic property [13], [14], and turbulent state [43], [44], [45] have been investigated in details. The Kolmogorov flow can be practically realized in the laboratory by means of the magnetohydrodynamics (MHD) [46], [47], [48] or in a soap film [49]. Furthermore, due to some additional physical factors such as the Coriolis term (i.e. the effect of the planetary rotation) [50], [51], the bottom friction [52], [53], [54], and the stratification [41], [55], Kolmogorov flows have also been widely studied in geophysical fluid dynamics.
In this investigation, we compare the CNS result of the two-dimensional Kolmogorov turbulence with that given by a traditional DNS using double precision. In this way, we can check the hypothesis p+δ=p in details, where the left-hand side is given by the DNS, and the right-hand side is given by the CNS, respectively. The mathematical model is described in Section 2, the detailed comparisons between CNS and DNS results are given in Section 3, and the discussions and conclusions are given in Section 4, respectively. As reported in Section 3, numerical noises as a kind of small-scale disturbances have huge influences on large-scale properties of the two-dimensional Kolmogorov turbulence not only in spatio-temporal trajectories and its spatial symmetry but also in various statistics. Especially, it is found that the hypothesis p+δ=p does not stand up for the two-dimensional Kolmogorov turbulence under consideration, say, the fine enough spatial grid spacing with small enough time-step cannot guarantee the validity of the DNS: it is only a necessary condition but not sufficient.

2. Mathematical model and numerical strategy

2.1. Governing equation and initial/boundary conditions

Consider an incompressible flow in a domain [0,L]×[0,L] under the so-called 'Kolmogorov forcing', which is stationary, monochromatic, and sinusoidally varying in space, with an integer nK describing the forcing scale and χ representing the corresponding forcing amplitude per unit mass of fluid [45]. Using the length scale L/2π and the time scale L/2πχ , the non-dimensional governing equation of this two-dimensional Kolmogorov flow in the form of stream function reads
??t(?2ψ)+?(ψ,?2ψ)?(x,y)?1Re?4ψ+nKcos(nKy)=0,
where x and y are horizontal and vertical coordinates with x,y[0,2π] , t denotes the time, the stream function ψ is defined by
u=??ψ?y,v=?ψ?x,
u and v represent horizontal and vertical velocities, Re is the Reynolds number defined by
Re=χν(L2π)32,
ν denotes the kinematic viscosity, ?2 is the Laplace operator with the definition ?4=?2?2 , and
?(a,b)?(x,y)=?a?x?b?y??b?x?a?y
is the Jacobi operator, respectively. Note that the stream-function ψ satisfies the periodic condition on the boundary
ψ(0,y,t)=ψ(2π,y,t),ψ(x,0,t)=ψ(x,2π,t).
Following Chandler and Kerswell [45], we choose the physical parameters in (1) as follows:
nK=4,Re=40,
with the initial condition (this is different from Chandler and Kerswell [45])
ψ(x,y,0)=?12{cos(x+y)+cos(x?y)},
corresponding to a sustained turbulent Kolmogorov flow.
Note that there exists a kind of spatial symmetry in the initial condition (7). According to the governing equation (1), the periodic boundary condition (5) and the initial condition (7) with the spatial symmetry, it is straight forward that the solution of the corresponding NS equations should be in the form of a series
ψ(x,y,t)=m=0+n=1+{am,n(t)cos(mx+ny)+bm,n(t)cos(mx?ny)},
where am,n(t) and bm,n(t) are unknown coefficients dependent upon the time t . Therefore, the exact solution of the Navier-Stokes equations have the spatial symmetry
ψ(x,y,t)=ψ(?x,?y,t)=ψ(2π?x,2π?y,t),
although the coefficients am,n(t) and bm,n(t) in (8) are unknown. Thus, the loss of such kind of spatial symmetry clearly indicates the large deviation from the true solution of the Navier-Stokes equations. This is exactly the reason why we choose the periodic boundary condition (5) and the initial condition (7) with the spatial symmetry!

2.2. Strategy of CNS

First of all, let us simply describe the strategy of CNS. The CNS is used to greatly decrease the background numerical noise, i.e. the truncation and round-off errors, to such a required tiny level that the corresponding 'false' numerical noise of the simulation is much smaller than, and thus negligible compared with, the 'true' physical solution of the considered Kolmogorov flow in an interval of time that is long enough for statistics. In this way, a convergent (reproducible) and reliable numerical result of the two-dimensional Kolmogorov turbulence can be obtained, which is used here as the 'clean' and 'true' benchmark solution. Briefly speaking, to decrease the spatial truncation error of the simulation to a small enough level, we discretize the spatial domain of the flow field by a uniform mesh Nx×Ny=256×256 , and adopt the Fourier spectral method for spatial approximation with the 3/2 rule for dealiasing. The corresponding spatial resolution is high enough for the Kolmogorov flow under consideration: the grid spacing is less than the minimum instantaneous Kolmogorov scale [4], as shown later in Section 3.5. Besides, to decrease the temporal truncation error to a required small enough level, we use the 60th-order (i.e. M=60 ) Taylor expansion with a time-step Δt=1×10?3 . Furthermore, in order to decrease the round-off error to a required small enough level, we use multiple precision with 200 significant digits (i.e. Ns=200 ) for all physical/numerical variables and parameters. In addition, the self-adaptive CNS strategy [32] is adopted to dramatically reduce the calculation amount. Similarly, we run another CNS using the same Fourier spectral method with the same uniform spatial mesh but having the even smaller temporal truncation error and round-off error by means of a higher-order (i.e. M=62 ) Taylor expansion with the same time-step (Δt=1×10?3 ) and the higher multiple precision with more significant digits (i.e. Ns=205 ). Comparing these two CNS results by means of the so-called 'spectrum-deviation' δs(t) whose definition is given in Section 2.3, it is found that δs<10?19 is satisfied throughout the whole simulation 0t1500 , say, these two CNS results have no distinct difference in an interval of time 0t1500 that is long enough for calculating statistics. This verifies the convergence and reliability of our CNS result in t[0,1500] given by means of M=60 , Δt=10?3 and Ns=200 , which can be regarded as a 'clean' benchmark result, since it is very close to the 'true' physical solution of the two-dimensional Kolmogorov turbulence under consideration. For details about the CNS algorithm of the Kolmogorov turbulence, please refer to Appendix A.
On the other hand, in the case of the same initial condition (7) and the same physical parameters (6), the same governing equation (1) is solved numerically in t[0,1500] by a traditional DNS algorithm, i.e. the fourth-order Runge-Kutta's method with the double precision using the time-step Δt=10?4 , whose numerical noise increases exponentially (and quickly) up to the same order of magnitude of the 'true' physical solution and thus can not be negligible thereafter. Therefore, the corresponding DNS result becomes a mixture of the 'true' physical solution and the 'false' numerical noise, which are mostly at the same order of magnitude. Then, comparing this DNS result with the CNS benchmark solution, we can investigate the influence of numerical noises as tiny artificial stochastic disturbances on the two-dimensional Kolmogorov turbulence in details.
For both of DNS and CNS, the grid spacing is fine enough, i.e. satisfying the spacing criterion [4], and the time-steps are small enough, i.e. satisfying the Courant-Friedrichs-Lewy condition [5], as shown in Section 3.5 in details.

2.3. Some measures of the flow

For the sake of simplicity, the definitions of some statistic operators are briefly described below. The spatial average is defined by
A=14π202π02πdxdy
and the spatiotemporal average (along the x direction) is defined by
x,t=12π(T2?T1)02πT1T2dxdt,
respectively, where T1=200 and T2=1500 are chosen in Section 3 corresponding to a long enough temporal interval, as well as a relatively stable state of the flow. Besides, the spatiotemporal average (over the whole field) is defined by
=14π2(T2?T1)02π02πT1T2dxdydt,
where different values of T1 and T2 correspond to different integral intervals of time.
For the two-dimensional Kolmogorov turbulence under consideration, vorticity is expanded as the Fourier series
ω(x,y,t)=?2ψ(x,y,t)=m=??Nx/3??Nx/3?n=??Ny/3??Ny/3?Ωm,n(t)exp(imx)exp(iny),
where m , n are integers, ?? stands for a floor function, i=?1 denotes the imaginary unit, and for dealiasing Ωm,n=0 is imposed for wavenumbers outside the above domain , respectively. Note that, when ω is a real number, Ω?m,?n=Ωm,n* must be satisfied, where Ωm,n* is the conjugate of Ωm,n . The enstrophy spectrum is defined by
Bl(t)=l?1/2m2+n2<l+1/2|Ωm,n(t)|2,
where the wave number l is a non-negative integer. In addition, for two different numerical results that give different enstrophy spectrum Bl(t) and Bl(t) , respectively, where Bl(t) corresponds to the simulation with smaller numerical noise, we define the so-called 'spectrum-deviation'
δs(t)=l=0|Bl(t)?Bl(t)|l=0Bl(t)
to measure the deviation of the simulation corresponding to Bl(t) , from the other corresponding to Bl(t) , at a given time t .
In Section 3, we will focus on the kinetic energy
E(x,y,t)=12[u2(x,y,t)+v2(x,y,t)]
and the kinetic energy dissipation rate
D(x,y,t)=12Reij[?iuj(x,y,t)+?jui(x,y,t)]2,
where i,j=1,2 , u1=u , u2=v , ?1=?/?x , and ?2=?/?y , respectively.

3. Influence of numerical noises as tiny disturbances

As mentioned in Section 2.2, comparing the DNS result, which is a mixture of the 'true' physical solution and the 'false' numerical noise that are mostly at the same order of magnitude, with the CNS benchmark solution, whose 'false' numerical noises are negligible and thus which is very close to its 'true' physical solution in a finite but long enough interval of time, we can investigate in details the influence of numerical noise as a kind of tiny artificial stochastic disturbances on the two-dimensional Kolmogorov turbulence under consideration. Some results of detailed comparisons are shown in this section.

3.1. Trajectory, spatial symmetry and statistics

As shown in Fig. 1 and 2, the DNS results of vorticity ω defined by (13) at different time t are compared with the corresponding 'clean' benchmark solution given by the CNS. Note that the flow states given by DNS and CNS are completely the same at the beginning times such as t=90 , as shown in Fig. 1(a) and (b), when the numerical noise of the DNS result has not yet increased to the same order of magnitude of the 'true' physical solution so that the macroscopic deviation between DNS and CNS cannot be observed by eyes. It should be emphasized that, for both results given by the DNS and the CNS up to t=90 , there exists a kind of spatial symmetry, as shown in Fig. 1(a) and (b). So, the DNS result in t[0,90] should be close enough to the 'true' physical solution. Unfortunately, the corresponding interval of time is too short for calculations of statistics. As the time increases, for instances at t=110 and t=120 , the deviation between the DNS result from the CNS benchmark solution becomes more and more obvious, as shown in Fig. 1(e)-(h). To clearly indicate this deviation, let us consider the normalized absolute error of the vorticity between DNS result and the CNS benchmark solution, defined by
δω=02π02π|ωDNS?ωCNS|dxdy02π02π|ωCNS|dxdy,
where ωDNS and ωCNS denote the vorticity given by the DNS and CNS, respectively. It is found that the normalized absolute errors of the vorticity, denoted by δω , are 0.006 at t=90 , 0.06 at t=100 , 0.49 at t=110 and 1.12 at t=120 , respectively. Note that, when t120 , the vorticity field given by the DNS completely losses the spatial symmetry, but the vorticity field given by the CNS benchmark solution always keeps such kind of spatial symmetry even at t=1500 , as shown in Fig. 2. Note that the 'false' numerical noise of the CNS benchmark result is much smaller than its 'true' physical solution and thus is negligible compared to the 'true' physical solution in the whole interval of time t[0,1500] . However, when t120 , the 'false' numerical noise of the DNS result is at the same order of magnitude as the 'true' physical solution, as shown in Fig. 3, say, thereafter the DNS result is badly polluted, which naturally leads to the loss of the spatial symmetry and a large-scale illustrative deviation in the flow field from the 'true' physical solution! The quantitative comparisons of the field geometrical structures will be given in Section 3.2.
Fig.1 Comparisons of the vorticity field ω of the two-dimensional Kolmogorov turbulence governed by (1) in the case of nK=4 and Re=40 with the same initial condition (7) given by the CNS (left) and DNS (right), respectively, at the different times. (a)-(b):t=90 ; (c)-(d): t=100 ; (e)-(f): t=110 ; (g)-(h): t=120 . For details, please see the corresponding video.
Fig.2 Comparisons of the vorticity field ω of the two-dimensional Kolmogorov turbulence governed by (1) in the case of nK=4 and Re=40 with the same initial condition (7) given by the CNS (left) and DNS (right), respectively, at the different times. (a)-(b): t=500 ; (c)-(d): t=1500 . For details, please see the corresponding video.
Fig.3 The normalized absolute error δω , defined by (18), and the spectrum-deviation δs , defined by (15), of the vorticity field given by the DNS of the two-dimensional Kolmogorov turbulence governed by (1) in the case of nK=4 and Re=40 with the initial condition (7), compared to the CNS benchmark solution. Red solid line: the normalized absolute error δω ; blue dash line: the spectrum-deviation δs .
As mentioned in Section 2.2, the exact solution of the Navier-Stokes equations considered in this paper have the spatial symmetry (9). Note that the result given by the CNS indeed has this kind of spatial symmetry in the whole interval of time t[0,1500] . But, unfortunately, the result given by the DNS loses such kind of spatial symmetry when t>120 , as shown in Figs. 1 and 2, because the numerical noises, which have no spatial symmetry, have been enlarged to the same order of magnitude as the 'true' solution and thus finally destroy this kind of spatial symmetry of the flow field. This is a clear and rigorous evidence that the numerical results given by the DNS are indeed “badly polluted” when t>120 .
How about the influences of the numerical noises as a kind of artificial tiny disturbances in statistics? First of all, let us consider the spatially averaged kinetic energy dissipation rate DA , where A is an operator of statistics defined by (10), and D denotes the kinetic energy dissipation rate defined by (17), respectively. As shown in Fig. 4, when t>150 , DA given by the CNS (with a approximate average 0.25) is obviously higher than that given by the DNS (with a approximate average 0.1), where the deviation of nearly 2.5 times in magnitude clearly reveals the large differences. Note that our DNS result is in accord with the DNS ones given by Chandler and Kerswell [45], which is also a mixture of the 'true' physical solution and the 'false' numerical noise that are mostly at the same order of magnitude.
Fig.4 Comparison of the spatially averaged kinetic energy dissipation rate DA of the two-dimensional Kolmogorov turbulence governed by (1) in the case of nK=4 and Re=40 with the initial condition (7) given by the CNS benchmark result (solid line in red), the DNS result (solid line in black), the previous study by Chandler and Kerswell [45] (solid line in blue).
Let us further consider the Fourier coefficient Ωm,n of the vorticity field ω defined by (13) at the time t=500 , corresponding to Fig. 2 (a) and (b). As shown in Fig. 5, when m is an odd number such as m=1 or m=3 , Re(Ωm,n ) given by the CNS benchmark solution is always zero when n is an even number, but Re(Ωm,n ) given by the DNS is not zero for all n , where Re( ) is an operator that takes the real part of a complex number. On the other hand, when m is an even number such as m=0 or m=2 , Re(Ωm,n ) given by the CNS is zero when n is an odd number, but Re(Ωm,n ) given by the DNS is not zero nearly for all n (except Ω0,0=0 ), as shown in Fig. 6. Besides, it is found that, when Re(Ωm,n ) given by the CNS are not zero, their values are usually an order of magnitude larger than those given by the DNS. In addition, the imaginary part Im(Ωm,n ) given by the CNS is always zero for all values of m and n , but this is not found in the DNS result, where Im( ) is an operator that takes the imaginary part of a complex number. These are mainly because the 'false' numerical noises of the DNS quickly increase to the same order of magnitude of the 'true' physical solution of the vorticity field, which certainly leads to the large-scale huge deviation from the CNS benchmark solution! All of these can clearly explain why the DNS result losses the spatial symmetry when t>120 .
Fig.5 Comparisons of (a) Re(Ω1,n) and (b) Re(Ω3,n) of the two-dimensional Kolmogorov turbulence at t=500 in the case of nK=4 and Re=40 with the the initial condition (7), where the Fourier coefficient Ωm,n of the vorticity ω defined by (13) is given either by the CNS benchmark solution (solid triangle in red) or the DNS result (odd n : empty circle in black, even n : solid circle in black). Here Re(a ) denotes the real part of the complex number a , |a| represents the absolute value of a , respectively.
Fig.6 Comparisons of (a) Re(Ω0,n) and (b) Re(Ω2,n) of the two-dimensional Kolmogorov turbulence at t=500 in the case of nK=4 and Re=40 with the initial condition (7), where Ωm,n of the vorticity ω defined by (13) is given either by the CNS benchmark solution (solid triangle in red) or the DNS result (even n : empty circle in black, odd n : solid circle in black).
Furthermore, the DNS results of the enstrophy spectrum Bl defined by (14) of the two-dimensional Kolmogorov turbulence have also huge deviation from the CNS benchmark solution, as shown in Fig. 7 and Fig. 8. The time histories of the enstrophy spectrums at some wave numbers, such as B1(t) and B50(t) as shown in Fig. 7 (a) and (b), respectively, present the apparent sensitivity to numerical noises as a kind of artificial tiny disturbances: the average value (about 1.13) of B1(t) given by the DNS is larger than that (about 0.60) given by the CNS, but B1(t) given by the CNS has larger variation range, which means that the flow given by the DNS contains more kinetic energy at the macroscopical scale and has a higher auto-correlation (that will be shown later in detail); On the contrary, the average value (about 6.0×10?12 ) of B50(t) given by the CNS are about 4 times larger than that (about 1.5×10?12 ) given by the DNS, indicating that the flow given by the CNS has more energy at the microscopical scale. These phenomena can be well explained by the obvious difference between the kinetic energy dissipation rates given by the CNS and DNS, as shown in Fig. 4. Besides, at some typical times such as t=500 and t=1500 as shown in Fig. 8 (a) and (b), respectively, as l increases, deviation between the spectrum Bl given by the DNS and CNS benchmark solution becomes larger and larger, and when l>10 the values of Bl given by the CNS are obviously higher than that given by the DNS. Note that these agree with the results in Fig. 5 and Fig. 6. It suggests that, for the relatively larger wave number (i.e. smaller turbulent structure), the flow given by the CNS should contain more kinetic energy than the DNS, which is reasonable considering the apparently higher kinetic energy dissipation rate of the CNS as shown in Fig. 4. In summary, the numerical noises as a kind of artificial tiny disturbances lead to huge deviations of the enstrophy spectrum given by the DNS from the CNS benchmark solution, on both of large and small scales, as shown in Fig. 7 and Fig. 8.
Fig.7 Comparisons of the enstrophy spectrums at certain wavenumbers, i.e. (a) B1 and (b) B50 , of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), given by the CNS (red) and the DNS (black), respectively.
Fig.8 Comparisons of the enstrophy spectrums Bl at (a) t=500 and (b) t=1500 of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), given by the CNS (triangle in red) and the DNS (circle in black), respectively.
The comparisons of the probability density function (PDF) of the kinetic energy E defined by (16) and the kinetic energy dissipation rate D defined by (17) given by the CNS and DNS are as shown in Fig. 9 (a) and (b), respectively. Note that in the whole paper x,y[0,2π] and t[200,1500] are used to calculate the PDFs. Near E1 , the DNS result has a larger probability of the distribution than the CNS benchmark solution, as shown in Fig. 9 (a). For large D , the CNS benchmark solution has a larger probability of the distribution than the DNS result, as shown in Fig. 9 (b). The comparisons of the vertical distributions of the spatio-temporal averaged kinetic energy Ex,t and the spatio-temporal averaged kinetic energy dissipation rate Dx,t given by the CNS and DNS, respectively, are as shown in Fig. 10, where x,t is an operator of statistics defined by (11). Note that the DNS result mostly has larger Ex,t than the CNS benchmark solution, but Dx,t given by the CNS is several times larger than that of the DNS. Besides, although the PDF of the horizontal velocity u(x,y,t) given by the DNS is close to that given by the CNS, the PDF of the vertical velocity v(x,y,t) given by the DNS is quite different from that given by the CNS, as shown in Fig. 11. In addition, the PDF of the vorticity ω(x,y,t) given by the CNS and DNS are also distinctly different: the vorticity ω given by the DNS is more concentrated on zero, as shown in Fig. 12.
Fig.9 Comparisons of the probability density function (PDF) of (a) the kinetic energy E(x,y,t) and (b) the kinetic energy dissipation rate D(x,y,t) of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), given by the CNS (red) and the DNS (black), respectively, where PDFs are integrated in x,y[0,2π] and t[200,1500] .
Fig.10 Comparisons of the vertical distributions of (a) the spatio-temporal averaged kinetic energy Ex,t and (b) the spatio-temporal averaged kinetic energy dissipation rate Dx,t of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), respectively, given by the CNS (red) and the DNS (black).
Fig.11 Comparisons of the PDF of (a) the horizontal velocity u(x,y,t) and (b) the vertical velocity v(x,y,t) of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), respectively, given by the CNS (red) and the DNS (black).
Fig.12 Comparison of the PDF of the vorticity ω of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), respectively, given by the CNS (red) and the DNS (black).
All of the above-mentioned comparisons between the CNS benchmark solution and the DNS result clearly indicate that tiny disturbances, resulting from the micro-level background numerical noise, can lead to huge deviations on both of large and small scales not only in spatio-temporal trajectory of the flow field but also in spatial symmetry of the vorticity field as well as statistics of many physical quantities, such as the enstrophy spectrum, the spatio-temporal average and PDFs of the velocity, vorticity, kinetic energy, and kinetic energy dissipation rate of the two-dimensional Kolmogorov turbulence considered in this paper.

3.2. Field geometrical structures

Based on the aforementioned investigation, the evolution of the flow field can also be observed through the changes in field geometrical structures. For the non-local statistics related to the topology of flow field, the main challenge lies in the quantitative (rather than illustrative) structure identification, which is relatively rare in the existing literature. To overcome this challenge, a novel attempt for the natural decomposition of scalar fields is the dissipation element (DE) analysis [56], [57], [58]. The corresponding DE structures, as a representation of the whole flow field, are space-filling and able to effectively capture and portray the flow dynamics across the entire field. Additionally, the utilization of the DE analysis approach can avoid the defect of the so-called scale-mixing effect [56], which is often encountered in traditional structure function analyses.
Specifically, starting from any spatial point in a (large enough) scalar field, along with its descending and ascending directions of the scalar gradient trajectory, one can inevitably reach a local minimum point and a local maximum point, respectively. The ensemble of spatial points whose gradient trajectories share the same pair of extremal points define a DE unit, and thus the DE structures consisting of several DE units can represent the whole flow field. Typically, the length scale of each DE unit can be parameterized with l , say, the linear distance between the two extremal points. For a deeper understanding of the properties and applications of the DE analysis, please refer to [56], [57], [58], [59].
Applying the DE analysis on both the CNS benchmark solution and the DNS result, without loss of generality, for the instantaneous vertical velocity field v(x,y) at t=500 , the difference between topological properties of these two flow fields can be revealed, as shown in Fig. 13: the local maximal and minimal points are represented by the red and blue dots respectively and the DE boundaries are shown in yellow solid lines. Note that here we choose the vertical velocity field rather than the horizontal one mainly for exhibiting an obvious deviation between the CNS benchmark solution and the DNS result, which corresponds to the results shown in Fig. 11. At this moment, the sizes of the DE structures for the DNS result are generally larger than those given by the CNS in the same fixed computational domain. Quantitatively, on the one hand, the evolution of the averaged length scale, i.e. l , of the DE units (obtained by means of each instantaneous vertical velocity field) is shown in Fig. 14(a), which indicates that the averaged length scales at different time moments given by the DNS are generally larger than those of the CNS benchmark solution after an initial period of time. On the other hand, Fig. 14(b) illustrates that the total number (denoted by Num ) of the DE units (obtained likewise via each instantaneous vertical velocity field) given by the DNS is less than that given by the CNS most of the time, which is a natural result corresponding to the difference between length scales mentioned above. Note that all of these results can be explained from a perspective of the kinetic energy dissipation rate (as shown in Fig. 4 and Fig. 9(b)): considering the definition of a DE unit, the lower kinetic energy dissipation rate given by the DNS usually corresponds to the smoother and larger DE structures.
Fig.13 Illustration of the dissipation element (DE) structures on the instantaneous vertical velocity field v(x,y) at t=500 for (a) the CNS benchmark solution or (b) the DNS result, where color represents the magnitude of v , the DE boundaries are shown in yellow solid lines, and the extremal points are marked by red (maximum) and blue (minimum) dots.
Fig.14 Comparison of the time histories of (a) averaged length scale l of the DE units or (b) total number Num of the DE units, where the corresponding DE structures are obtained by means of the instantaneous vertical velocity field given by the CNS (red) or the DNS (black).
In summary, the numerical noises as a kind of artificial tiny disturbances leads to quantitatively large deviations on the field geometrical structures of the two-dimensional Kolmogorov turbulence under consideration.

3.3. Intermittent stability

According to the time histories of EA and DA in t[700,1200] given by the CNS benchmark solution, as shown in Fig. 15(a) and (b), there exists a short period of time, for example, 870<t<970 , in which the two-dimensional Kolmogorov turbulence is in a relatively stable state, corresponding to an intermittent stability of turbulence [60], [61], [62]. Note that the vorticity fields ω given by the CNS benchmark solution at t=900 , t=920 , t=940 , and t=960 are almost unchanged, as shown in Fig. 16(a)-(d). By contrast, the corresponding vorticity fields given by the DNS has no such kind of intermittent stability of turbulence in the whole period of simulation!
Fig.15 Time histories of EA and DA of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7) given by the (a)-(b) CNS and (c)-(d) DNS in t[700,1200] , respectively, where the period of time 870<t<970 corresponds to the intermittent stability of Kolmogorov turbulence.
Fig.16 Vorticity field ω of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7) given by the CNS in the state of the intermittent stability at (a) t=900 , (b) t=920 , (c) t=940 , and (d) t=960 .
It should be emphasized that the DNS result quickly becomes a mixture of the 'false' numerical noise δ and the 'true' physical solution p , which are at the same order of magnitude mostly, i.e. δp . On the contrary, the numerical noise δ of the CNS benchmark solution is much smaller than the 'true' physical solution p , i.e. |δ|?|p| , and thus negligible, so that the CNS benchmark solution is very close to the 'true' physical solution p . As illustrated above in this section, the DNS result of the two-dimensional Kolmogorov turbulence has significant deviations from the CNS benchmark solution not only in trajectory but also in statistics of many physical variables.
Note that the intermittent stability of turbulence in 870<t<970 can be successfully found by the CNS, but is unfortunately lost by the DNS. This indicates once again that CNS should be a powerful tool to investigate turbulent flows more accurately than DNS. Hopefully, CNS as a new method could bring us some new discoveries in turbulence.

3.4. Scale-to-scale energy flux

To delve into the underlying physical mechanisms, we employ the Filter-Space-Technique (FST) to extract the scale-to-scale energy flux, denoted as Π[l] (see definition below). FST, initially developed for large eddy simulation in the 1970s [63], involves applying a low-pass filter to the velocity field. Mathematically, this is expressed as:
f[l](x,t)=G[l](x?x)f(x,t)dx,
where f=u(x,t) represents a vector of the two-dimensional velocity field, and the Gaussian kernel G[l](r)=exp(?r2/2l2) is commonly used as the filter [64], [65], [66]. For the incompressible Navier-Stokes equation, the scale-to-scale energy flux can be derived analytically as:
Π[l](x,t)=?i,j=1,2[(uiuj)[l]?ui[l]uj[l]]?ui[l]?xj
where u1=u , u2=v are velocity components and x1=x , x2=y denote spatial coordinates. This term emerges from filtering the nonlinear terms in Navier-Stokes equations, capturing the interaction between scales removed by the filter and those that remain. Crucially, the sign of Π[l] reveals the direction of energy transfer: a positive value indicates a cascade from larger scales (>l ) to smaller scales (<l ), while a negative value signifies the reverse. Thus, Π[l] provides valuable insights into both the direction and intensity of energy transfer across scales in turbulent flows.
To quantify the energy transfer across different scales, we apply FST to the velocity field in this investigation, extracting Π[l] for scales ranging from 0 to 2.45. Fig. 17(a) depicts the averaged scale-to-scale energy flux Π[l] as a function of the scale l for two distinct time periods: 50t100 and 110t250 , where the spatio-temporal average is defined by (12). For the two-dimensional Kolmogorov turbulence investigated in this paper, the consistently negative values of Π[l] across all scales confirm the presence of an inverse energy cascade that is a characteristic of two-dimensional turbulence. Notably, while the CNS curves remain unchanged between these two time periods, a marked deviation is observed between the DNS curves, which corresponds to the emergence of the above-mentioned qualitative difference of the flow field: when t120 , the DNS result is badly polluted by numerical noises, as shown in Fig. 1 and Fig. 2. Furthermore, FST reveals a significant discrepancy in Π[l] at scales around llf/2 (where lf=2π/nK=1.57 denotes the forcing scale) between the results obtained by the CNS and DNS, even in the early stages (such as 50t100 ) when the visible difference between these two cases' velocity fields is absent, which is particularly intriguing and of course caused by the numerical noise.
Fig.17 (a) Averaged scale-to-scale energy flux Π[l] and (b) the corresponding derivative dlΠ[l] versus the scale l of the two-dimensional Kolmogorov turbulence investigated in this paper, where the vertical dash line in black corresponds to the forcing scale lf=2π/nK=1.57 and the horizontal one denotes the zero value of vertical coordinate. These averaged results are integrated in t[50,100] (solid lines, i.e. '-1') and t[110,250] (dash and dot lines, i.e. '-2') by means of the CNS benchmark solution (red and blue lines) and the DNS result (green and orange lines).
Considering that the CNS has greatly reduced the background numerical noise so that the corresponding result can be regarded as a 'clean' benchmark solution, these findings, coupled with the observed inverse energy cascade, provides rigorous evidence that the deviation caused by the numerical noise of the DNS first manifest at the smallest scales and subsequently propagate upwards to larger scales along with the inverse energy cascade in this two-dimensional Kolmogorov turbulence. This amplification process eventually leads to the emergence of significant differences at larger scales as time progresses.
To delve deeper into the intricacies of the energy cascade loop, we turn our attention to the local balance of energy flux at a given scale l , as described by the following equation [1]:
dlΠ[l]=?ν(l)?Ein(l).
In this equation, dl stands for the derivation with respect to l , ?ν(l) represents the energy dissipation density function due to fluid viscosity, while Ein(l) denotes the energy injection density function resulting from external forces. It is worth noting that ?ν(l) must be greater than or equal to zero, and it approaches zero for scales significantly larger than the Kolmogorov scale (such as its minimum ηmin0.042 mentioned in Section 3.5). This allows us to simplify the equation as follows:
Ein(l)?dlΠ[l].
Fig. 17(b) presents the measured dlΠ[l] versus the scale l in this two-dimensional Kolmogorov turbulence, providing insights into energy injection and extraction processes. The results reveal that kinetic energy is injected into the system within the scale ranges of 0.1?l?0.7 for the CNS and 0.1?l?1.0 for the DNS. However, when scales exceed these thresholds (0.7 for the CNS and 1.0 for the DNS), kinetic energy is extracted by external forcing. This extraction is attributed to the alignment of the velocity field and the external force with two opposite phases at large scale. It is worth noting that this pattern of energy injection and extraction has also been observed in other two- and quasi-two-dimensional systems of turbulence, which is a topic that we plan to explore in greater detail in the future.
The inverse energy cascade has been observed in various (quasi) 2D flow systems, including Kraichnan turbulence [65], [67], bacterial turbulence [68], oceanic flows [69], [70], [71], and the dynamics of Jupiter's and Saturn's weather layers [72]. Surprisingly, these 2D systems retain intricate dynamics despite dimensionality reduction. Our findings here challenge the perception of reduced complexity in lower dimensions and suggests that energy transfer pathways can be more intricate than often assumed. For example, an external force may not always inject energy; it could pump energy into the cascade at an intermediate scale while simultaneously extracting energy at larger scales to balance the inverse cascade. A similar pattern was observed in oceanic surface flows under the Coriolis force [71]. Ultimately, the final status may depend on intricate small-scale interactions, as our results demonstrate. A systematical study of other two- and quasi-two-dimensional systems of turbulence will be reported in future.

3.5. Validity of direct numerical simulation (DNS)

Traditionally, DNS results of Navier-Stokes equations are widely regarded as 'reliable' benchmark solutions of turbulence, as long as grid spacing is fine enough (i.e. less than the minimum Kolmogorov scale) [4] and time-step is small enough, say, satisfying the Courant-Friedrichs-Lewy condition (Courant number < 1) [5]. Is this definitely true?
According to our CNS result, the maximum kinetic energy dissipation rate Dmax=5.095 appears at t=659.6 , corresponding to the minimum Kolmogorov scale [4], say,
ηminRe?3/4(Dmax)?1/4=0.042.
Note that Δx=Δy=2π/Nx0.0245 . Thus, the famous criterion on the grid spacing [4], i.e.
Δx=Δy<ηmin,
is satisfied for the CNS. It should be emphasized that, for the DNS result, the maximum kinetic energy dissipation rate Dmax is about 4.372, corresponding to the minimum Kolmogorov scale ηmin=0.043 , so that the grid spacing criterion (24) is also satisfied. Therefore, the spatial resolution adopted in this paper is fine enough for both of the CNS and DNS of the two-dimensional Kolmogorov turbulence under consideration.
Besides, it is well-known that a small enough time-step is needed for DNS of Navier-Stokes equations, say, the so-called Courant-Friedrichs-Lewy condition, i.e. Courant number < 1, must be satisfied [5]. For the two-dimensional Kolmogorov turbulence under consideration, the corresponding Courant number is 0.16 for CNS and 0.016 for DNS: both of them satisfy Courant-Friedrichs-Lewy condition! Therefore, the time-step used for both of DNS and CNS are small enough, too.
Unfortunately, even so, the DNS result has a huge deviation on both of small and large scales from the CNS benchmark solution not only in trajectory and spatial symmetry but also even in statistics, as reported above. Especially, as mentioned in Section 2.1, the true solution of the NS equations under consideration has the spatial symmetry (9), but unfortunately the DNS loses this kind of spatial symmetry when t>120 , while the CNS remains this kind of spatial symmetry in the whole interval of time t[0,1500] : this clearly indicates the large deviation of the DNS results (when t>120 ) from the true solution of the NS equations! This highly suggests that even the fine enough grid spacing satisfying the spacing criterion (24) with small enough time-step satisfying the Courant-Friedrichs-Lewy condition [5] could not guarantee the reliability and correction of DNS results: it is only a necessary condition but not a sufficient condition for the validity of the DNS! This also provides us a theoretical evidence that the hypothesis p+δ=p does not stand up for the two-dimensional Kolmogorov turbulence under consideration, where denotes an operator of statistics, p+δ is the result of the NS equations given by the DNS, p is the 'true' physical solution of the NS equations, respectively.
In addition, it is traditionally assumed that DNS results are correct and reliable if they are obtained using different time-step and/or grid spacing but match well each other in statistics. However, it is found that the statistic results given by the DNS using different time-steps match perfectly with each other, as shown in Fig. 18. But unfortunately, even such kind of perfect match still cannot guarantee the validity of the DNS results: they all quickly lose the spatial symmetry (9) and are quite different from those given by the CNS benchmark solution even in statistics! All of these challenge our traditional beliefs.
Fig.18 Comparisons of the probability density function (PDF) of (a) the kinetic energy E(x,y,t) and (b) the kinetic energy dissipation rate D(x,y,t) of the two-dimensional Kolmogorov turbulence in the case of nK=4 and Re=40 with the initial condition (7), given by the CNS (red line), the DNS (black line) with Δt=10?4 and the DNS' (blue symbol) with Δt=2×10?4 , respectively, where PDFs are integrated in x,y[0,2π] and t[200,1500] .

4. Concluding remarks and discussions

Traditionally, results given by the direct numerical simulation (DNS) of Navier-Stokes equations are widely assumed as 'reliable' benchmark solutions of turbulence, as long as grid spacing is fine enough (i.e. less than the minimum Kolmogorov scale) [4] and time-step is small enough, say, satisfying the Courant-Friedrichs-Lewy condition (Courant number < 1) [5]. Is this definitely true?
To answer this important fundamental question, a two-dimensional sustained turbulent Kolmogorov flow driven by an external body force is solved by means of the two numerical methods with detailed comparisons: one is the traditional 'direct numerical simulation' (DNS), the other is the 'clean numerical simulation' (CNS). The results given by the DNS are a kind of mixture of the 'false' numerical noise and the 'true' solution of the NS equations, which are mostly at the same order of magnitude, as shown in Fig. 1, Fig. 2 and Fig. 3. On the contrary, the numerical noise of the result given by the CNS is much smaller than the 'true' solution of the NS equations in a finite but long enough interval of time, i.e. t[0,1500] , so that the CNS result is very close to the 'true' solution and thus can be used as a benchmark solution of the NS equations for comparison. It is found that numerical noise as a kind of small-scale disturbances can lead to huge deviations in both of large and small scales on the two-dimensional Kolmogorov turbulent flow, not only qualitatively but also quantitatively (even in statistics), including the spatial symmetry of the flow field, the enstrophy spectrum as well as many statistics of physical quantities such as the velocity, vorticity, kinetic energy, kinetic energy dissipation rate and so on, as reported in Section 3.1. Besides, the field geometrical structures are also sensitive to the background numerical noise as a kind of artificial tiny disturbance, as illustrated in Section 3.2. In addition, it should be emphasized that the so-called intermittent stability of turbulence can be successfully found by the CNS, but is unfortunately lost by the DNS, as shown in Section 3.3. It is found in Section 3.4 that the averaged scale-to-scale energy flux given by the DNS has a large deviation from that given by the CNS benchmark solution. Furthermore, it is illustrated in Section 3.5 that a fine enough spatial grid with a small enough time-step alone cannot guarantee the correction/reliability and validity of the DNS: it is only a necessary condition but not sufficient. This also provides us a theoretical evidence that the hypothesis p+δ=p does not stand up for the two-dimensional Kolmogorov turbulent flow under consideration, where denotes an operator of statistics, p+δ is the result given by the DNS that might be “badly polluted”, p is the 'true' solution of the Navier-Stokes equations, respectively. So, unfortunately, DNS results of sustained turbulent flows, especially those having energy exchange with external systems, might have huge deviation from the 'true' solution of Navier-Stokes equations. Hopefully, CNS as a new tool to investigate turbulent flows more accurately than DNS could bring us some new discoveries. Of course, more investigations for various types of turbulent flows are needed in the future.
It should be emphasized that all results reported above is for the case of the Reynolds number Re=40 . Similarly, we investigated the two-dimensional Kolmogorov turbulent flow under the same initial condition (7) in ten different Reynolds numbers Re that are randomly chosen in the domain Re[35,50] , and obtained the similar results described above. Thus, all of our conclusions mentioned above have general meanings.
In summary, using a two-dimensional sustained turbulent Kolmogorov flow driven by an external body force as an example, we illustrated that a DNS result of the Navier-Stokes equations is not definitely correct/reliable even in statistics, even if the grid spacing and time-step are small enough! This finding might challenge some traditional beliefs and assumptions about the NS equations and turbulence, although they are worldwide accepted nowadays. Our finding highly suggests that one had better check very carefully DNS results of some complicated turbulent flows governed by the Navier-Stokes equations, especially those having energy exchange with external systems.
Certainly, a lots of investigations should be done in future. Below are some of our viewpoints for discussions.
First of all, it should be emphasized that the CNS result has a kind of spatial symmetry (9), which comes from the spatial symmetry of the initial condition (7) and the governing equation (1). As shown in Figs. 1 and 2, the flow field keeps such kind of spatial symmetry (9) for CNS in the whole interval of time t[0,1500] , but for the DNS only in a short interval of time from the beginning, i.e. t[0,120) , beyond which the random numerical noises of the DNS, which have no spatial symmetry at all, have been enlarged to the same order of magnitude as the “true” solution of the NS equations and thus finally destroy this kind of spatial symmetry (9) of the flow field. So, the loss of the spatial symmetry (9) of the DNS clearly indicates its large deviation from the “true” solution of the NS equations. In mathematics this provides us a rigorous evidence that small disturbances of the NS equations might be increased to a macroscopic level, say, at the same level of magnitude as the macroscopic physical variables under considerations. This confirms once again the previous conclusions reported by Lin, Wang & Liao [17] and Qin & Liao [12] for a Rayleigh-Bénard convection turbulent flow. However, this is contrary to the general belief that small disturbances of turbulent flow should be dissipated by the viscosity of fluid: it is based on this general belief that NS equations generally do not contain the unavoidable tiny noises such as thermal fluctuations, environmental disturbances, and so on. However, if NS equations are indeed an excellent model for turbulence, then our results given by the CNS reveals that its spatiotemporal trajectories are very sensitive to small disturbances, so that these small disturbances, no matter they are artificial or physical, must be considered, say, should not be neglected. Obviously, this leads to a logical paradox! Therefore, it should be wrong in physics for NS equations to neglect unavoidable small disturbances!
Note that nearly all DNS results of NS equations quickly become badly polluted. This is obvious from the well-known fact that, even for exactly the same NS equations with the same boundary/initial conditions and the same physical parameters, all DNS results quickly become quite different in spatiotemporal trajectories, although their animations of flow field look quite similar. Here, the obvious difference of spatiotemporal trajectories implies the 'badly' polluted, denoted by p+δ . However, if the statistics of the corresponding turbulent flows are stable in statistics, say, the hypothesis p+δ=p holds, then the badly polluted DNS results p+δ agree well with the true solution in statistics! This is the reason why many DNS results agree well with the experiments: the key point is the statistical stability, which is more important than the accuracy of numerical simulations of DNS.
Secondly, how can we know in general cases whether a DNS result of turbulent flow governed by the NS equations is correct/reliable in statistics or not? Although the two-dimensional sustained turbulent Kolmogorov flow considered in this paper is just a simple example, it highly suggests that the statistics stability, say, p+δ=p , should be very important for turbulent flows. This is exactly the reason why Liao [16] proposed the so-called 'modified fourth Clay millennium problem':
“The existence, smoothness and statistic stability of the Navier-Stokes equation: Can we prove the existence and smoothness of the solution of the Navier-Stokes equation with physically proper boundary and initial conditions, whose statistics are stable (or unstable), i.e., insensitive (or sensitive) to small disturbances? ”
Note that the statistic stability is emphasized here, since unavoidable physical disturbances (such as thermal fluctuations) always exist in practice. Statistic stability of turbulent flows is also the first problem of “Eight Key Open Questions in Ocean Engineering” [73] issued in 2023 by State Key Laboratory of Ocean Engineering, China. The key point is to find out the necessary and sufficient condition of such kind of statistic stability. From practical view-point of CFD, statistic stability of NS equations is more important than the existence and smoothness of its solution.
Thirdly, it should be emphasized that, from physical point of view, the NS equations (as a model of turbulence) neglect the influence of the unavoidable small physical disturbances such as thermal fluctuations and environmental noises. Note that physical thermal fluctuation is often larger than the artificial numerical noise of the DNS. Thus, if artificial numerical noises could be regarded as physical thermal fluctuations (see Eckmann and Ruelle [74]), our results imply that small but unavoidable physical disturbances sometimes might have huge influences on the long-term statistics of turbulent flow and thus should not be neglected in physics. Therefore, from physical point of view, it should be better to use the Landau-Lifshitz-Navier-Stokes (LLNS) equations [75], [76], which consider the influence of unavoidable thermal fluctuations, instead of the NS equations, to model turbulent flows. So, although the NS equations might be physically wrong since small unavoidable physical disturbances are not considered, artificial numerical noises in computational fluid dynamics (CFD) practically play a role like physical thermal fluctuations [74], which might “correct”, more or less, this kind of physical mistakes of the NS equations. Therefore, from physical viewpoint, artificial numerical noises might have positive meanings for modelling of turbulence, although they have no mathematical meanings at all. It should be emphasized that the NS equations are only a mathematical model of turbulent flows in physics: they are quite different things. Possibly, the unavoidable but nasty numerical noises might be very important in correctly modelling turbulence. It is a pity that the statistic stability and the influences of artificial numerical noises on modelling turbulence have been neglected.
Finally, let us point out the differences and relationships between the CNS and DNS. First of all, CNS contains an important concept, namely the “critical predictable time” Tc , within which numerical noise of numerical simulation (given by CNS) is much smaller than its true solution and thus is reliable (i.e. “clean”) and can be used as a benchmark solution. Secondly, the reliability of result in t[0,Tc] given by CNS with a numerical noise s1 can be verified by another result given by CNS with a smaller numerical noise s2<s1 . However, different from the CNS, DNS has not the concept “critical predictable time” Tc at all, since it is a traditional belief that small disturbances of NS equations should not enlarge to a macroscopic level. In other words, Tc of DNS is assumed to be infinite from the viewpoint of CNS. Unfortunately, our results reported in this paper indicate that this assumption is wrong! This is a fundamental difference between CNS and DNS. On the other side, the traditional DNS that uses data in the double precision and the 4th-order Runge-Kutta method in temporal dimension can be regarded as a special case of the CNS that adopts the multiple-precision and high-order Taylor expansion in temporal dimension. Note that the DNS result reported in this paper loses the spatial symmetry (9) and thus deviates largely from its true solution when t>120 , mainly because its “critical predictable time” Tc is only about 120, say, Tc120 , which is too short for calculating statistics. However, unlike DNS, the CNS can greatly enlarge the so-called “critical predictable time” Tc by reducing the numerical noises. For example, the “critical predictable time” Tc of the CNS result reported in this paper is greater than 1500, say, Tc>1500 , by means of the multiple precision with 200 significant digits for all physical/numerical variables and parameters and the 60th-order Taylor expansion with the time-step Δt=10?3 . In this way, one can gain a 'clean', reliable, benchmark solution of NS equations by means of the CNS. Shortly speaking, DNS can be regarded as a special case of the CNS: the CNS has much less numerical noises and thus much larger “critical predictable time” Tc than the DNS. In other words, CNS is more general than DNS. Indeed, for the very first time, CNS provides us a way to rigorously verify DNS results: any a CNS result with a “critical predictable time” Tc can be verified by another CNS result with even smaller numerical noises. Without doubt, CNS as a more general and more precise tool than DNS can provide us a totally new way to investigate turbulent flows.

Funding

This work is partly supported by National Natural Science Foundation of China (Grant No. 12302288; 12272230), Shanghai Pilot Program for Basic Research - Shanghai Jiao Tong University (No. 21TQ1400202), Postdoctoral Fellowship Program of CPSF under Grant Number GZC20231594, and State Key Laboratory of Ocean Engineering.
Data Availability Statement
The data that support the findings of this study are available on request from the corresponding author.

Author ORCID

Shijie QIN, https://orcid.org/0000-0002-0809-1766; Yongxiang HUANG, https://orcid.org/0000-0002-2011-9215; Shijun LIAO, https://orcid.org/0000-0002-2372-9502

Author contributions

Qin performed the analysis and wrote the draft. Yang performed the analysis. Mei and Wang analysed the results especially related to the field geometrical structures. Huang analysed the results, especially related to the scale-to-scale energy flux and the energy injection/extraction processes. Liao conceived/designed the analysis, and wrote the manuscript.

Declaration of Competing Interest

The authors report no conflict of interest. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. The CNS algorithm for two-dimensional Kolmogorov flow

In 2020, in order to overcome the huge computation requirements of the CNS algorithm in spectral space, an efficient CNS algorithm [17], [18], [19] was proposed in physical space for spatio-temporal chaos. Recently, Qin and Liao [12] have successfully applied the CNS algorithm in physical space to a two-dimensional turbulent Rayleigh-Bénard convection. Here, the basic idea of this kind of efficient CNS algorithm is briefly described by using the two-dimensional Kolmogorov flow as an example.
We rewrite the governing equation (1), which is based on Navier-Stokes equations, as following:
??t(ψxx+ψyy)=ψy(ψxxx+ψxyy)?ψx(ψxxy+ψyyy)+1Re(ψxxxx+2ψxxyy+ψyyyy)?nKcos(nKy),
where x and y as subscripts denote the spatial derivatives corresponding to horizontal and vertical directions, respectively.
The computational domain x,y[0,2π] is discretized by means of Nx×Ny equidistant, i.e.
xj=2πjNx,yk=2πkNy,
where j=0,1,2,...,Nx?1 and k=0,1,2,...,Ny?1 .
To reduce truncation errors in the dimension of time, the high-order Taylor expansions are adopted, i.e.
ψ(xj,yk,t+Δt)m=0Mψ[m](xj,yk,t)(Δt)m,
where Δt is the time-step, M is the order of Taylor expansion, with the definition
ψ[m](xj,yk,t)=1m!?mψ(xj,yk,t)?tm.
Here, the order M should be large enough so as to reduce the truncation errors to a required tiny level.
Differentiating (m?1) times both sides of (A.1) with respect to t and then dividing them by m! , we obtain the governing equation of ψ[m] :
ψxx[m](xj,yk,t)+ψyy[m](xj,yk,t)=1m{Fm+1Re[2ψxxyy[m?1](xj,yk,t)+ψxxxx[m?1](xj,yk,t)+ψyyyy[m?1](xj,yk,t)]+r=0m?1ψy[r](xj,yk,t)[ψxxx[m?1?r](xj,yk,t)+ψxyy[m?1?r](xj,yk,t)]?r=0m?1ψx[r](xj,yk,t)[ψxxy[m?1?r](xj,yk,t)+ψyyy[m?1?r](xj,yk,t)]},
where m1 and
Fm={?nKcos(nKyk),m=10,m>1.
Note that there exist some spatial partial derivatives denoted by subscripts in (A.5), such as ?s1+s2ψ[r]/(?xs1?ys2) with r,s1,s20 . In order to approximate these spatial partial derivative terms with high computational efficiency and precision from the known discrete variable ψ[r](xj,yk,t) , we adopt the spatial Fourier series
ψ[r](x,y,t)nx=??Nx/3??Nx/3?ny=??Ny/3??Ny/3?Ψ[r](nx,ny,t)exp(inxx)exp(inyy),
where nx , ny are integers, ?? stands for a floor function, i=?1 represents the imaginary unit, for dealiasing Ψ[r]=0 is imposed for wavenumbers outside the above domain , and
Ψ[r](nx,ny,t)=1NxNyj=0Nx?1k=0Ny?1ψ[r](xj,yk,t)exp(?inxxj)exp(?inyyk),
is determined by the known ψ[r](xj,yk,t) . Then, we can obtain the rather accurate approximations of the spatial partial derivative terms in (A.5), say,
?s1+s2ψ[r]?xs1?ys2(xj,yk,t)is1+s2nx=??Nx/3??Nx/3?ny=??Ny/3??Ny/3?(nx)s1(ny)s2Ψ[r](nx,ny,t)exp(inxxj)exp(inyyk).
Here, the fast Fourier transform (FFT) algorithm is used to increase computational efficiency. In this way, the spatial truncation error can be decreased to a required tiny level, as long as the mode numbers Nx and Ny are large enough.
Note that, if the order M of the Taylor expansion (A.3) is large enough, temporal truncation errors can be reduced to a required tiny level. Besides, if the spatial discretization is fine enough, say, the mode numbers Nx and Ny are large enough, spatial truncation errors in Fourier expression (A.7) and the corresponding spatial derivative terms in (A.5) can be decreased to a required tiny level. However, this is not enough, since there always exist some round-off errors that might be larger than temporal or spatial truncation errors. Therefore, all physical/numerical variables and parameters are expressed in multiple precision (MP) with a large enough number Ns of significant digits, so that the round-off errors can also be reduced to a required tiny level. In this way, the background numerical noise, i.e. truncation errors and round-off errors as a whole, can be reduced to a required tiny level. This is different from most of the other numerical algorithms, including the general DNS which has double precision at best. In addition, note that the results given by the CNS are useful only in a limited (but long enough) interval of time t[0,Tc] , in which the numerical noise can be negligible.

Appendix A. Supplementary materials

Supplementary Data S1. Supplementary Raw Research Data. This is open data under the CC BY license http://creativecommons.org/licenses/by/4.0/

[1]
U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov. Cambridge University Press (1995)

[2]
S. Peter, Explaining Chaos. Cambridge University Press (1998)

[3]
T. Bohr, M.H. Jensen, G. Paladin, A. Vulpiani, Dynamical Systems Approach to Turbulence. Cambridge University Press (2005)

[4]
S.B. Pope, Turbulent Flows. IOP Publishing (2001)

[5]
R. Courant, K. Friedrichs, H. Lewy. Math. Ann., 100 (1) (1928), pp. 32-74, DOI:

[6]
E.N. Lorenz, J. Atmos. Sci., 20 (2) (1963), pp. 130-141

[7]
G. Chen, X. Dong. From Chaos To Order. Co. World Scientific Pub. Co., Singapore (2014)

[8]
C.E. Leith, R.H. Kraichnan. J, Atmos. Sci., 29 (6) (1972), pp. 1041-1058

[9]
G. Boffetta, S. Musacchio. Phys. Fluids, 13 (4) (2001), pp. 1060-1062

[10]
G. Boffetta, S. Musacchio. Phys. Rev. Lett., 119 (5) (2017), p. 054102 DOI: 10.1103/PhysRevLett.119.054102

[11]
Y.C. Li, R.D.J.G. Ho, A. Berera, Z.C. Feng. J, Fluid Mech., 904 (2020), p. A27

[12]
S. Qin, S. Liao. J. Fluid Mech., 948 (2022), p. A7

[13]
N. Platt, L. Sirovich, N. Fitzmaurice. Phys. Fluids A, 3 (4) (1991), pp. 681-696

[14]
F. Feudel, N. Seehafer. Phys, Rev. E, 52 (4) (1995), p. 3506

[15]
S. Liao, Tellus Ser. A-Dyn. Meteorol. Oceanol., 61 (4) (2009), pp. 550-564 DOI: 10. 1111/j.1600-0870.2009.00402.x

[16]
S. Liao, Clean Numerical Simulation. Chapman and Hall/CRC, Boca Raton (2023)

[17]
Z. Lin, L. Wang, S. Liao. Sci. China-Phys. Mech. Astron., 60 (1) (2017), pp. 1-13

[18]
T. Hu, S. Liao. J. Comput. Phys., 418 (2020), p. 109629

[19]
S. Qin,S. Liao. Chaos Solitons Fractal., 136 (2020), p. 109790

[20]
J.C. Sprott, Elegant Chaos: Algebraically Simple Chaotic Flows. World Scientific, Singapore (2010)

[21]
W.-K. Lee, A.G.L. Borthwick, P.H. Taylor. J. Hydraul. Res.,52 (2) (2014), pp. 219-227 DOI: 10.1080/00221686.2013.855950

[22]
S. Gao, L. Tao, X. Tian,J. Yang. J. Fluid Mech., 851 (2018), pp. 687-714 DOI: 10.1017/jfm.2018.526

[23]
S. Liao, Chaos Solitons Fractal., 47 (2013), pp. 1-12

[24]
S. Liao, Commun. Nonlinear Sci. Numer. Simul., 19 (3) (2014), pp. 601-616

[25]
P. Oyanarte, Comput. Phys. Commun., 59 (2) (1990), pp. 345-358

[26]
S. Liao, P. Wang. Sci. China-Phys. Mech. Astron., 57 (2) (2014), pp. 330-335 DOI: 10.1007/s11433-013-5375-z

[27]
S. Liao, X. Li. Int. J, Bifurcat. Chaos, 25 (09) (2015), p. 1530023

[28]
X. Li, S. Liao. Sci. China-Phys. Mech. Astron., 60 (12) (2017), p. 129511

[29]
X. Li, Y. Jing, S. Liao. Publ. Astron, Soc. Jpn., 70 (4) (2018), p. 64

[30]
S. Liao, X. Li, Y. Yang. New Astronomy, 96 (2022), p. 101850

[31]
S. Liao, S. Qin. Adv. Appl, Math. Mech., 14 (4) (2022), pp. 799-815 DOI: 10. 4208/aamm.oa-2021-0364

[32]
S. Qin, S. Liao. Adv. Appl. Math. Mech., 15 (5) (2023), pp. 1191-1215

[33]
T. Xu, J. Li, Z. Li, S. Liao. Phys. Fluids, 33 (3) (2021), p. 037111

[34]
S. Qin, S. Liao. J. Fluid Mech., 960 (2023), p. A15

[35]
Y. Yang, S. Qin, S. Liao. Chaos Solitons Fractal., 167 (2023), p. 113037

[36]
B. Zhang,S. Liao. Physica D, 455 (2023), p. 133886

[37]
V.I. Arnold, L.D. Meshalkin, Usp. Mat. Nauk, 15 (247) (1960), pp. 20-24

[38]
L.D. Meshalkin, Y.G. Sinai. J, Appl. Math. Mech., 25 (6) (1961), pp. 1700-1705

[39]
D.N. Beaumont, J. Fluid Mech., 108 (1981), pp. 461-474

[40]
A. Thess, Phys. Fluids A, 4 (7) (1992), pp. 1385-1395

[41]
N.J. Balmforth, Y.-N. Young. J,Fluid Mech., 450 (2002), pp. 131-167 DOI: 10.1017/s0022111002006371

[42]
C. Marchioro, Commun. Math. Phys., 105 (1) (1986), pp. 99-106

[43]
J.S.A. Green. J. Fluid Mech., 62 (2) (1974), pp. 273-287

[44]
G.I. Sivashinsky, Physica D, 17 (2) (1985), pp. 243-255

[45]
G.J. Chandler, R.R. Kerswell. J,Fluid Mech., 722 (2013), pp. 554-595 DOI: 10.1017/jfm.2013.122

[46]
N.F. Bondarenko, M.Z. Gak, F.V. Dolzhanskii, Akad. Nauk SSSR, 15 (1979), pp. 1017-1026

[47]
A.M. Obukhov, Russian Math. Surveys, 38 (4) (1983), pp. 113-126

[48]
J. Sommeria, J. Fluid Mech., 170 (1986), pp. 139-168

[49]
J.M. Burgess, C. Bizon, W.D. McCormick, J.B. Swift, H.L. Swinney. Phys, Rev. E, 60 (1) (1999), p. 715

[50]
E.N. Lorenz, J. Atmos. Sci., 29 (2) (1972), pp. 258-265

[51]
E. Kazantsev, Nonlinear Process. Geophys., 5 (4) (1998), pp. 193-208 DOI: 10.5194/npg-5-193-1998

[52]
A.J. Manfroi, W.R. Young. J, Atmos. Sci., 56 (5) (1999), pp. 784-800

[53]
Y.-K. Tsang, W.R. Young, Phys. Fluids, 20 (8) (2008)

[54]
Y.-K. Tsang, W.R. Young. Phys. Rev. E, 79 (4) (2009), p. 045308 DOI: 10.1103/PhysRevE.79.045308

[55]
N.J. Balmforth, Y.-N. Young. J,Fluid Mech., 528 (2005), pp. 23-42 DOI: 10.1017/S002211200400271X

[56]
L. Wang,N. Peters. J. Fluid Mech., 554 (2006), pp. 457-475 DOI: 10.1017/S0022112006009128

[57]
L. Wang,N. Peters. J. Fluid Mech., 608 (2008), p. 113138

[58]
L. Wang, N. Peters. Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci., 371 (1982) (2013), p. 20120169 DOI: 10.1098/rsta.2012.0169

[59]
D. Denker, A. Attili, J. Boschung, F. Hennig, M. Gauding, M. Bode, H. Pitsch. J. Fluid Mech., 905 (2020), p. A4

[60]
F.D. Costa, O.C. Acevedo, J.C.M. Mombach, G.A. Degrazia. J, Atmos. Sci., 68 (8) (2011), pp. 1714-1729

[61]
G. Gogia, J.C. Burton. Phys. Rev, Lett., 119 (17) (2017), p. 178004 DOI: 10.1103/PhysRevLett.119.178004

[62]
G. Gogia, W. Yu, J.C. Burton. Phys, Rev. Res., 2 (2) (2020), p. 023250

[63]
A. Leonard, Adv. Geophys., 18 (1975), pp. 237-248

[64]
S. Chen, R.E. Ecke, G.L. Eyink, X. Wang, Z. Xiao. Phys. Rev. Lett., 91 (21) (2003), p. 214501

[65]
G. Boffetta, R.E. Ecke. Annu. Rev, Fluid Mech., 44 (2012), pp. 427-451, DOI:

[66]
Q. Zhou, Y. Huang, Z. Lu, Y. Liu,R. Ni. J. Fluid Mech., 786 (2016), pp. 294-308 DOI: 10.1017/jfm.2015.673

[67]
R.H. Kraichnan, Phys. Fluids, 10 (1967), pp. 1417-1423, DOI:

[68]
L. Wang, Y. Huang. Phys, Rev. E, 95 (2017), p. 052215, DOI:

[69]
G.K. Vallis, Atmospheric and Oceanic Fluid Dynamics. Cambridge University Press (2017),

[70]
S. Lovejoy. Weather, Macroweather, and the Climate:Our Random Yet Predictable Atmosphere. Oxford University Press (2019),

[71]
D. Zhang, J.J. Zhang, Y. Gao, Y. Peng, J.Y. Hu, F.G. Schmitt, Y.X. Huang, Front. Marine Sci. (2024) ( accepted)

[72]
P.L. Read, Annu. Rev. Fluid Mech., 56 (2024), pp. 271-293 DOI: 10.1146/annurev-fluid-121021-040058

[73]
S.K.L. of Ocean Engineering, J. Ocean Eng. Sci., 8 (6) (2023), pp. iv-vi, DOI:

[74]
J.-P. Eckmann, D. Ruelle. Rev. Mod, Phys., 57 (1985), p. 617

[75]
L.D. Landau, E.M. Lifshitz, Course of Theoretical Physics: Fluid Mechanics ( Vol. 6). Addision-Wesley, Reading (1959)

[76]
D. Bandak, N. Goldenfeld. Phys, Rev. E, 105 (2022), p. 065113

Outlines

/