CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Attempting to validate force code with 2D cylinder

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   December 17, 2007, 20:33
Default Indeed, you need a much finer
  #21
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
Indeed, you need a much finer mesh. Check out this thread[1].

http://www.cfd-online.com/OpenFOAM_D...tml?1192703387
msrinath80 is offline   Reply With Quote

Old   December 17, 2007, 23:15
Default Thanks for the advice, I re-ra
  #22
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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).
andrewburns is offline   Reply With Quote

Old   December 17, 2007, 23:48
Default Here are images of the new mes
  #23
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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
andrewburns is offline   Reply With Quote

Old   December 18, 2007, 02:02
Default Nice. Very nice. You certainly
  #24
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
Nice. Very nice. You certainly do know what you're doing 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.
msrinath80 is offline   Reply With Quote

Old   December 18, 2007, 17:37
Default Thankyou. I'm going to try the
  #25
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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.
andrewburns is offline   Reply With Quote

Old   December 18, 2007, 17:55
Default Careful with the Y+ though. If
  #26
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
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).
msrinath80 is offline   Reply With Quote

Old   December 18, 2007, 18:05
Default I've just started it with k-ep
  #27
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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?
andrewburns is offline   Reply With Quote

Old   December 19, 2007, 18:39
Default Hmm I changed my mesh to be mo
  #28
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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.
andrewburns is offline   Reply With Quote

Old   December 19, 2007, 18:43
Default Yup your boundaries are messed
  #29
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
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.
msrinath80 is offline   Reply With Quote

Old   December 19, 2007, 19:02
Default Thanks for the suggestion, I d
  #30
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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.
andrewburns is offline   Reply With Quote

Old   December 20, 2007, 18:35
Default Just checking in to say my Re
  #31
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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
andrewburns is offline   Reply With Quote

Old   December 20, 2007, 19:02
Default Neat. I wonder whether this ss
  #32
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
Neat. I wonder whether this ssimpleFoam solver you refer to derives its force calculation routines from liftDrag or someone wrote it from scratch.
msrinath80 is offline   Reply With Quote

Old   December 20, 2007, 19:06
Default I think someone took the liftD
  #33
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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).
andrewburns is offline   Reply With Quote

Old   December 20, 2007, 19:08
Default How many non-orthogonal correc
  #34
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
How many non-orthogonal correctors do you have? Also what does time-step continuity error show?
msrinath80 is offline   Reply With Quote

Old   December 20, 2007, 19:38
Default Currently (running stable) I h
  #35
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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).
andrewburns is offline   Reply With Quote

Old   December 20, 2007, 19:48
Default Alright. As long as the boundi
  #36
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
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;
};
msrinath80 is offline   Reply With Quote

Old   December 22, 2007, 12:45
Default Regarding ssimpleFoam, I confi
  #37
Member
 
Rosario Russo
Join Date: Mar 2009
Location: Trieste, Italy
Posts: 56
Rep Power: 8
ariorus is on a distinguished road
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
ariorus is offline   Reply With Quote

Old   January 2, 2008, 18:22
Default Ok so I'm back from holiday an
  #38
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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
andrewburns is offline   Reply With Quote

Old   January 2, 2008, 21:24
Default paraView keeps crashing when I
  #39
Member
 
Andrew Burns
Join Date: Mar 2009
Posts: 36
Rep Power: 8
andrewburns is on a distinguished road
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?
andrewburns is offline   Reply With Quote

Old   January 2, 2008, 21:27
Default Of course transient. Turbulent
  #40
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
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.
msrinath80 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
drag force in flow over a cylinder student Main CFD Forum 1 December 13, 2008 17:59
Urgently Need the code of Lift force and VM force Kai Yan Main CFD Forum 0 July 16, 2008 07:07
Code for flow over a cylinder Abhi Main CFD Forum 3 July 14, 2006 01:49
Error attempting to open X display pradeep PANTANGI CD-adapco 5 September 14, 2004 14:25
the question about validate code Bin Li Main CFD Forum 5 October 10, 2003 16:37


All times are GMT -4. The time now is 06:44.