Study of the Three Dimensional Eular Equations

One of the most difficult problems in computational fluid dynamics that has been addressed is whether there is evidence for or against the existence of a singularity of the three-dimensional incompressible Euler equations. Since two publications in 1987 in Physical Review Letters (Ash- urst and Meiron; Pumir and Kerr), the favored approach to this problem has been direct, Eulerian simulation of an initial value problem with two prescribed vortices. One of the major points of these letters was that while Lagrangian vortex methods, are adequate for representing the initial phase of development due to vortex core deformation, they are inadequate for representing the late phases when the question of singularity arises. Since then, there has developed a growing body of evidence from a variety of calculations that suggest there is no singularity and that at late times the growth in the peak vorticity is, at most, exponential. The most detailed investigation is a series of inviscid simulations of antiparallel vortices by Pumir and Siggia. Inviscid calculations with random initial conditions in a periodic box with 256 resolution could also be interpreted as supporting exponential growth. In contrast to these inviscid calculations, there have been a few viscous calculations that have suggested a singularity." The major difficulty with using a series of viscous calculations to infer trends toward the limit of zero viscosity is that the range of Reynolds numbers that can be used is too small. In addition, for at least one similar problem, whether there is a singularity of the inviscid, two- dimensional, incompressible, magnetohydrodynamic equations, it has been demonstrated that the trend toward a singularity indicated by viscous calculations could not be reproduced by inviscid calculations. Therefore, due to the current weight of evidence against a singularity, any claim for a singularity must match the computational power of the calculations finding exponential growth, be able to reproduce the exponential behavior for similar initial condi- tions, and from this suggest a reason for at least two types of behavior, exponential and singular. In addition, the growth of various quantities such as the peak vorticity, that is, the norm of the vorticity, and peak strains should be consistent with analytic constraints for the Euler equations Evidence for a singularity in an inviscid calculation that meets these requirements will be presented in this paper. The result is consistent with the conclusion of Kerr and Hussain that with their initial conditions, viscous calculations gave signs of singularity behavior. In order to investigate the question of a singularity of Euler, it is necessary to have a well-resolved numerical method and initial conditions that are predisposed toward possible singular behavior The first two sections of this paper contain a brief discussion of the numerical method— details appear in the Appendix—and a complete discussion of the initial conditions, including subtle, but significant, differences with earlier work. Following this introduction, the analysis of the calculation is broken into four major parts. First, symmetry plane contour plots and three- dimensional isosurface plots of the vorticity will be used to show the basic development. Second singular development will be demonstrated by plotting in several ways peak vorticity, peak axial strain, enstrophy, and enstrophy production as functions of time. These plots will show that the trends are consistent with the analytic limits12-14 and that this behavior is maintained far beyond the time of similar behavior in Pumir and Siggia. That is, singular behavior is maintained after the time when the calculations of Pumir and Siggia shift to exponential, nonsingular growth. Included is an analysis of the values and positions of the important components of the strains that demonstrates self-similar behavior consistent with singular growth. This analysis is designed to meet the objection that, even for this initial condition, the flow might eventually shift to expo- nential growth and filtering the initial condition just provides an ingenious way to delay the inevitable behavioral shift. In the next section there is a detailed structural analysis via contour plots and three-dimensional surface plots to show how the curvature in three-dimensions could be sustaining the singular behavior via Biot-Savart interactions. In the final analysis section, there is additional discussion of resolution limits using curves of peak vorticity and axial strain at different resolutions and spectra. Dynamical information contained in spectra is also discussed.

