I've recently been trying to v
 

I've recently been trying to verify the force outputs of the 'ssimpleFoam' solver I got from another thread on this forum, it's essentially the simpleFoam solver which computes and outputs forces in 3 axes to a file during the simulation.

The problem is a 0.5m diameter 2D cylinder in a flow field of 10m/sec and an nu of 0.05, am I correct in calculating a reynolds number of 100 for this situation?

Re = (10*0.5)/0.05

Taken from a pdf I downloaded off the Nasa technical report server I have a graph of experimental Cd values for infinite cylinders Vs. Re. For a cylinder of Re 100 in steady translation through a viscous fluid I should be getting a Cd of between 1.3 and 1.5 or there-abouts.

At an Re of 100 flow around the cylinder should be laminar, however I ran a laminar simulation with ssimpleFoam and ended up with a drag of only roughly 15N for this cylinder. Using the equation:

Cd = Drag/0.5*rho*v^2*S

Where:

S = 0.5m
rho = 1.225 Kg/m^3
Drag = 15N

Cd = 0.49 Which is far far too low.

I turned on turbulence with a k-epsilon model from this point and ran another 2000 iterations. I noticed a very slow oscillation in the Fx force, I've attached a plot of Fx vs iteration as an example but you can't really see the oscillation as its not really severe.

http://i76.photobucket.com/albums/j2..._iteration.png

The mesh was made in gmsh so it's unstructured. I tried to tighten up the mesh close into the surface of the cylinder and it's fairly detailed so I don't think it's a mesh problem but I can't be certain. Here are pics of the mesh as you zoom in.

http://i76.photobucket.com/albums/j2...whole_mesh.png
http://i76.photobucket.com/albums/j2...mesh_close.png
http://i76.photobucket.com/albums/j2...ltra_close.png

The average final value of drag from the turbulent model solution was about 23.5N, which using the above equation for Cd still only gives a Cd of 0.767 which is far too low.

Does anyone have any ideas as to why my Cd are so far under what they should be? Other than nu in transport properties are there any other ways to specify the fluid properties like density that I'm missing?

Andrew, What is the depth o
 

Andrew,

What is the depth of your model the z-direction?
If it is large, it can introduce error in your force estimates. This is due to a really large aspect ratio of the cells closed to your cylinder, and the procedure used to compute forces. (If you do a search for liftDrag, you will find a thread where Dr. Jasak mentions this).
Try making you z-depth something like 1% - 10% of the diameter of your cylinder.

If the solver outputs forces then your drag will be:

Cd = (Force per unit span)/(1/2 * rho V^2 * chord )=(Drag/(z-depth))/(1/2 * rho * v^2 * S).

I mean your S, which I assume is the diameter of your cylinder.

Alternatively, you can compute CD as:

CD = (Force)/(1/2 * rho V^2 * wettedArea)=(Drag)/(1/2 * rho * v^2 * (z-depth)*pi*S).

By S I mean your S.

I hope this helps,

Alessandro

Additional resources: http
 

Additional resources:

[1] http://www.cfd-online.com/OpenFOAM_D...es/1/2726.html

[2] http://www.cfd-online.com/OpenFOAM_D...es/1/2320.html

[3] http://www.cfd-online.com/OpenFOAM_D...es/1/1678.html

Andrew, I forgot to mention
 

Andrew,

I forgot to mention the following:

nu = mu / rho

so your value of viscosity already includes information about density.

This means that:

Cd = (Force per unit span)/(1/2 V^2 * chord )=(Drag/(z-depth))/(1/2 * v^2 * S)

CD = (Force)/(1/2 * V^2 * wettedArea)=(Drag)/(1/2 * v^2 * (z-depth)*pi*S).

In other words OpenFOAM forces are really Force/density.

Alessandro

The current cylinder depth is
 

The current cylinder depth is 100% cylinder diameter, so yes the aspect ratio is fairly far off close into the cylinder. However yesterday I accidently tried running with a mesh thickness 1% of the cylinder diameter and I could not get convergence, in fact the solution pressure would explode after only a few iterations, I think this is because while the prisms up close to the cylinder were fine, the ones out at the boundaries of the solution were flat as pancakes.

I will try 10% diameter depth and see how that goes however. One thing I am still confused with is how OpenFoam deals with dynamic and kinematic viscosities. Just for reference, which of the two is nu in transport properties?

Andrew, I think nu stands f
 

Andrew,

I think nu stands for kinematic viscosity.

Alessandro

Andrew, Do you enforce the
 

Andrew,

Do you enforce the following boundary conditions?

z-constant plane towards you (type empty)
z-constant plane away from you (type empty)

