Energy conservation uncertainly due to the use of an implicit time integration method clarified by a multiplexed TGV inviscid flow

This study attempts to elucidate the conservation errors induced by the implicit time integration method within a multiscale flow. Using the multiplexed Taylor analytical solution as the initial field, three distinct multiscale turbulence fields were formulated, each characterised by different wavenumber attributes. A two-dimensional inviscid flow field, replicated in a computational domain with a side of 2π and subject to periodic boundary conditions, was used to illuminate kinetic energy conservation errors. The implicit time integration method - specifically the second-order accurate Crank-Nicolson method - was compared with explicit methods, such as the second-order accurate Adams-Bashforth and third-order accurate Runge-Kutta methods. The study also examines the differences between the analytical Taylor solution and a random flow field, particularly in their satisfaction of the Navier-Stokes equations and wavenumber composition, while highlighting the need for preliminary analysis for random flow fields prior to validation calculations.


Introduction
Understanding fluid phenomena is a widespread need in engineering.Experimental measurements in wind tunnels and numerical calculations are widely used for this purpose.In wind tunnel measurements, multi-scale turbulence is analysed using fractal grids in addition to square grids.In addition, studies to clarify turbulence phenomena using numerical analysis are also widely carried out.Numerical analysis is widely used because of its advantages, such as the ability to measure physical quantities that are difficult to measure experimentally.However, numerical analysis requires knowledge of the accuracy of the simulations because of the analytical errors involved.In order to achieve highly accurate analyses, conservation laws should be established with a high degree of accuracy in numerical analysis [1].For incompressible fluids, in addition to the mass and momentum conservation laws as governing equations, the kinetic energy conservation laws derived from these equations should also be satisfied.Generic fluid analysis solvers (e.g., [3]) are widely used in engineering applications because of their ability to analyse complex geometries, however these conservation laws may not be established with sufficient accuracy.Therefore, additional knowledge is required about the accuracy of the analysis when errors occur in the conservation laws.
This study focuses on the energy conservation error [2].This error can occur even when the continuity equation and the Navier-Stokes equations are solved discretely.Schemes have been developed to explicitly remove this error and allow the flow field to evolve.Initial value problems in inviscid flows are used as a way to validate this analytical approach.The time evolution of the initial flow field under inviscid conditions is studied.Note that under inviscid conditions the turbulent kinetic energy is not dissipated at all and should be constant in time.The initial flow field used is a Taylor analysis solution [4], based on the Taylor-Green vortex [5], or a random field [2].Here the random field does not satisfy the Navier-Stokes equations, but the continuity equation does.Using this flow field, the conservation properties of the analysis code based on the finite volume method have been investigated.
This study focuses on the difference in the properties of these two initial fields.There are differences between the flow field from the analytical Taylor solution and the random flow field.First, the Taylor analytical solution satisfies the Navier-Stokes equations, whereas the random flow field does not.For this reason, a random flow field is expected to require an analysis to develop the flow field prior to the validation calculations.This random flow field consists of waves of different frequency periods.In particular, the spectrum of the random flow field is constant with respect to the wavenumber.On the other hand, the analytical Taylor solution consists of a wave with a single wavenumber.The actual turbulent flow field consists of waves with wavenumbers that satisfy the Navier-Stokes equation and span a wide bandwidth.Therefore, these two initial fields differ in some respects from the properties of the actual turbulence field.
This study addresses this point by multiplexing the analytical Taylor solution to include waves with wavenumbers over a wide bandwidth in the initial flow field.This flow field can satisfy the Navier-Stokes equations.The aim of this study is to use this flow field to investigate energy conservation errors.Specifically, energy conservation errors occur when the implicit integration method is used in time, as opposed to the explicit integration method.