To provide an overview of the structural development Figs. 3(a)-3(h) show three-dimensional isosurface plots of the vorticity squared. The x-z plane in front of each isosurface plot is through the symmetry plane. None of the isosurface figures shows the entire domain. At early times, Figs. 3(a)-3(c), t=0-9, a larger subdomain containing the entire region of the initial perturbation is shown. At late times, Figs. 3(d)-3(h), only the region of maximum vorticity is shown. Figures 3(c) and 3(d) show the vorticity at the same time t=9, in the different domains, to allow a comparison between the two volumes used. Off the symmetry plane, Fig. 3(a) shows the appearance of more anomalous regions of vorticity that result from the initial filtering. Based on this figure, the maximum vorticity of the most obvious leg is about 0.2. At later times contour plots show that the positions of these anomalous regions do not change with respect to the major vortex and have about the same values, showing that once the peak vorticity has grown sufficiently these relics of the initial conditions do not affect the development in the vicinity of the peak vorticity. This is supported by Fourier calcu- lations with less filtering in the initial conditions that do not have these anomalous regions, but still show singular development. To illustrate the development in the symmetry plane more clearly, Figs. 4(a)-4(c) show a series of late times with the axes stretched in a manner similar to Pumir and Siggia.5 The head and tail in the symmetry plane that come from the dynamics, not from the initial conditions, develop from a tilted, flattened region into a dimple in Fig. 4(a) at t—12.0. At later times this dimple becomes a corner separating the head and tail. The definition of the "head" will be the region extending above the vorticity peak just behind the leading edge of the vortex. The "tail" will be the vortex sheet extending behind the peak vorticity. By stretching the coordinates to match the degree of collapse in the different directions at a given time, then following the structure in time, it should be possible to determine whether one direction becomes progressively more flattened or stretched as the flow evolves and whether the rates of collapse are the same in different directions. For singular behavior, one expects that the rates of collapse in all three directions will be the same. By stretching z it can be seen that both the head and tail are vortex sheets and remain so, where in earlier work only the tail seemed sheet like. As the flow evolves the region of peak vorticity concentrates into the region where the vortex sheets of the head and tail meet at a sharp angle. In Pumir and Siggia, these features are not observed. Instead, the head of their structures flattens into the vortex sheet of the tail and the highest contours are not as concentrated. In viscous calculations, the peak vorticity was located in the head. Figures 4(a)-4(c) indicate equal rates of collapse in x and z up to t= 17, which shows that as this flow evolves the only changes are that the structures become better defined and more localized. Contour plots beyond t= 17 show a continuation of these trends in the immediate vicinity of the peak vorticity, but there is too much numerical noise outside that region to justify inclusion in this paper. Three-dimensional isosurface plots in Figs. 3(e)-3(h) also show the flattening of the initial vortex tubes with the appearance of a head and tail by t= 15. These figures show tilting in y toward the dividing plane, with the tilted vortex coming in at the top of these figures from where it connects to the main horizontal vortex. Some focusing toward the tip of the vorticity in the symmetry plane, indicated by arrows in Figs. 3(g) and 3(h), occurs with the head seeming to increase in radius as it moves away from the symmetry plane.

The structure of the head in the symmetry plane cannot be identified at the final two times in Figs 3(g) and 3(h) due to the resolution of the graphics. Filaments coming off the back of the tail should have no influence for the development around the peak vorticity.

Evidence has been presented that a finite-time singularity of the incompressible Euler equations develops from a smooth initial condition in a finite geometry. The strongest evidence for a singularity comes from the time dependence of the peak vorticity and the vorticity production between t—15 and t= 17. Extensive resolution comparisons indicate that these quantities, along with the peak axial strain eyy.p and the enstrophy are well converged and consistent with singular behavior up to the final time calculated, even though at a fairly early stage, by t= 12.5 spectra begin to be under-resolved. To complement this primary evidence, indications of self- similar behavior are presented. Equal rates of collapse in all three directions are indicated by following ratios and locations of maxima and minima of components of the strain and by comparing structures in two- and three-dimensional plots at different times. Singularities in axisymmetric flows, which have recently been claimed, are not related to the singular behavior reported here, in part because those flows cannot have collapse in all directions as a property. It is this self-similar behavior with strong curvature in all directions that allows the strain to increase along with the vorticity through Biot-Savart interactions, and prevents the strain saturation that has led to exponential growth in other investigations of Euler singularities.5 Without showing this self-similar behavior the possibility that the structures could evolve toward a nonsingular form would exist. Resolution comparisons show that the basic self-similar structure is not affected by the possible resolution limitations indicated by spectra in Figs. 18(a)-18(c). Consistent with algebraic, power-law, singular behavior for the peak vorticity, reasons why the