Since you problem is 2-D, the pancake like elements far away from the cylinder should not matter too much as the bad-aspect-ratio side is being neglected (if you set them to type empty).

Not super sure about this last statement. But make sure you make your frontBack planes as type empty.

Alessandro

Previously the Z planes toward
 

Previously the Z planes towards and away from me were set as symmetryplane.

I've now set them to empty and I reduced the height of my mesh to 10% of the cylinder diameter and I'm re-running the solution, which appears to be running correctly. If it does work and my answer is different again I may try reducing the cell heights to 1%.

Still getting a Cd that's too
 

Still getting a Cd that's too high, approximately 0.75.

Sorry I meant the Cd is too sm
 

Sorry I meant the Cd is too small, not high.

"I've recently been trying to
 

"I've recently been trying to verify the force outputs of the 'ssimpleFoam' solver I got from another thread on this forum, it's essentially the simpleFoam solver which computes and outputs forces in 3 axes to a file during the simulation. "

Which thread are you referring to?

http://www.cfd-online.com/Open
 

http://www.cfd-online.com/OpenFOAM_D...es/1/5067.html

About half way down that page.

I used the same mesh and increased Re to approx 500000 and ended up with a Cd of around 0.9, when it should be around 0.3.

I don't think the separation point is being accurately simulated, which would explain why I'm getting low Cd's for detached cases and high Cd's for attached cases.

Pack your case and post it her
 

Pack your case and post it here or email it to me. I will look into it. My email address is available in my profile.

Hi Andrew, To check force o
 

Hi Andrew,

To check force outputs, I wouldn't try the cylinder case with a Re of 100. If I remember well, this is about the starting point for vortex sheding, which might explain the oscillations you noticed.
Do you have results to compare with for other Re numbers?
Then, I'm not sure it's right but I would computed Cd as:
Cd = Drag/(0.5*v^2*(domaineDepth*D))
with D being the diameter of the cylinder and domainDepth the depth of your domain (1%, 10% of D)

Hope this helps a little,

Vincent

Why don't you start off with a
 

Why don't you start off with a static 2D cylinder at Re=200 using icoFoam. Just laminar vortex shedding which is well documented and you can test if your forces are calculated correctly.

Frank

Srinath Madhavan, I'll pack up
 

Srinath Madhavan, I'll pack up the case and send it to you tomorrow when I'm back at work, thanks for the help.

Currently I'm trying the simplefoam case because the only solver I have that outputs forces is ssimplefoam, I haven't yet compiled liftdrag into my version of openfoam.

Srinath, did you get a chance
 

Srinath, did you get a chance to look at the case I sent to you? I sent it to your gmail.

Yup, got it. Sorry I have not
 

Yup, got it. Sorry I have not had a chance to look at it. I am at work right now. I'll give it a looksie once I'm at home. By the way I noticed that your Re is far from the laminar regime. I almost certainly can tell you to expect discrepancies when comparing dimensionless force coefficients (such as lift/drag) with experiments or DNS data. As you can see, there are folks still working on validating the liftDrag tool for turbulent flows. However, the one regime where I can be confident that you will get good agreement with exp. data is the laminar steady or even unsteady regime. I am willing to try out your case for that Re if you are interested.

Also, what Frank suggests is v
 

Also, what Frank suggests is very logical. I would also recommend that you start off with low Re, validate your Cd/Cl calculation and then progress to higher Re. Frank has performed comprehensive tests on a circular cylinder[1]. Be sure to check that thread for pointers. Other things worth investigating are:

a) Use of a convective outlet B/C. See [2].
b) Use of a Fixed mean value pressure outlet B/C. See [3].


References:

[1] http://www.cfd-online.com/OpenFOAM_D...tml?1152126462
[2] http://www.cfd-online.com/OpenFOAM_D...tml?1190761096
[3] http://www.cfd-online.com/OpenFOAM_D...tml?1197751034

I just tried the cylinder with
 

I just tried the cylinder with an Re = 50, this Reynolds number should be too low for vortex shedding so a steady state solution should do, Cd should be around 1.5.

I'm getting a Cd of 0.68 after 1500 iterations, pressure has converged to 1*10^-3 and the force is fairly stable, slightly decreasing.

The force readout is:
Total pressure Force = (0.012949 7.69589e-05 6.39884e-21)
Total viscous Force = (0.00678759 -5.97422e-05 -2.00328e-21)
Total turbulent Force = (0 0 0)
Total Force = (0.0197366 1.72167e-05 4.39556e-21)

