Vortex shedding is a natural phoenomenon we can experiment constantly in our life. The symmetrical solution appears at very low Re number but it is not stable when Re is greater than a threshold. Such kind of bifurcation is very typical of the non-linear dynamics, that is when the convective part of the momentum transport has a magnitude order greater than the diffusion counterpart. Vortex shedding is still laminar and deterministic up to a certain Re number, then becomes turbulent and a wider range of characteristic scales appear.

A numerical solution, wherein the scheme has no relevant artificial viscosity and the grid is fine to solve correctely the boundary layer, simply reproduces the instabily of the flow problem.

Note that using half domain and symmetric conditions forces the solution to be symmetric, a fact that is not physical, at least if you do not solve the statistically averaged NS equations.

The question is indeed do you have true symmetry, or are numerical errors the cause. I've never tried, but I suspect that when resolving half the domain and using a symmetry boundary condition at the centerline, you won't see vortices. We did a similar run with 2-phase microfluidics once, about trapping a bubble in a diverging flow, in a T-junction. In case the junction was fully modeled, the droplet would never get stuck at the impingement point due to numerical differences which would explode; it would always drift to one side. Only when modeling half the domain with a symmetry plane in the center, the droplet would be truly stuck. (edit: I agree with reasons below about onset of vortex shedding, the above was merely speculation about what might happen if you enforce the flow to be symmetric)

Dear Mohammad, von Karman vortices are the results of flow separation around an obstacle for appropriate Reynolds number values; I don't think that they depends on inequality in roughness between "lower and upper surfaces". They appear also in numerical simulations, where, in general, geometry has not "roughness", because NS equations are in quite a good mathematical model of flow motion, therefore their numerical resolution can reproduce real shape.  Gianluca

Vortex shedding is a natural phoenomenon we can experiment constantly in our life. The symmetrical solution appears at very low Re number but it is not stable when Re is greater than a threshold. Such kind of bifurcation is very typical of the non-linear dynamics, that is when the convective part of the momentum transport has a magnitude order greater than the diffusion counterpart. Vortex shedding is still laminar and deterministic up to a certain Re number, then becomes turbulent and a wider range of characteristic scales appear.

A numerical solution, wherein the scheme has no relevant artificial viscosity and the grid is fine to solve correctely the boundary layer, simply reproduces the instabily of the flow problem.

Note that using half domain and symmetric conditions forces the solution to be symmetric, a fact that is not physical, at least if you do not solve the statistically averaged NS equations.

The question was why vortex shedding appears in numerical simulation. It seems that with raising the Reynolds number (inertial forces become more important with respect to the viscous one) the symmetric solution becomes less and less stable. (Imagine a circular plate clamped to the wall by a nail. If the nail is above its center of gravity the plate is in stable equilibrium. By lowering the nail position the stabilization moment is lower, and when the nail is below the center of gravity the equilibrium is labile, and any small perturbation will turn the frame in stable position).

In a real flow around a cylinder, any imperfection in: geometry, roughness, fluid homogeneity, velocity profile uniformity, external vibration and so on can be a trigger to force the flow to oscillate.

In the numerical approach the round off error is one kind of trigger. Also due to Navier-Stokes equation non-linearities, the discretized equations are non-linear and never solved exactly. Also the systems of linearized equations are usually solved by an iterative method, and the iteration procedure usually stops before the residuals are reduced to zero. All these numerical errors can initiate oscillations in the solution.

I believe that the symmetrical solution is intrinsically unstable (the higher the Reynolds number, the smaller perturbation is needed to initiate oscillations). Even if you try to force the solution symmetry at high Reynolds number (i.e. to compute the flow in one half of domain), you will see that residuals will not converge toward zero but they will oscillate around certain value. In my opinion, these oscillations clearly indicate that the steady state solution no exists.

I agree. A simple experiment about the dependecy of a symmetric solution on the Re number is possible using the simple 2D lid-driven flow. Start at Re=0 and find the steady state solution. Look at the symmetry of the flow. Then increases the Re number, 1, 10, 100, 1000. For increasing effect of the non linearity (convection) over the linearity (diffusion) terms, the solution will be steady but no longer symmetric,

The ultimate driver of vortex shedding is of course that the shedding solution  it satisfies the second thermodynamic principle better than the symmetric solution.

In all bifurcation problems there are at least two paths possible, and it is the Second Principle that must be invoked to pick the most likely to occur. . .