axial strain at the peak vorticity will scale as l/(tc—t) are discussed, but the theoretical prediction12

for the peak vorticity is less constrained, where Reasons are given why one might expect thatand plots of the inverses of emp, and are used to demonstrate this A singular time tc is predicted by fitting the peak vorticity , peak axial strain eyyp, and enstrophy production to the form 1/[tc—t]. The evidence suggests that while the singularity of the peak vorticity occurs at a point the singularity is not on a set of measure zero because the enstrophy also blows up. But the enstrophy blowup is logarithmic, so this is not the singularity predicted by spectral closures, which predict that the enstrophy scales as 1/t2. For the enstrophy to blow up, the kinetic energy spectrum must be no steeper than k -3 at the singular time. From the time dependence of the slope of energy spectra and the observation of logarithmic enstrophy blowup, it is suggested that the spectrum at the singular time might be identically k-3. The time dependence of the peak axial strain eyy,p shifts to the correct asymptotic behavior in Fig 6(b) only after Z=17. This shift does not appear to have any effect upon the overall development, but the possibility exists that,

FIG. 18. Spectra or Chebyshev distributions in (a) x, (b) y and (c) z, and (d) time dependence of slopesTimes are 2 (0),4 9 ( + ), 12.5 (X), and 18

when starting from more general initial conditions, the flow could follow several stronger false asymptotes before developing the correct singular behavior. Indications of nonsingular behavior in development from random initial conditions and Taylor-Green flow might be because those flows are following incorrect initial asymptotes. At least in Taylor-Green, after a period of skewness saturation, there is an indication that when the skewness begins to rise again that a new type of behavior occurs that could be singular. One way to investigate whether there are singularities in these flows would be to compare with viscous calculations of vortex reconnection. In Kerr and Hussain, it was concluded from a series of viscous calculations using an initial conditions similar to the one used here that their solutions would show a singularity in a finite time in the inviscid limit. Reasons have been given why a series of viscous calculations cannot provide a reliable prediction of singular behavior, but with the completion of this calculation, inviscid support for their conclusion has been provided. Their calculation therefore provides us with a picture of how viscosity would regularize the singular behavior found here and what the development after the time of an inviscid singularity may be. These effects should be looked for in other viscous calculations showing singular behavior, and work is in progress showing the relationship of the singular behavior found here to events at late times in calculations that were thought not to show singular behavior. How does this simulation compare with what is known about homogeneous, isotropic numerical turbulence in a periodic box, the turbulent flow we know the most about and a paradigm for many turbulent theories? Without a priori knowledge of decaying and forced simulations of isotropic homogeneous turbulence, the simulations reported here would predict that when starting from random initial conditions a singularity will develop, although it might not dominate the initial behavior. Furthermore, if it were assumed there is a finite probability of these events occurring homogeneously in space and time in the inviscid limit, then time averages of properties such as enstrophy production and fourth moment of the vorticity in recon- nection calculations should produce many of the statistics of isotropic calculations. Based upon the time dependence of the enstrophy production and fourth moment of vorticity as the singular time is approached, one would then expect Reynolds number dependence of both the velocity derivative skewness and flatness, with the skewness dependence being weaker. This is what is reported experimentally (see the references in Van Atta and Antonia). In addition, because the components of the strain have the same time dependence as the vorticity as the singular time is approached, one would expect that all flatnesses of fourth-order combinations of velocity derivatives would have the same Reynolds number dependence. The magnitudes of these derivative flatnesses could be very different, just as the magnitude of the peak vorticity is much larger than the magnitudes of the strain components, but they should scale with Reynolds numbers in the same way. Few of these properties are found in simulations of isotropic, homogeneous turbulence. There is no evidence from the inviscid calculations6 for a singularity, and this seems associated with vortex flattening, saturation of the strain, and saturation of the velocity derivative skewness near a value of 0.5. This is consistent with a long history of decaying and forced simulations of isotropic, homogeneous turbulence, now reaching on a 512 mesh, where the steady-state value for the skewness is about 0.5 and independent of Reynolds number. A prediction of the simulation of forced turbulence is that different fourth- order derivative flatnesses should have different Reynolds number dependencies. One's first prejudice would be that these differences are because the Reynolds number of the isotropic simulations is not large enough. But given the enormous consistency of calculations of homogeneous, isotropic turbulence over the years, it seems unlikely that low Reynolds number is the only reason for these differences. It seems that there is some important physics that we do not yet understand in going from singularities and vortex re- connection to fully developed turbulence. The evidence would suggest that in homogeneous, isotropic numerical turbulence in a box the probability of nearly singular events occurring is zero. In addition to the singular solution, another initial condition with nonsingular, exponential growth has been discussed. Numerous other calculations have also observed exponential growth. Which of the two classes of solutions is more relevant to understanding turbulence? It is possible that both are important. The mechanism that suppresses the singularity for the unfiltered initial condi- tion might have a role in suppressing the singularity in homogeneous turbulence. The dynamics of the unfiltered initial conditions might also be relevant to viscous boundary layers, where a region of negative vorticity between the primary vortex and the wall, similar to those appearing from the unfiltered initial conditions, does form and move in front of the primary vortex, leading to what is known as the bursting phenomena. The resolution tests in this study generally support the numerical accuracy of these calculations but some areas of concern have been noted, in particular, power-law spectra reaching the highest wave numbers at an early stage. It would be desirable to reproduce and possibly extend the singular behavior with a calculation that uses an entirely different numerical method, for example, a nested mesh finite-difference code.5 Hopefully, vorticity would develop the same asymptotic behavior and the signs of self-similar growth could be improved. The major difficulty with a finite-difference code could be finding an initialization procedure that included filtering in physical space. The success of wavelet filtering of initial conditions in two dimensions suggests that as a possible procedure. If the initialization problem can be overcome, the initial objective would be over-resolution during the phase equivalent to the early stages of this calculation. Once this has been accomplished, it could then be discovered whether the singular behavior at late times found here can be reproduced.

