Turbulence could be a reason but I am not able to understand how does it is getting tripped. The inlet has a constant uniform flow. Is it numerical errors ?

Turbulence is not due to numerical error. One can describe it by saying that
when the Reynolds number becomes large, small sctuctures appear in the solution. Their size are smaller and smaller as Re is large. The Reynolds number is Re=UL/\nu, where U is the order or magnitude of the velocity, L is the typical length of the system, and \nu is the viscosity. Here there is no formal viscosity, but you can replace it by the numerical viscosity, which is proportional to the mesh size. Thus the more you refine, the more you should see some fluctuations.

At the outlet, I wanted to implement the Neumann B.C. for pressure that’s why I used the above code. Also I have imposed the no-slip boundary condition on the top, bottom, front and back (ID=9,10,11,12) surfaces of the domain along with the curved surface of cylinder (ID = 7).

I will check with a larger time step to see if the flow around the cylinder looks correct. Otherwise, I will try with a simplier 2D domain.

This is not a Neumann BC for pressure. There is a natural pressure release condition that arises from integrating the pressure gradient term by parts, but this is not a Neumann-type condition.

I think that to write no boundary term is the correct BC at outlet. It means that p=0 on this boundary, ie stress free.
On the other boundaries except the inlet, you have put the slip condition (ie no penetration u.n=0), which is good.
At the inlet you should put only u1=1 i.e. u.n=1, you are not allowed to put u2=0, u3=0.

Yes, it is correct with 1e10.
for the lateral boundaries which are flat, instead of this complicate penalty you can use the standard on( , ui=0) with ui the normal velocity.
The penalty is useful only on the curved boundary.
Also it is better to use the definition of normalapp that I gave 2 days ago.

Your code looks correct. I don’t understand well your image. Anyway it is a good way to well understand the 2d first. Note however that turbulence properties in 2d and 3d maybe very different.

where pb is a new unknown (Lagrange multiplier) and qb its associated test function. They should be taken in the same space as u, i.e. P2.
It seems to be stable with the implicit convection

Hii Everyone !
Update: I modified my governing equations. I did some reading and found about artificial compressibility method. I implemented it in 2D for the cylinder case and I got satisfactory results except at the stagnation point. I have attached the plot of pressure distribution along the cylinder surface varied with the artificial compressibility parameter and number of points on the surface.

I have used the same boundary conditions as we had discussed earlier. Dirchlet at inlet, no b.c. on outlet and rest of the walls with slip b.c. I will try increasing the domain and check whether it helps.