Have you test the size of computational domain?  Perhaps your computational domain is too small

Dear Renan Thomes,

I have gone through your pictures on IMGUR and really appreciate the effort that you have put into specifically explaining your problem. 

I do not have any experience with simulating flows around spheres but have validated CD and CL values for circular cylinders using FLUENT.

I think there is some confusion regarding the transition Re for flow over a sphere. Is it the case that the flow is susceptible to transition for Re>50? In that case, I would probably recommend an unsteady simulation (at least) for Re=80.

About your query regarding the wake inducing higher pressure at the back of the sphere...that pressure is only higher relative to the suction induced at the middle of the sphere (due to the Bernoulli effect occurring from a reduction in flow area). It is not the wake that is inducing higher pressure...wakes are generally (localized) low pressure regions in a flow. The higher pressure is due to the Bernoulli effect (yet again). As the fluid "leaves" the sphere, the area available for "flowing" increases which locally induces a "diffuser-like" effect reducing the flow velocity. By Bernoulli's principle, the pressure has to increase at this point which is why you observe relatively higher pressure at the back-side of the sphere. 

About the mesh, I think the level of concentration within the central "box" holding the sphere is sufficient....however, in my experience, the mesh should be coarsened very gradually away from this box in all directions. Hence, the cell size on either end of the "box" boundary should be almost equal and then very gradually coarsened away from it. I am attaching a picture of a 5.9 million hex mesh that I had generated for the 3D cylinder problem.

But then, I may be wrong about the mesh :) because you do seem to have gradually coarsened the mesh behind the sphere :) (?!?) And I think for the boundary layer, I had placed the first grid point about 0.1 mm away from the surface of the cylinder and gradually increased cell size away from the cylinder surface over 20 layers. 

In your simulations you report errors in both pressure as well as viscous drag, however, the error in viscous drag seems to be greater. And since your CD values are being over-predicted (especially at lower Re), could it be due to some form of numerical diffusion in the solution? I think you are using a SOU-type (linear) scheme for momentum (velocity) advection...you could try using a higher-order (quadratic) scheme for estimating U,V,W at the faces (?!?). 

Further, if you are able to query the viscous FORCE on the sphere from OPENFOAM, you could try comparing your numerical value with Stokes' Law (at low Re). This link might be of some help... https://en.wikipedia.org/wiki/Stokes%27_law ...then, you can be certain that the viscous force is being over-estimated in the simulation. 

I do not believe my answer is of much help here...however, if you do find a solution to the problem do let me know...I inquire just out of interest :)

Best regards.

Shaswat Saincher

Shiuh-Hwa Shyu,

Thanks for the reply. I don't think the domain is small since it seems like the velocity gradient far away from the sphere is approximately 0. But maybe...

Shaswat Anurag Saincher,

Thank you so much for the prompt and explanatory reply.

I'm reading an article (available at www.dtic.mil/get-tr-doc/pdf?AD=ADA494935) about the simulation of the flow past a sphere for several Reynolds numbers (20 until  10^6) and the CD results showed good agreement with the experimental data (especially for laminar flow). The article states that the regime of the flow is: Steady axisymmetric for 20 < Re < 210; Steady planar-symmetric for 210 < Re < 270; Unsteady planar-symmetric for 270 < Re < 400; etc...

About the mesh, I've thought of that too, but I didn't do that because it would increase considerably the number of nodes xD. I thought my mesh had too much nodes for a simple case (600 k), but since you said that yours had 5.9 million nodes(?) I'm starting to think mine is kind of coarse. I will try that.

" because you do seem to have gradually coarsened the mesh behind the sphere :) (?!?)" Yes, I did :) My first grid point is 0.35 mm away from the surface of the sphere (a 100 mm-diameter sphere).

"I think you are using a SOU-type (linear) scheme for momentum (velocity) advection." I think that's what OpenFOAM calls divSchemes, but I'm not sure. I'm using the upwind.

As I understood from Stoke's Law, there would be only viscous drag for very low Reynolds number. But I'm fiding both viscous and pressure drag in my simulations.

I compared the simulation's with the analytical's velocities beginning at the top of the sphere and ending at the boundary of the domain (for Re = 0.05) (see attached pictures). The error is around 20%, just like the CD error. And the error of the velocity u (x axis) of the first node is 21.4%. Perhaps, the shearStress (=viscosity*du/dz) is following this error.

Observations:

• The freestream velocity is 1 m/s and the velocity far away from the sphere is 1.03 m/s. Is that normal?
• I guess the analytical solution is not 100% correct, because the velocity will only be equal as the freestream velocity when z -> infinity;
• u = velocity in the x axis (freestream direction)

w = velocity in the z axis.

• You can find the equations I used in this link https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-090-introduction-to-fluid-motions-sediment-transport-and-current-generated-sedimentary-structures-fall-2006/course-textbook/ch3.pdf

https://ocw.mit.edu/courses/earth-atmospheric-and-planetary-sciences/12-090-introduction-to-fluid-motions-sediment-transport-and-current-generated-sedimentary-structures-fall-2006/course-textbook/ch3.pdf

