[TPMann]

Hello,
After many hours trying to figure this out, I'm at a loss. I will preface this with the fact that my knowledge of aerodynamics is not as solid as it should be, and I'm hoping to at least be pointed in the right direction.

I am attempting to simulate a single Selig-1223 wing element at 20 degrees AOA and 15m/s (Re>100,000) in a mesh generated in Ansys, so that I can prove similarity of results. I was able to do this easily with Spalart-Allmaras, but the kOmegaSST model is providing a far more difficult fight.

I started with the default conditions from the motorBikeSteady case and loaded my mesh in, which resulted in a brief runtime ending in a floating point error. From the iteration readout I can see that leading up to the crashes, the GAMG solver on p starts reaching 1000 iterations for several timesteps and then crashing. Attached file is the output from one example of this.

I have tried dozens of boundary conditions for k and omega, including but not limited to fixedValue, zeroGradient, and turbulentIntensityKineticEnergyInlet.

I am at a loss as to how to fix this. If anyone can offer some guidance on how to bring this back into the realm of stability I would greatly appreciate it.

[Tobermory]

Some typical approaches to improve stability:

1. You are running SIMPLE and so the fields need relaxing - check your relaxation factors and adjust if necessary.

2. Start your simulation with stable solver settings (eg upwind divSchemes, more conservative settings for laplacianSchemes and snGradSChemes)

3. Check your mesh quality, and if poor then improve!

4. Check your freestream turbulence levels in the domain (i.e. plot them before the crash) - have they decayed to nothing? Turbulence adds diffusion, which can add stability to the solution.

Good luck!

[TPMann]

1. You are running SIMPLE and so the fields need relaxing - check your relaxation factors and adjust if necessary.

Following this, I attempted to change relaxation factors. Additionally, I tried defining a relaxation factor for p in case it needed one, but it either didn't care of didn't help.
U: 0.9 -> 0.7
k: 0.5 -> 0.3
omega: 0.5 -> 0.3
p: Does Not Exist -> 0.3

2. Start your simulation with stable solver settings (eg upwind divSchemes, more conservative settings for laplacianSchemes and snGradSChemes)