Numerical Methods
The present analysis is carried out for a two-dimensional incompressible flow field.The computational domain is a square with one side of 2π [2,4].The boundary condition for velocity and pressure is set to periodic.In the present analysis, the flow field is divided into 32 sections in the xand y-directions respectively.The governing equations are the 2D incompressible flow continuity equation and the Navier-Stokes equations.In the Navier-Stokes equations, the viscous terms are removed and an inviscid field is formed.This method does not involve numerical viscosity.The differential method for the convection terms is of the skew-symmetric type [2], a method that naturally conserves kinetic energy.The velocity-pressure coupling method uses the fractional step method.This flow field is a multi-scale turbulent flow field containing multiple wavenumbers and satisfying the equations.
A multiplexed Taylor-Green vortex is used for the initial field.In previous studies [6], multiscale turbulence grids with multiplexing of a single wavenumber are used to generate decaying turbulence in a wind tunnel experiment.In the present study, a multiscale initial field is similarly generated by multiplexing a single-scale analytical solution.The analytical solution used is Taylor's analytical Here, x, y, t are nondimensional coordinates and time, respectively, and Re in the equation is the Reynolds number.By setting these values, multiscale turbulence with three different wavenumber characteristics is generated.The multiplexed Taylor analytical solution used in this study is also inspired by the previous studies of wind tunnel experiments.The fractal turbulence grids used to generate decaying turbulence are obtained by multiplexing a square turbulence grid [6].The basic elements that give rise to fractal turbulence grids are not only square turbulence gratings, but also square-shaped gratings.Based on previous experiments on such multiplexing, this study generates an initial field close to the actual turbulence field by multiplexing the analytical Taylor solution.
The explicit method is a method in which the calculation of the next time step is performed from the current time when the calculation has been completed.The explicit method has the advantage of small computational capacity due to the simplicity of each individual calculation, but has the  disadvantage that the time step should be set small for a stable analysis, thus increasing the computational time.The implicit method uses the current value and the value of a next time step when calculating a next time step from the current time.The implicit method does not require as small a time step width as the explicit method and is therefore considered to have a shorter overall computation time.
In the present analysis five time steps of t = 0.03, 0.01, 0.001, 0.003 and 0.0003 are used with a time evolution up to t = 30.The magnitude of the initial field is normalised so that the spatial mean of the kinetic energy is unity.The results are compared with the Adams-Bashforth method, which is an explicit time integration method, and their values are verified using the Runge-Kutta method, which is an explicit time integration method with third-order accuracy.

