What is the mystery of turbulence? One of the great unsolved problems in physics is turbulence but I'm not too clear what the mystery is. Does it mean that the Navier-Stokes equations don't have any turbulent phenomena even if we solve it computationally? Or does it mean we simply don't have a closed-form solution to turbulent phenomena?
 A: Turbulence is not one of the great unsolved problems in physics. Physics tells us exactly how turbulence emerges as a direct consequence of local mass and momentum conservation. We can create multiparticle computer models such as lattice gas automata that generate turbulence at large length and time scales. We can write down the equations that govern turbulence. These are the Navier-Stokes equations. 
From a fundamental physics perspective, turbulence is a solved problem that has entered the engineering realm a long time ago.
So what is the unsolved problem associated with turbulence? In short, turbulence is an unsolved problem not in physics but in mathematics. The point is that mathematicians struggle to answer the question if the Navier-Stokes equation always allows for solutions that at fine enough length and time scales are well behaved. In fact, if you manage to conclusively answer this question ("yes" or "no"), you will win a math prize that comes with a handsome cheque of $ 1,000,000.
In case you want to give it a try, the precise question is:

Prove or give a counter-example of the following statement: In three
  space dimensions and time, given an initial velocity field, there
  exists a vector velocity and a scalar pressure field, which are both
  smooth and globally defined, that solve the Navier–Stokes equations.

The math difficulties have to do with the fact that turbulence emerges when the highest derivative term in the Navier-Stokes equations (the viscosity term) becomes small compared to the other terms. You can take almost any non-linear partial differential equation, and get mathematicians to cringe simply by multiplying the highest derivative term with a factor $\epsilon$ and ask about the limiting behavior of the equation when $\epsilon$ approaches zero. 
Fundamental physicists shrug and continue studying real mysteries such as quantum gravity.
A: Turbulence is indeed an unsolved problem both in physics and mathematics. Whether it is the "greatest" might be argued but for lack of good metrics probably for a long time.
Why it is an unsolved problem from a mathematical point of view read Terry Tao (Fields medal) here.
Why it is an unsolved problem from a physical point of view, read Ruelle and Takens here.
The difficulty is in the fact that if you take a dissipative fluid system and begin to perturb it for example by injecting energy, its states will qualitatively change.
Over some critical value the behaviour will begin to be more and more irregular and unpredictable.
What is called turbulence are precisely those states where the flow is irregular.
However as this transition to turbulence depends on the constituents and parameters of the system and leads to very different states, there exists sofar no general physical theory of turbulence.
 Ruelle et Takens attempt to establish a general theory but their proposal is not accepted by everybody.
So in answer on exactly your questions : 
yes, solving numerically Navier Stokes leads to irregular solutions that look like turbulence
no, it is not possible to solve numerically Navier Stokes by DNS on a large enough scale with a high enough resolution to be sure that the computed numbers converge to a solution of N-S. 
A well known example of this inability is weather forecast - the scale is too large, the resolution is too low and the accuracy of the computed solution decays extremely fast.
This doesn't prevent establishing empirical formulas valid for certain fluids in a certain range of parameters at low space scales (e.g meters) - typically air or water at very high Reynolds numbers. These formulas allow f.ex to design water pumping systems but are far from explaining anything about Navier Stokes and chaotic regimes in general.
While it is known that numerical solutions of turbulence will always become inaccurate beyond a certain time, it is unknown whether the future states of a turbulent system obey a computable probability distribution. This is certainly a mystery.
A: Re: What is the mystery of turbulence?
In 1941,RA Bagnold, discussing the long-crested simple harmonic waves that arise in wind or water flows during transition – and which remain and increase in amplitude at turbulent flow rates – philosophied thus: "instead of finding chaos and disorder, the observer never fails to be amazed at a simplicity of form, an exactitude of repetition and a geometric order.” DG Thomas (Science 1964) using a layer of tiny glass beads along a cylinder in water flows found that the beads formed simple harmonic waves at transition, persisting at turbulent flow rates, relating them to Bagnold's sand waves, feeling there was a similar Physics cause. 
Both resemble the accumulation of particles in a standing waves sound field as in the Kundt's tube experiment of high school physics. Bagnold's photographs show sand particles ejected perpendicular to the flow, being deposited at shallow angles on the crests, just as one mught expect if a simple harmonic standing wave sound field had developed during transition and persisted in turbulence.
A similar wave pattern is seen, outlined by periodic dilations and narrowing of the compliant walls of arteries during rapid injections of radiopaque x-ray "dyes" during rapid injections (arteriographic standing waves), as if a stationary simple harmonic sound field was created by the injection shear forces.
In 1867, Tyndall found that specific simple harmonic sounds would cause turbulence to erupt in laminar jets, concluding that the sound waves were superimposed on similar sound waves ("vibrations") created by fluid shear forces along the tube walls, amplifying them, triggering turbulence at lower flow rates.  Tyndall believed that this solved the mystery of transition to turbulence. I agree with Tyndall.
A: Transition to turbulence is similar for air and water. Simple harmonic long crested (SHLC) shear waves develop as water, an incompressible liquid, flows during transition along a smooth flat plate.  Each LCSH oscillation-containing boundary lamina must have identical oscillations in the two laminae abutting it. Similarly, every adjacent water lamina, layer upon layer, must form similar SHLC waves. Any variation in amplitude of adjacent water laminae (increase or decrease) would cause bands of compression (impossible with liquids) or decompression (impossible without cavitation). Therefore, all boundary layer water laminae must display identical in-phase sinusoidal SHLC waves – the boundary layer flutter (BLF) shear waves of transition. 
Also, the water lamina closest to the boundary that displays these waves cannot converge on the boundary without compression, nor diverge from it without cavitation. Therefore, there must be SHLC water waves on the boundary, under the BLF wave crests (sub-BLF waves).
An oscillation (vibration) of any mass in a fluid creates a sound wave and SH oscillations in a mass of fluid, flowing along a smooth flat plate during transition, must create SH sound waves. Thus, SH oscillations (vibrations) of water flowing along a flat plate must be associated with SH sound waves, and must be reflected from the boundary transversely into the flow.  This analysis varies considerably from the accepted understanding of the fluid shear waves of transition, based on the 1941 analysis the SHLC laminar velocity oscillations found by Schubauer and Skramstad. Velocity oscillations are not shear waves at all, but are secondary effects – graphical representations of the velocities as SH laminar oscillations pass by hot wire anemometer sensors (Hamilton G, Simple Harmonics, 2015). 
By feeding on the dynamics of the flow, the SH fluid oscillations (vibrations) – and the sound waves they produce – create high energy trans-laminar transverse oscillation of molecules transmitting the sound – initially triggering spots of boundary layer transverse laminar freezing.  The focal areas of abrupt freezing shift the flow resistance to the boundary, ripping chunks out of the SHLC wave fronts, as random head-over-heels vortices (“turbulent spots”). Further increase in flow rate results in a sudden onset of established turbulence with many random turbulent spots and noise. In cylinder flows, the onset of transverse freezing of the laminae (laminar interlocking) changes the previously parabolic isovelocity profile of laminar flow to the flattened isovelocity profile of turbulence, with the resistance to flow – that bore a linear relationship to the velocity in laminar flow – now becoming related to the square of the velocity in turbulence.
When an edge juts into the boundary layer normal to the flow during transition, all nascent turbulent spots are triggered to emerge in unison along a SHLC wave front, yielding the SH sound of edge tones. In water flow in a shallow stream, a transverse linear deformity in the streambed similarly aligns all emerging turbulent spots, causing them to emerge simultaneously, creating SHLC stationary waves. 