1. W. Ashurst and D. Meiron, "Numerical study of vortex reconnection," Phys. Rev. Lett. 58, 1632 (1987). 2. A. Pumir and R. M. Kerr, "Numerical simulation of interacting vortex tubes," Phys. Rev. Lett. 58, 1636

(1987).

3. A. Chorin, "The evolution of a turbulent vortex," Commun. Math. Phys. 83, 517 (1982). 4. E. D. Siggia, "Collapse and amplification of a vortex filament," Phys. Fluids 28, 794 (1985). 5. A. Pumir and E. D. Siggia, "Collapsing solutions to the 3-D Euler equations," Phys. Fluids A 2, 220

(1990).

6. J. R. Herring and R. M. Kerr, "Development of enstrophy and spectra in numerical turbulence," submitted to Phys. Fluids A (1992). 7. R. M. Kerr and F. Hussain, "Simulation of vortex reconnection," Phys- ica D 37, 474 (1989). 8. O. N. Boratav, R. B. Pelz, and N. J. Zabusky, "Reconnection in orthogonally interacting vortex tubes Direct numerical simulations and quantification in orthogonally interacting vortices," Phys. Fluids A 4, 581

(1992).

9. D. L. Marcus and J. B. Bell, "Vorticity intensification and transition to turbulence in the three-dimensiona Euler equations," Commun. Math. Phys. 147, 371 (1992). 10. S. A. Orszag and C. M. Tang, "Small-scale structure of MHD turbulence," J. Fluid Mech. 90, 129 (1979). 11. U. Frisch, A. Pouquet, P. L. Sulem, and M. Meneguzzi, "The dynamics of two-dimensional ideal MHD," J M6c. Theor. Appl. D 2, 191 (1983). 12. J. T. Beale, T. Kato, and A. Majda, "Remarks on the breakdown of smooth solutions of the 3-D Euler equations," Commun. Math. Phys. 94, 61 (1984). 13. A. J. Majda, "Vorticity, turbulence and acoustics in fluid flow," SIAM Rev. 33, 349 (1991). 14. G. Ponce, "Remark on a paper by J. T. Beale, T. Kato and A. Majda," Commun. Math. Phys. 98, 349