Attempted the following changes to divSchemes, received a Foam Fatal Error for div(phi,U) (attempting to read past EOF, seems to happen when it doesn't like the spacing between name and definition). Restored to default to see if that was the only problem, runtime returned to crashing after about two dozen timesteps.
default none;
div(phi,U) bounded Gauss linearUpwindV grad(U); -> Gauss linearUpwind
div(phi,k) bounded Gauss upwind; -> Gauss upwind
div(phi,omega) bounded Gauss upwind; -> Gauss upwind
div((nuEff*dev2(T(grad(U))))) Gauss linear;

3. Check your mesh quality, and if poor then improve!

This exact mesh ran without issues with SpalartAllmaras on OpenFOAM and in Ansys with k-omega SST. I tried to upload the mesh to this reply but got a missing security token error(?). In a previous setup that managed to survive a short time, I noticed strange fluctuations in omega and nut that appear near regions of mesh that were split in the geometry to set controls for bias, but in p the boundary at the wing appears noisy, with sharp jumps in value between neighboring cells.

4. Check your freestream turbulence levels in the domain (i.e. plot them before the crash) - have they decayed to nothing? Turbulence adds diffusion, which can add stability to the solution.

The conditions in /case/0 have internalField definitions, but no freestream defined. I'm not sure if I'm interpreting this instruction correctly, so I apologize if it's not what was intended.

[Tobermory]

For relaxation factors, beware that there are two sets - one for fields and one for equations. If you are running SIMPLE (ie have "consistent no;" in the SIMPLE{} block in fvSolution) then you will need something like:

relaxationFactors
{
   fields
   {
      p   0.3;
   }
   equations
   {
      U   0.7;
     "(k|epsilon|omega)"   0.7;
   }

you can even drop the p factor to 0.2, to slow down the changes to the solution. If you are running SIMPLEC ("consistent yes;") then you don't need to under-relax p and something like the following should work as a starting point:

relaxationFactors
{
   equations
   {
      U   0.9;
     "(k|epsilon|omega)"   0.9;
   }

On the divSchemes, try "bounded Gauss upwind;" on U, k & omega (the "bounded" is optional, but is relevant because you are running a steady solver). The 1st order scheme (linear uipwind is 2nd order) is more stable.

On the mesh - what does checkMesh report? What's the max aspect ratio and the max non-orthogonality? If the a/r is too high, then the solver can struggle; if the orthogonality is bad then you'll need some non-orthogonal correctors to manage it.

Finally, for "freestream" - yes, I meant the values out in the main domain, ie. the internalField values as you noted. Check the values (k, omega and nut) at the inlet are sensible, that the starting field (t=0) is sensible, and that they convect through the domain properly (t>0) without nut dying away to nothing.

I hope that the above help you get closer. Otherwise, trying uploading the case and maybe I will be able to find some time to get it running.

[TPMann]

you can even drop the p factor to 0.2, to slow down the changes to the solution.

Relaxing this seemed to stabilize everything somewhat, everything except p maintained 1-2 iterations while p bounced from 70-1000 iterations over several timesteps. This could have also been related to the switch to bounded Gauss upwind. It was unfortunately not stable enough to get past ~70s still, which means I've likely missed something else important.

On the mesh - what does checkMesh report? What's the max aspect ratio and the max non-orthogonality?

The checkMesh output is attached to this reply. The only alarm that appears is the # severely non-orthogonal faces, of which there are 32.

I hope that the above help you get closer. Otherwise, trying uploading the case and maybe I will be able to find some time to get it running.

I cannot overstate how much this helps, I am beyond appreciative that you have taken the time to help diagnose this. I apologize for my misunderstandings, I'm drinking from one firehose of server management and another of aerodynamics/fluid simulation lately and the fluid dynamics has been the harder of the two.

Due to size, I dumped the entire case into DropBox so that it can be accessed whenever needed. Please let me know if there are additional files or information needed, I can try to scrape information out of Ansys if more is needed.

https://www.dropbox.com/scl/fi/bu4ry...=0a7d5v6n&dl=0

[Tobermory]

Okay - I think that we are getting somewhere. OpenFOAM CAN be overwhelming at the beginning, especially since it is much less forgiving than the main commercial codes that I spent decades using. On a positive note, there is lots of good training material out there (eg Joel Guerrero at Wolf Dynamics - I strongly suggest you google him; he has a YT channel as well), but you have to actively search for it.

The non-orthogonality is perhaps an issue - you should add in some correctors for this (change nNonOrthogonalCorrectors in the SIMPLE block in fvSolution to 3, for example.

I will see if I have some time later on today to take a look at your dropbox case ... incidentally, which version of OF are you running?

[TPMann]

OpenFOAM CAN be overwhelming at the beginning

I'm enjoying some of it so far, the sheer volume of useful tools is truly amazing.

Joel Guerrero at Wolf Dynamics

I actually stumbled across some of his powerpoints online, they've been a guiding light for many of the previous issues I had trouble with before.

The non-orthogonality is perhaps an issue

The default from the motorBikeSteady case is 10, I've tried replacing this value with 0,3,4,5,8, and 10 with only a slight appearance of improved stability around 4, although this could have been any number of changes I tried

which version of OF are you running?

My systems are running OpenFOAM dev (build dev-d118afa2bd5b) and openmpi 4.1.2 on Ubuntu 22.04.5. Please let me know if any hardware/hypervisor information is needed for troubleshooting.

[Tobermory]

I have had a play with your test case - this is the thinking pattern that I applied.

I started by checking the fvSchemes/fvSolution set up, and changed a few things:

• set consistent to no (ie run SIMPLE not SIMPLEC)
• added 1 nonOrthogonalCorrector
• adjusted the relaxation factors slightly
• changed relTol on pressure to 0.1 (i.e. work less hard on each iteration, and just do more iterations, to speed things up)

The simulation still blew up (continuity exploding). I then checked the boundaries. Your setup was probably okay, but I used the tutorial example airFoil2D to replace your inlet/outlet with freestream boundaries ... which essentially do the same thing. I also adjusted the inlet values to "more sensible" values. Had no effect; still crashed.

I then plotted the fields in paraFoam (having adjusted controlDict to save every 5 iters) and the nut field looked very unhappy - a great big hole around the airfoil. The omega values are also enormous around the airfoil.

So I extracted the yPlus values using foamPostProcess and this was illluminating: for time t=0, the y+ values on the airfoil vary from 2.5e-4 to 5.4e-4 ... i.e. are REALLY small - too small in fact. These should be order of 1. I imagine that this is the reason why the omega field goes crazy and the nut field collapses?

So, in conclusion - please check your setup. The near-wall mesh seems to be exceedingly fine - is that intended? Seems too fine for the viscosity and velocity that you are running with. I don't imagine that you'll able to stabilise the kOmegaSST solver with such low y+ values, whatever you try.

I hope that helps ...

[Tobermory]

Just as an addendum, I increased the viscosity by a factor of 1000 (to bump up the y+ values), and the case seemed to run okay, which indicates that this is indeed the cause of your solution blowing up.

PS. you can run for 1000 iterations with an inflated viscosity to get the simulation up and running, and then restart with the correct viscosity value. That worked for me (with low Re wall boundaries - see my next post).

[TPMann]

the y+ values on the airfoil vary from 2.5e-4 to 5.4e-4 ... i.e. are REALLY small - too small in fact.

The intent here was to get a fine enough mesh at the boundary to get a velocity gradient modeled correctly. If I recall correctly, the bias factor near the wing is somewhere close to 90,000 in an attempt to get the refinement needed.

viscosity by a factor of 1000 (to bump up the y+ values), and the case seemed to run okay

I would not have expected that. Maybe something happened in the fluentMeshToFoam process that messed up a scale factor somewhere? I'm confused as to how Ansys was able to run this exact mesh without issues. It didn't converge very nicely due to the stall conditions that are expected at 20 degrees, but it certainly didn't blow up.

Is there something either in fluentMeshToFoam or in the way the simulation is run that would cause such different behavior in OpenFOAM?

[Tobermory]

... and as an addendum to the addendum - it just occurred to me that you are using basic wall functions: nutkWallFunction on nut, omegaWallFunction on omega, kqRWallFunction on k.

Remember that for these to work properly, you need to ensure that the yplus is somewhere between 30 and 300. For resolved boundary layer simulations, dump the wall functions and use zeroGradient for k & nut and fixedValue of 1e-10 for omega at the wall. Running like this, and using the above viscosity hack, I got a plausible simulation with an average y+ of around 0.1 which is still a little low.

[Alczem]

Originally Posted by TPMann

I would not have expected that. Maybe something happened in the fluentMeshToFoam process that messed up a scale factor somewhere? I'm confused as to how Ansys was able to run this exact mesh without issues. It didn't converge very nicely due to the stall conditions that are expected at 20 degrees, but it certainly didn't blow up.

Is there something either in fluentMeshToFoam or in the way the simulation is run that would cause such different behavior in OpenFOAM?

Hey,


Since you brought it up first, I will mention that I have to scale down the meshes I import from Ansys, they are systemically a thousand time larger after import, so make sure to double check the bounding box dimensions reported by checkMesh, and you will probably see an issue. Although the fact that the y+ is so small makes me think you will have to scale up perhaps?

[Tobermory]

Originally Posted by Alczem

Hey,
Since you brought it up first, I will mention that I have to scale down the meshes I import from Ansys, they are systemically a thousand time larger after import, so make sure to double check the bounding box dimensions reported by checkMesh, and you will probably see an issue. Although the fact that the y+ is so small makes me think you will have to scale up perhaps?

Yes - good idea, I had the same thought, but the wing chord seems sensible in Paraview (~1m), so I don't think that's it.

It looks like getting the simulation started in this stalled condition is the challenge - bumping up the viscosity for the first 500 iterations seems to help stabilise things. Although my simulation did eventually go nuts - it seems to really struggle at the trailing edge. Looking at the mesh - it is pretty bad there, with a sudden transition between slanted prism cells aligned with the airfoil to orthogonal cells. As ever, improving the mesh always helps improve simulation stability.

Finally, thinking more about the flow case (stalled wing, strong separation off the back face, turbulent shedding etc.), then I think that the best approach for you might well be to model it as a transient case - especially with such a fine grid, where you will be resolving much of the turbulence shed from the wing. Transient solvers are inherently more stable than steady (no need for under-relaxation) ... BUT typically take much longer to run, and you will need to time average to get the mean lift, mean drag etc., so it's not a simple/easy suggestion. One to consider though if you can't get the steady approach to work.

[arjun]

Originally Posted by TPMann

I'm confused as to how Ansys was able to run this exact mesh without issues.

Because they are not openfoam (or any other solver). The point is that they are not doing exactly same thing when the mesh is non perfect.

The same thing would run in Wildkatze also just fine. Now we also do lots of things in different way than openfoam. For example, just the gradient procedure is at least two times more expensive than openfoam and makes sure that things work even if mesh has very bad cell quality.

Per iteration cost then expected to be more than openfoam. So on perfect mesh OF might seem faster too.




PS: In recent tests though, it seems wildkatze is 20 percent faster than OF per iteration.

[Alczem]

Commercial softwares seem to tune the numerics on the fly depending on the mesh quality of course. In OpenFoam, it is up to the user, but there are many options to consider (although I have never managed to get a faster solution with OpenFoam, iteration wise, but I am not as knowledgeable as others).


Regarding the gradient computation and its ability to deal with poor mesh quality, the default is usually Gauss linear, but I have seen iterativeGauss linear being used with relative success, and leastSquares as well. That might help.


Using the corrected scheme for laplacian and snGradSchemes is also quite ambitious considering the max non orthogonality of 85. I would go with uncorrected first, run the case, then switch to limited 0.2 and increase up to 0.5 if the case is stable. But I second Tobermory about the mesh quality being the most sensible factor to improve.

[TPMann]

Thank you all for your help. I am starting work on making fresh meshes, although the Ansys meshing algorithm seems hellbent on making very open, rough meshes with many non-tets and throwing errors while doing so, but that's a problem for a different thread.

I had hoped that my boundary conditions were the main problem and would be a silver bullet to straighten out these crashing simulations, but as it seems with most CFD the mesh is the root of the problems here. I'll do some digging to see if I can find any good sources on how to avoid bad mesh features, and from there I should be on a better track.

Again I greatly appreciate your help in figuring this out.