Results and Discussion
In this analysis, the visualisation results of the flow field at initial time are shown first.Figure 1(a) shows the distribution of the kinetic energy at the initial time, including all wave numbers.The kinetic energy is distributed in an X-shape, with the distribution centred on the central axis.The all wavenumbers analysis includes waves with high wavelengths, so parts of the flow field show abrupt changes in velocity.Figure 1(b) shows the distribution of kinetic energy at the initial time when the wavenumber is halved and Figure 1(c) when the wavenumber characteristic is reduced to 1/4.As in Fig. 1(a), the kinetic energy is distributed in an x-shape with the distribution centred on the central axis.For the full wavenumber case, the distribution of the high kinetic energy identified is a narrow x-shape.For the half wavenumber case, the high kinetic energy distribution is wider than in Figure 1(a).Similarly, the width of the distribution is wider when the wavenumber is a quarter.The narrower distribution in the high kinetic energy region indicates that the velocity difference in the flow field occurs rapidly.These results indicate that the velocity difference in the flow field is more abrupt at full wavenumber and that decreasing the wavenumber results in a slower velocity difference in the flow field.The rapid velocity difference also indicates the presence of high wavenumber waves in terms of wavenumbers.Therefore,   although high wavenumber waves are present at all wavenumbers, the low wavenumbers indicate that high wavenumbers are decreasing and low wavenumbers are increasing.
To verify the accuracy of the calculations at each wavenumber, time series of the values of the continuity equations at each wavenumber are shown.Figure 2 shows the time series of the value of the continuity equation for each time integration method for each wavenumber condition.Figures 2(a), 2(b) and 2(c) show the results for the full, half and quarter wavenumber conditions respectively.For Figure 2(a), the values of the continuity equation are of the order of the -11th power of 10 for all time integration methods.This value is considered sufficiently small to indicate that the calculations are satisfactory for all wavenumber conditions.The value of the continuity equation is always small from the initial time, indicating that the calculation converges from the initial time.This can be attributed to the fact that the initial field satisfies the equations.Similarly, in the case of Figures 2(b) and 2(c), the values of the continuity equation are of the order of the -11th power of 10 for all time integration methods, indicating that the calculations have been carried out accurately.These results verify the validity of the calculations for all time integration methods for the three wavenumber conditions.
Figure 3(a) shows the kinetic energy distribution for the Adams-Bashforth method and Figure 3(c) for the Runge-Kutta method.Figure 3(b) shows the distribution of kinetic energy for the Crank-Nicolson method, an implicit time integration method.Figures 3(a) and 3(c) show a rotation of the flow field.On the other hand, the kinetic energy distribution obtained using the implicit time integration method, the second-order Crank-Nicholson method, is significantly different from these results.In Figure 3(b), the distribution of most of the energy near the central region as seen in Figures 3(a) and 3(c) is absent.
The time evolution of the kinetic energy for each of the time integration methods is then shown in Figure 4.As the flow field in the present analysis is inviscid, there is no viscous dissipation and the kinetic energy analytically retains its initial value.The initial value of the kinetic energy is normalised to unity.Figure 4 shows the results of the time evolution of the kinetic energy for each wavenumber condition.These results establish that the implicit time integration method reduces the kinetic energy conservation property and that this tendency is not affected by the wavenumber property of the initial field.
Figure 5 shows the dependence of the slope of the magnitude of the error in the analytical solution of the kinetic energy at t = 0 on the time step width.The slope of the error at t = 0 is obtained from a polynomial approximation of the time series with a fourth order expression.The slope of the fourth order approximation of the time series at t = 0 is the coefficient of the linear term.The slope at t = 0 is calculated by a fourth order approximation of the time series under each condition.The figures show the results for each wavenumber condition.From Figure 5(a) the results obtained using the second-order accurate Adams-Bashforth method are proportional to the square of the time step width.The results obtained using the third-order Runge-Kutta method are also proportional to the cube of the time interval.These results confirm that the order of accuracy of the respective time integration methods coincides with the order of the time steps, and this trend is also maintained in Figures 5(b) and 5(c).Then the results obtained using the Crank-Nicolson method are examined.From Figure 5(a) it can be seen that the results obtained using the second-order accuracy Crank-Nicolson method are of the order of the time increment range.

Conclusion
To form the multiscale flow field, Taylor's analytical solution was multiplexed and set as the initial field.Three different multiscale turbulence fields were formed by setting three different wavenumber characteristics.An inviscid flow field was used to clarify the kinetic energy conservation error.The flow field was two-dimensional, and the flow field was reproduced in a computational domain with one side of 2π.The periodic boundary condition is used for each side.The governing equations are discretised by the central difference method so that the numerical viscosity is not included and the kinetic energy is analytically conserved.
In this study, the differences in the initial field due to different wavenumber characteristics were first visually confirmed.The flow field containing higher wavelengths was visually confirmed to have more rapid velocity differences in the flow field, visually confirming the differences in the flow field due to wavenumber characteristics.The identified multi-scale flow field was time evolved using three different time integration methods.The accuracy of the calculations was verified using the continuity equation.For all flow fields and time integration methods, the values of the continuity equation are small and the calculations are found to be sufficiently accurate.

Figure 1 .
Figure 1.Dependence of the multiplexed initial flow field on the computational domain and spatial resolution.

Figure 2 .
Figure 2. Temporal development of the continuity-equation error.

Figure 3 .
Figure 3.Time evolution of the flow fields from the initial field presented.

Figure 4 .
Figure 4. Developing deviation of the kinetic energy from the initial value.

Figure 5 .
Figure 5.The deviation of the initial derivative of the kinetic energy.