(1985).

15. R. D. Moser, P. Moin, and A. Leonard, "A spectral numerical method for the Navier-Stokes equations with application to Taylor-Couette flow," J. Comput. Phys. 52, 524 (1983). 16. R. M. Kerr, "Evidence for a singularity of the three-dimensional incompressible Euler equations," in Topological Aspects of the Dynamics of Fluids and Plasmas, Proceedings of the NATO-ARW Workshop at the Institute for Theoretical Physics, University of California at Santa Barbara, edited by H. K. Moffatt, G. M Zaslavsky, M. Tabor, and P. Comte (Kluwer-Academic, Dordrecht, The Netherlands, 1991), p. 309. 17. M. V. Melander (private communication). 18. M. V. Melander and F. Hussain, "Cross-linking of two antiparallel vortex tubes," Phys. Fluids A 1, 633

(1989).

19. A. Leonard, "Vortex methods for flow simulation," J. Comput. Phys. 37, 289 (1980). 20. M. E. Brachet, D. I. Meiron, S. A. Orszag, B. G. Nickel, R. H. Morf, and U. Frisch, "Small-scale structure of the Taylor-Green vortex," J. Fluid Mech. 130, 411 (1983). 21. M. E. Brachet, "Direct simulation of three-dimensional turbulence in the Taylor-Green vortex," Fluid Dyn Res. 8, 1 (1991). 22. A. Pumir and E. D. Siggia, "Vortex dynamics and the existence of solutions of the Navier-Stokes equations," Phys. Fluids 30, 1606 (1987). 23. J. A. Domaradzki and R. S. Rogallo, "Local energy transfer and nonlocal interactions in homogeneous isotropic turbulence," Phys. Fluids A 2, 413 (1990). 24. R. M. Kerr, "Velocity, scalar and transfer spectra in numerical turbulence," J. Fluid Mech. 211, 309 (1991). 25. R. Grauer and T. Sideris, "Numerical computation of three dimensional incompressible ideal fluids with swirl," Phys. Rev. Lett. 67, 3511 (1991). 26. A. Pumir and E. D. Siggia, "Finite time singularities in the axisymmet- ric three dimensional Euler equations," Phys. Rev. Lett. 68, 1511 (1992). 27. A. Pumir and E. D. Siggia, "Development of singular solutions to the axisymmetric Euler equations," Phys Fluids A 4, 1472 (1992). 28. M. Lesieur, Turbulence in Fluids (Martinus Nijhoff, Dordrecht, The Netherlands, 1987). 29. C. W. Van Atta and R. A. Antonia, "Reynolds number dependence of skewness and flatness factors of turbulent velocity derivatives," Phys. Fluids 23, 252 (1980). 30. S. Chen, G. D. Doolen, R. H. Kraichnan, and Z.-S. She, "On statistical correlations between velocity increments and locally averaged dissipation in homogeneous turbulence," Phys. Fluids A 5, 458 (1993). 31. R. M. Kerr, "Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence," J. Fluid Mech. 153, 31 (1985). 32. V. J. Peridier, F. T. Smith, and J. D. A. Walker, "Vortex-induced boundary-layer separation. Part 1. The unsteady limit problem Re-, oo," J. Fluid Mech. 232, 99 (1991). 33. P. R. Spalart, R. D. Moser, and M. M. Rogers, "Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions," J. Comput. Phys. 96, 297 (1991). 34. K. S. Breuer and R. M. Everson, "On the errors incurred calculating derivatives using Chebyshev polynomials," Preprint (1991).