https://www.dtic.mil/get-tr-doc/pdf?AD=ADA494935

Have you compare two different sizes of domain to validate the independency

https://www.youtube.com/watch?v=ylZQXiEdYYo&t=94s

Shiuh-Hwa Shyu,

I increased 20% of the domain size in all directions (20% for  -X, X ...), and the CD results changed 0.41% :)

Shaswat Anurag Saincher,

I've gradually coarsened the mesh away form the box in all directions (just like I did behind the sphere) and the CD results changed 0.10%.

As you said, I was using a first order interpolation scheme for the advection term (Gauss upwind). So I did some research in the OpenFOAM User Guide and chose the Gauss linearUpwind, which is a second order scheme. And that was the problem indeed! After I did that I'm having a 1,3%-error for 20 < Re < 200. Unfortunately, for very low Reynolds numbers (Re = 1 and 0.5) the results are still the same.

Thank you all very much for your help.

Dear Renan Thomes,

I am happy that (atleast) some of my advice proved useful. From your (I must say pretty extensive) searches and simulation attempts, we are now confident that it is NOT the rate of mesh/grid extension/expansion/coarsening that is the problem. Further, you stated that you had already tried refining the mesh which didn't lead to major changes in the results. So at this stage we can (probably) assume that it is NOT the grid that's the problem.

My experience in CFD tells me that ONE computational framework (or solver setting) cannot conquer all cases of a given problem....in your case....different values of Re of flow over the sphere.

One has to understand the balance between inertia and viscous forces in a given simulation (hoping that those are the only two things we need to worry about :))...my Ph.D. senior taught me that (♡). Hence, for 20 < Re < 200, the flow is comparatively more advection (inertia) dominant than when Re < 20. Hence, there is a risk of artificial/numerical diffusion contaminating the solution at larger Re...so, a higher order (LinearUpwind in your case) scheme is necessary. On the other hand, because the low Re flow would be diffusion (viscous) dominant, a (pure) first order scheme would probably suffice and using a higher order scheme there may involve a risk of over-predicting velocities (which higher-order schemes do...by the way). 

Regarding the free-stream velocity being 1.03 m/s far away from the sphere...the free-stream velocity will increase AROUND THE SPHERE due to flow separation since the domain has to accommodate the same flow-rate (across any section) regardless of separation bubbles getting formed. But far away from the sphere...a 3% increase in free-stream velocity might be of some concern...hope you are not using WALL/SLIP/NO-PENETRATION conditions at the domain boundaries....What you can do is carry out a flux-balance between the inlet and the outlet (an option available in FLUENT but not sure about OpenFOAM)....the inlet-outlet mass flow rate difference should be around 1e-15 kg/s (or m3/s) or less to ensure that there is no volume dilatation in the domain. 

Regarding your attempts towards validating the analytical solutions...I think the expressions are in terms of radial (Ur) and tangential (Uθ) velocities...so I hope you used some kind of transformation to convert them from curvilinear (r,θ,z) to Cartesian space (x,y,z). Regarding the differences between simulation and analytical values....the analytical solution neglects all acceleration effects occurring during interaction with the sphere (as stated in the MIT-OCW document). For instance, analytically, the vertical velocity (W) is consistently zero above the sphere which is not possible in practice (observable in your results). Moreover, the U-velocity above the sphere ought to be more in practice (also observable from your results) because of the "crowding effect" (also mentioned in the MIT-OCW document).

But the fact that the free-stream velocity will actually decrease for creeping flows (Re

You said the low Reynolds no is still bad. But how bad is it

Shiuh-Hwa Shyu, the CD error is about 20%.

what are your boundary conditions on the domain ?

Martin Pauthenet,

Inlet: Velocity uniform fixed value, [1 0 0] m/s

Outlet: Pressure fixed value, 0

Walls of the domain: free-slip wall

Sphere: no-slip wall

You should first check if your setups are correct, including the schemes, BCs, etc. You can also try other BCs. I read a paper (Lorena A. Barba, Reproducible and replicable CFD: it’s harder than you think) saying that correct but inappropriate BCs may cause strange results.

In fact I saw someone else having similar problems with yours. His solution is changing the computational domain to a cylindrical or spherical one (not rectangular). Another reason may be the boundary layer mesh, since the mesh there has large aspect ratio (lenth:width). You can try a mesh without the boundary layer.

Hello Renan Thomes, I am having a similar issue for Re=0.05. The theoretical CD = 24/Re = 480. But I am getting 520 when I am using a domain box with the size of 20d*20d*20d. So I am wondering whether you have fixed the issue. Your valuable comments would be greatly appreciated!!

Thanks and best regards,

Min

Hi Min Zhang

I have the same problem. My drag in Re = 1 is 26, while it should be 24 according to Stokes relation or 28.5 according to modified relations.

I used 20d*8d*d domain with 1.5e6 cells.

Did you find the reason?

Thanks

Ali