Greetings again,
I'm simulating the flow past a I-beam using k-omega SST turbulence model. I could get the solutions for grids with max y+ of 60, 50 and 40. However, when I simulate with y+ 30 or below, the calculations just blow up (I've tested with and without wall functions).
I've already set the number of pressure correctors up to 3 (PISO) and number of outer correctors up to 2 (PIMPLE), used upwind schemes for the convective terms, enlarged the computational domain, and the results were the same. (I'm using Crank-Nicolson 0.5 for the time derivatives, Courant number = 1, I'm running in parallel with 4 processors).
I've atteched some figures:
1, 2 - The mesh;
3 - Pressure field for time = 4.4 s;
4 - Initial residuals for the y+ = 30 grid;
5 - Drag coefficient for the y+ = 30 grid;
6 - Drag coefficient for the y+ = 40 grid.
Does anybody have any clue what could be causing this?
Thanks in advance.
Dear Renan, in numerical simulations of fluid flows a blow-up of solutions vs grid increasing usually (but not ever) is due to some incompatible boundary condition. What are the BCs which you have used? Gianluca
Gianluca, thanks for the replie.
For velocity (U):
internalField uniform (1 0 0);
boundaryField
{
Inlet
{
type uniformFixedValue;
uniformValue constant (1 0 0);
value uniform (1 0 0);
}
Outlet
{
type zeroGradient;
}
Free-slip-wall
{
type slip;
}
Beam
{
type uniformFixedValue;
uniformValue constant (0 0 0);
value uniform (0 0 0);
}
}
For pressure/rho (p):
internalField uniform 0;
boundaryField
{
Inlet
{
type zeroGradient;
}
Outlet
{
type fixedValue;
value uniform 0;
}
Free-slip-wall
{
type zeroGradient;
}
Beam
{
type zeroGradient;
}
}
For k:
internalField uniform 0.00015;
boundaryField
{
Inlet
{
type fixedValue;
value uniform 0.00015;
}
Outlet
{
type inletOutlet;
inletValue uniform 0.00015;
value uniform 0.00015;
}
Free-slip-wall
{
type fixedValue;
value uniform 0.00015;
}
Beam
{
type fixedValue;
value uniform 1e-10;
}
}
For omega:
internalField uniform 7.5;
boundaryField
{
Inlet
{
type fixedValue;
value uniform 7.5;
}
Outlet
{
type inletOutlet;
inletValue uniform 7.5;
value uniform 7.5;
}
Free-slip-wall
{
type fixedValue;
value uniform 7.5;
}
Beam
{
type fixedValue;
value uniform 1e-10;
}
}
For nut:
internalField uniform 0;
boundaryField
{
Inlet
{
type calculated;
value uniform 0;
}
Outlet
{
type calculated;
value uniform 0;
}
Free-slip-wall
{
type calculated;
value uniform 0;
}
Beam
{
type nutUSpaldingWallFunction;
value uniform 0;
}
}
The inletOutlet means that if the velocity vector at the outlet points out of the domain, then the boundary condition will be of the Neumann type, ie, zeroGradient; and if the velocity vector points into the domain, then we can imagine that the outflow is no longer an outflow, but an inflow, and that we therefore would like to specify a Dirichlet boundary condition, ie, with constant value.
I'm not very sure if the inletOutlet condition for k and omega is consistent, maybe it could be the source of the problem.
Dear Renan Thomes,
Hello again. I don't have much experience with using the k-w SST model.
But I do have some general advice. Since I am into algorithm development, tracking solution blowup is a simpler matter in the sense that it is possible to interrogate the code about two pieces of information:
(a) from where does the blowup initiate?
(b) which variable (U,V,p,k,w) is involved?
I know that this is slightly difficult to manage in FLUENT...however, if you could possibly generate some surface monitors in "vulnerable regions" of the domain (probably near-wall region of the beam) that may be prone to divergence....I think you may be able to narrow down the problem.
One more thing...(you probably might have used this but still) I am attaching a link to a y+ calculator on CFD online...
https://www.cfd-online.com/Tools/yplus.php
With few trials on the calculator....one could figure out that the "estimated wall distance" decreases with decreasing "desired y+". In that case, if Δt (time step size) is kept constant whilst decreasing "desired y+"....the solution may diverge...the Crank-Nicholson scheme is not immune to divergence when flow is involved...this link may also prove useful...
https://en.wikipedia.org/wiki/Crank-Nicolson_method
You may consider lowering Δt for lower y+ values....I reckon that with the same (large) Δt, the simulation will not diverge immediately...but would do so after the flow separates from the I-beam....that's when larger velocities and shear would be induced on the beam.
I hope this is helpful.
Best regards.
Shaswat Saincher
Hello again, Shaswat Anurag Saincher,
Thanks a lot for your replies.
I don't know if I completely understood your questions (especially the "where"), but...
(b) from the residuals, it seems to me that the first variable to show some problem is the p.
(a) the first thing I noticed to be wrong was the unphysical pressure field in the region around the downstream finer grid (figure 3).
I didn't understand much the "surface monitors", is it something like monitoring some variable on some surface?
In addition to that, I'm writing the fields every 0.1 s (flow time), but for some reason OpenFOAM didn't write some time steps, which seemed very strange to me.
Concerning the delta t, the algorithm guarantees that Courant number will be 1, so it automatically changes the time step size to ensure that.
By the way, how do you type the delta as a Greek letter?
Some skills possibly helpful.
1. Try to use another type of mesh, like O-mesh or unstructured mesh.
2. Significantly reduce your time step, starting the simulating with CFL=0.1 or something like that.
As I experienced, the problem always come from the mesh. For your case, mostly from the mesh around the charp corners.
Hi there,
You are not using wall functions for omega and k. Are you using yplus = 30? If yes, you have to use wall functions, but you have re-circulation and it is a good idea using automatic wall functions to deal with y+ lower than 30.
The only function implemented in Openfoam(vanilla) for omega is already an automatic wall function. When dealing with SST it is a good idea to use this omega wall function even using low-Re model because it is incorrect to set omega as zero in near wall. Since Openfoam 2.x.x there is an automatic wall functions for k as well. I do not know what is the implementation of the k automatic wall function in Openfoam but at least it is a automatic wall function and can deal in re-circulation zones(lower y+). I hope so.
Good luck
Dear Renan Thomes,
Regarding the surface monitors, I wrote in context to FLUENT....I forgot that you were using OpenFOAM....not sure whether such an option is available in OpenFOAM or not.
Regarding the Courant number, I support Xinliang Tian's suggestion of starting with a lower value such as 0.1. It is also true that divergence problems come from the mesh....however, since it is a BLOWUP in your case....I am still concerned about what exactly is causing the solution to blowup.
I am also not sure whether wall functions lead to solution blowup in later stages (?!?).
I would once again suggest that (if possible) try interrogating the solution for the exact place and variable from which divergence initiates and then try figuring out a way to resolve the issue. I know this suggestion is of a very general nature but it is due to my lack of experience with implementing wall functions in a RANS framework.
About the delta symbol....I just copy it from google...Δ....:) :)) :)))
Best regards.
Shaswat Saincher
Hi Shaswat Anurag Saincher,
"I am also not sure whether wall functions lead to solution blowup in later stages (?!?)."
I agree with you. The residual monitors are showing that the turbulence quantities are not reason of the divergence. I was only correcting the turbulence wall bc. I did not mentioned it. My fault.
Any way, i recommend to change the wall bcs to:
For velocity (U):
{
Beam
{
type fixedValue;
value uniform (0 0 0);
}
}
For pressure/rho (p):
{
Beam
{
type zeroGradient;
}
}
For k:
{
Beam
{
type kLowReWallFunction;
}
}
For omega:
{
Beam
{
type omegaWallFunction;
}
}
For nut:
{
Beam
{
type nutUSpaldingWallFunction;
value uniform 0;
}
}
It is strange to use uniformFixedValue in the wall but probably works.
Thank you all for your replies, they're being really helpful.
Xinliang Tian,
I'm restricted to use quad/hexahedral mesh, and I don't understand how a O-grid could help in this case. Could you explain better what advantages it could bring?
I'm simulating with CFL = 0.3 right now, if I get good results I'll let you know.
I do believe to have some problem with the sharp edges (see attached figure 1), but wouldn't it be worse with coarser meshes (y+ 40, 50, 60 ...)? I didn't find a way to overcome this grid problem.
I. A. Accordi,
I've tested the simulation with wall functions too, and the calculations still blew up. I used the same BC that you mentioned, except for k, that I used kqRWallFunction. However, from my researches, I understood that wall functions were not applicable for re-circulations zones. Now, I'm really confused. Why do those wall functions are suitable for the re-circulation?
Shaswat Anurag Saincher,
I'm really concerned about the pressure field in the region across the "flange line" (attached figure 2), that's the first thing to indicate the divergence, I guess. I'll try to use some relaxation factors.
Hi Renan,
I am not sure about the mathematical formulations of wall functions used in openFOAM since I work mostly with Fluent. But, normally speaking, SST k- omega models work well and accurate if is is resolving the viscous sub-layer. Now to capture this layer you should maintain the y+ value below 5.
Also, as it is can been seen your residuals shoots up after certain iterations. I would like to suggest one more area which you can look at. When the number of cells increases, even though the accuracy of the results increases, there is a possibility of longer wave length errors diminishing at a lower pace. This is applicable if you are applying an iterative solver.
To compensate for this, normally solvers use multigrid methods, say for example Algebraic Multigrid method. As far as I know openFOAM should also be using some multigrid methods. Please check the fvsolution file for the type of solver you are using for different variables.
Thanks
Dhyan
I fully agree with Shankar's reply. When this happens I usually change the approach to a less stiff approach, for example transient. It helps a lot!!. What usually happen is that as you reduce the mesh you reduce the numerical dissipation, that contrary to what people believe they usually help you to stabilize the solution. It makes sense because is dissipation at the end of the day. Therefore as you reduce the mesh more frequencies are involved and therefore more fluctuations. In RANS I'm not that sure how this affect. Nevertheless your solution must adecualtely capture the viscous stresses when you use the law of the wall so you better check the y plus for your model.
Thanks for the replies.
Dyhan, I am actually using Generalised geometric-algebraic multi-grid for solving the pressure field :)
Good news, I've managed to get soluitons in two different ways:
1. I set 3 outer correctors (PIMPLE) (atteched figure 1). However, the computational time increased so much. (I changed from upwind to linear upwind scheme at 80 s, that's why the CD changed like that)
2. I increased the number of cells (reduced the size) downstream and at the top and bottom by reducing the growth rate from 1.05 to 1.025. See figure 2 for the grid and 3 for the CD. But, I've noticed a very strange line in the pressure field (figure 4).
Would you like to make some remarks? Thanks.
- Badges
- Science topic