(note I'm running a laminar only solution). One thing to note is viscous force is about half that of pressure force, could it be that the mesh close to the surface of the cylinder is not fine enough to accurately capture the effects of viscous force? I thought that at such a low Re viscous forces would account for a large percentage of drag.

Indeed, you need a much finer
 

Indeed, you need a much finer mesh. Check out this thread[1].

http://www.cfd-online.com/OpenFOAM_D...tml?1192703387

Thanks for the advice, I re-ra
 

Thanks for the advice, I re-ran the simulation for Re = 50 while tightening up the mesh quite significantly around the surface of the cylinder. I ran the simulation until all the momentum residuals were around 1*10^-4 and the pressure residual was down to around 5*10^-4 (about 700 iterations) and the forces were fairly steady.

However this time the boundary cells must have been enough because the final force readout was:

Total pressure Force = (0.00936247 -0.000599356 1.53853e-22)
Total viscous Force = (0.0282414 1.3681e-05 3.44915e-22)
Total turbulent Force = (0 0 0)
Total Force = (0.0376038 -0.000585674 4.98768e-22)

As you can see the viscous force is much larger than before. Using this force I calculate a Cd of 1.504, which is close enough to the ~1.5 value I took from a graph. Additionally my mesh was only around 280000 cells as it was very coarse in the far-field (not attempting to model vorticies).

Here are images of the new mes
 

Here are images of the new mesh if you're interested. It was made with gmsh using un-structured prisms, obviously I'd much rather a hex mesh around the cylinder but I don't know if gmsh is capable of that in this case.

Far field:
http://i76.photobucket.com/albums/j2...6/new_far1.png

Close into the cylinder:
http://i76.photobucket.com/albums/j2...new_close1.png

Boundary layer:
http://i76.photobucket.com/albums/j2...w_boundary.png

Nice. Very nice. You certainly
 

Nice. Very nice. You certainly do know what you're doing http://www.cfd-online.com/OpenFOAM_D...part/happy.gif Dig around the gmsh posts in this forum, you may pick up some useful pointers for a hex boundary layer. Good Luck and keep us posted if you compare higher Re solutions with data in the literature.

Thankyou. I'm going to try the
 

Thankyou. I'm going to try the same mesh with an Re of 500000 today to see if I can validate for turbulent flow (Cd should be around 0.3). I may have to refine the boundary layer cells more however as it stands to reason the boundary layer will be much more compact for high Re flows.

Careful with the Y+ though. If
 

Careful with the Y+ though. If you use wall functions, there are limitations. If on the other hand you are running pure Navier-stokes, refinement is the way to go (assuming you are doing DNS).

I've just started it with k-ep
 

I've just started it with k-epsilon, I don't really know enough about running other turbulence models and I'm limited to those supported by the ssimplefoam solver, until I can get the liftdrag utility working.

I checked my Y+ and it's actually a lot better than I thought it was going to be, average of around 0.5 and maximum about 0.8.

Suggestions?

Hmm I changed my mesh to be mo
 

Hmm I changed my mesh to be more detailed in the entrance and wake regions in front and behind the cylinder and somehow my simulation ended up looking like this (can you spot the problem).

http://i76.photobucket.com/albums/j2..._06/whoops.png

Clearly the boundaries are set incorrectly but I don't know how, they're properly set in gmsh and named in the same order in my boundaries file, I'm at a loss as to why it's all suddenly decided to change on me.

Yup your boundaries are messed
 

Yup your boundaries are messed up. It's like a lid-driven cavity problem now with three moving walls and one outlet? Check if your boundaries are what you think they should be by opening the mesh in paraFoam and selecting each boundary individually. They way it is at present, it looks like you have three patches sharing the inlet patch.

Thanks for the suggestion, I d
 

Thanks for the suggestion, I didn't know you could isolate patches in paraview. I just checked and the top and bottom faces of the boundary are 'inlets' while the end (which is supposed to be the outlet) is symmetry plane and the inlet is the outlet.

Very strange though, I've always had the patches in the polymesh 'boundary' file in the same order as they appear in the gmsh .geo file, I don't know why they're so out of order this time.

Just checking in to say my Re
 

Just checking in to say my Re = 500000 case is coming along smoothly. I initialised the solution with potentialFoam then started it in ssimplefoam with a laminar turbulence model. After 100 iterations I switched to k-epsilon turbulent model and left it running up to about 900 iterations.

At this point the drag was about 0.45N, for the correct Cd of 0.3 I'd expect the drag of the cylinder to be 0.75N, however as you can see from the graph bellow, the drag was still increasing linearly with number of iterations. I switched the scheme to a more aggressive linear one (from upwind) and tightened some relaxation factors I had low to keep the solution from blowing up initially. I'm still having problems with epsilon blowing up with the new relaxation factors so I'll probably have to revert at least one, however I hope the solution will take less time to reach a steady force now.

Here's a graph of drag vs iterations, you can clearly see where I switched from laminar to turbulent.

http://i76.photobucket.com/albums/j2...k_converge.jpg

Neat. I wonder whether this ss
 

Neat. I wonder whether this ssimpleFoam solver you refer to derives its force calculation routines from liftDrag or someone wrote it from scratch.

I think someone took the liftD
 

I think someone took the liftDrag equations and integrated them into the solver rather them being a separate routine. I don't have the separate liftDrag utility properly installed so it can't be using that.

Epsilon is still being touchy, I don't know why it's causing me so many issues unless I tread lightly (which will take a long time to converge, oh well).

How many non-orthogonal correc
 

How many non-orthogonal correctors do you have? Also what does time-step continuity error show?

Currently (running stable) I h
 

Currently (running stable) I have 3 non-orthagonal correctors and my time step continuities are:

sum local = 1.08838e-08
global = -6.28549e-10
cumulative = 2.52733e-08

Solution looks stable using relaxation factors of 0.3 for k and epsilon and using upwind for epsilon (linear for everything else).

Alright. As long as the boundi
 

Alright. As long as the bounding epsilon messages don't keep surfacing I think it's ok. If you would like the sum local, global and cumulative to drop further, converge the pressure more tightly (i.e. in fvSolution, change the tolerance for pressure from 1e-06 to 1e-08. Also try and see if the GAMG solver improves solution time. This is the format for using the GAMG solver for the pressure in fvSolution.

p GAMG
{
tolerance 1e-08;
relTol 0;

// DICGaussSeidel is slightly more expensive than GaussSeidel
smoother DICGaussSeidel;
// smoother GaussSeidel;
nPreSweeps 0;
nPostSweeps 2;
nFinestSweeps 2;

cacheAgglomeration on;
nCellsInCoarsestLevel 100;
agglomerator faceAreaPair;
mergeLevels 1;
};

Regarding ssimpleFoam, I confi
 

Regarding ssimpleFoam, I confirm it is exactly the same as liftDrag. It is useful for me having forces written at run time.

By the way liftDrag utility is not correct in computing turbulent forces near the wall (and so ssimpleFoam neither). As Hrvoje suggests it is better to use nuEff.


So I'd tell you to change the force computation in this way:

pressureForce += gSum
(
p.boundaryField()[patchI] *p.mesh().Sf().boundaryField()[patchI]
);

viscousForce += gSum
( -nu.value()*U.boundaryField()[patchI].snGrad()
*p.mesh().magSf().boundaryField()[patchI]
);


turbForce += gSum
( -turbulence->nuEff()().boundaryField()[patchI]*
U.boundaryField()[patchI].snGrad()*
mesh.magSf().boundaryField()[patchI]
);


Info << nl << "Pressure Force = "<<pressureForce;
Info << nl << "Viscous Force = "<<viscousForce;
Info << nl << "Viscous + turbulent Force = "<<turbForce;
Info << nl << "Total Force = "<<pressureForce+ turbForce;



Ciao.

Rosario

Ok so I'm back from holiday an
 

Ok so I'm back from holiday and while I was gone I ran my Re = 500k cylinder for 6000 iterations.

The forces I came back to are:

Pressure Force = 0.681694
Viscous Force = 0.0104461
Turbulent Force = -0.0162484
Total Force = 0.675891

This is with the original liftdrag routine, I haven't attempted to make any changes to it yet.

This gives a Cd of around 0.27, which is fairly close to the 0.3 it should be which is good. The problems come when you look further though.

Firstly the Fy force, which should be almost 0 is 0.340623, which tells me there must be some large assymetry in my flow field.

Secondly even after 6000 iterations my pressure residual is still only 1.5*10-3, I'd want it at least an order of magnitude or two smaller than this, this tells me that the solution isn't really converging at this point.

Finally I made some excel plots of the cylinder forces, you can see the large instability in the Fx and Fy forces. I'm making an animation of the velocity field in paraView now but it's going to take a while, I'll upload it as soon as it's done.

http://i76.photobucket.com/albums/j2...k_forces_1.png

paraView keeps crashing when I
 

paraView keeps crashing when I try to make an animation of my solution, so here are two pics instead.

The first is the velocity field around the cylinder, the second are vector glyphs of the velocity field.

http://i76.photobucket.com/albums/j2...h_re_cyl_1.png
http://i76.photobucket.com/albums/j2...h_re_cyl_2.png

I'm thinking at this point my lift/drag values are not converging because I'm trying to run a steady state solution on a flow which clearly has large detached vortices, do you think I should run a transient solution instead?

Of course transient. Turbulent
 

Of course transient. Turbulent vortex shedding is clearly evident. And please don't even try to expect a good agreement with *anything* when using turbulence models. Your best bet is using LES.