Uncontrolled blowup in simpleFoam
 

Hello everyone.

I am a new user to OpenFoam and I would like to use this tool to obtain a base flow in a axisymmetric jet flow at Reynolds 10^4.

I have started with the case detailed in : https://turbmodels.larc.nasa.gov/jetsubsonic_val.html that I refer to as ASJ. Using NASA's mesh in an incompressible setting, with significantly slower flow and a velocity driven inlet, I have obtained results for simpleFoam with a SA model that look reasonable (see step 500 attached). This case is at Re=10^5.

My case file is attached, I use the mesh available at : https://turbmodels.larc.nasa.gov/Jet..._129.p3dfmt.gz.

Copying the p3dfmt file inside the case and using ./Allrun should work.

However, when I try my own mesh made with gmsh (nozzle case) the calculation always blows up in about 10 iterations (floating point error, arbitrary velocities, there's a picture at the first calculation step attached). I have made minimal changes compared to the above case - running a diff --recursive between the two provided case directories should return 4 entries only.

• Two changes are related to the inlet and coflow boundary conditions that have a different speed.
• The other two are related to handling the mesh from gmsh properly (running gmshToFoam, replacing boundary types, and moving directly to patchSummary instead of calling autoPatch and createPatch).

Clearly the biggest difference in the two cases is moving from a NASA structured polygon mesh to a gmsh unstructured tetrahedron mesh - the geometry is not identical but very similar.

Changing mesh refinement (my mesh is less adaptative than the NASA one but finer by one order of magnitude), time step or viscosity doesn't seem to help. I cannot share my mesh file because of attachment size constraint, but I can share the .geo file.

I'm running gmsh v4.7.1 and openfoam v2112. I'm out of ideas.

Scaled mesh
 

Writing the above message there was a scale difference in the geometry between the cases. It's removed in the attached geo file, my problem remains.

Implementing a continuation method for viscosity (incrementally reducing viscosity, see Viscouscont script) and a structured mesh (see attached pictures, the .geo file is also included) allows me to get more results, but I still get a blowup before getting to interesting regions.

Worse, the flows I am getting look laminar (see the attached picture, meant to be a Re=10^3 flow right before blowup, the inlet wall is a 2D baffle) - which raises the question of why to go for a RANS model to obtain something a Newton solver could have spewed...

I am pretty disappointed with OpenFoam. It is a very opaque tool (I once had a gmshToFoam error because of excessive refinement on a boundary that was more misleading than human-readable), the documentation is spread across very shallow guides, gigabytes of tutorials and a wiki that's more designed for supporters of the software than end user and often presents a single line of explanation for the meaning of anything.

Besides, my previous posts on this forum (very old-school) have gone unanswered despite almost 80 views. I feel like OpenFoam is only open-source in name and designed for paid users - when confronted to any practical problem with it, one has to know an expert to get anything done...

Without full understanding of your case... You have a disputable mesh: High aspect ratio, sudden changes of element sizes come into view. It looks unnecessary complex for me.

My recommendation: Start with an setup more straightforward, less nodes, lower Re, stay with blockmesh. Use first simple 2d and not a wedge. And look what happens.

> I am pretty disappointed with OpenFoam. It is a very opaque tool

I don't sign this. In contrary, nothing is hidden here. But you need to understand what you do. OpenFoam is not pocket calculator.

Thanks for the input !

As said in my first post, my progression went from a case with a provided mesh (attached, it looks complex with elongated structures to me) to the same physical equations being solved with my geometry (created with gmsh, that I knew how to use).

I understand how the entire process looks debatable from the outside, but I think my mistake instead comes from thinking it would be simple and easy to extrapolate from this case to my case. I've been at it for months now.

I'm pretty sure part of the blame lays with me, but I still can't help but feel it shouldn't be this hard to move from one mesh to another. On the whole, OpenFoam seems extremely sensible to mesh quality.

Nevertheless, you need a point where to start. If changing the mesh is not an option that change the physics. Start with potential flow, change to laminar and then to low Re flow. And look what happens.

Often, boundary conditions are unphysical or at least introducing instability.

A sketch of your case with the b.c. would be helpful. So wen don't need to analyze your files in the first step.

I'll take the advice on rolling back and starting simple, even if it feels a lot like going backwards.

In that case, I'll stick to the same solver, that will reduce tedious adaptation tasks. Adapting the case to pure 2D with a simpler mesh is very straightforward, and my continuation method has low Re in there to begin with.

Not sure how to quantify "less nodes", and I'm definitely not diving into blockMesh. I have made a simpler mesh (picture attached, as well as the whole case) and calculations run until about the same viscosity. As before, there's a pressure blowup at the top right that gradually breaks everything down (again, picture attached).

I agree about the boundary conditions - they are clearly so critical I had to include them in full. However, a drawing can't hurt. Attached is the plane I'm interested in (it's axisymmetric, I control the viscosity to get the high Re I want). Notice that the wall is an infinitely thin baffle.

Maybe there's something wrong with these but I have difficulty to see how they produce a blowup... I'll stress again I don't observe this in the NASA case.

On the other hand, this could have something to do with the "relative" small dimension in R ? I'll think about it.

I don't understand your b.c. for U. You have three inlets with differen velocities. Are there regions where the meet? Most probably. What happens in the corner elements?

Not exactly. In my sketch you can see all velocity conditions in x (the others are all 0 or zeroGradient) the inlet and coflow regions that have different but similar profiles. These are my two inlets, that are both velocity-driven in x (with u_y=u_z=0). The top boundary, far, is more of a no-slip condition ; it does not bring fluid to the domain, merely drives a weak parallel coflow to the jet axis.

The two inlet profiles meet at the wall where all three BCs (inlet, wall, & coflow) should be compatible and go to u<-(0 0 0).

As for far & coflow, I believe these two are compatible with u<-(0.1 0 0), this corresponds to the maximum of the hyperbolic tangent profile at r=r_1. Again, at the x,r=0,r_1 corner, I expect u_y=u_z=0.


The question for the top-right corner (x,r=L,r_1) is more interesting - I didn't think a problem could arise from a fixed Value going into a zeroGradient.


I'm not entirely sure what happens at the bottom with the wedge condition to be perfectly honest.

Now in the case files there are three codeStreams because of the initial condition that is taken to be the extension of the boundary condition at x=0 all across the domain. In other words at t=0, there is a cylindrical double hyperbolic tangent with u=0 at the shear layer in place of a jet (it's very clear at step 0 of the Viscouscont script).

This could be an issue I suppose, but my blowup occurs so late that I think it's a stretch. I think of this initial condition as more of a shortcut to get results faster.

> The two inlet profiles meet at the wall where all three BCs (inlet, wall, & coflow) should be compatible and go to u<-(0 0 0).

yes, I see.

I recommend calculation to time shortly before exploding. Then reduce time step by far and restart calculation. In the p and U graphics you should where the irregulariyte starts.

I'm not sure the timestep changes anything in a steadyState computation - it certainly had no influence when I tried changing it...
Yeah I will take that piece of advice concerning the origin of the unsteadiness (it could be originating from somewhere else as the top right corner). This (and refining / complexifying the mesh) has been my process so far with limited success. It's a very time-consuming process, too...

Refining the mesh ist not the first thing I would do. Rather omitting the tanh b.c. firts end replacing by a wall.

You may change to pimple Foam and use my timestep advice. I don't use simpleFoam very often because pimpleFoam finishes with steady state too. But may be you may stop the simulation at at certain point and re-run it "slower", means slower convergence. It is important to see where the problem ist.

I see problems with the mesh. But of course, I am not sure.

Please check the mesh quality with "checkMesh -allGeometry -allTopology", and see any differences in the mesh metrics for the meshes you have. Please do use the same boundary conditions you have used for the NASA case as a first step (just try to replicate the NASA setup as closely as possible.)

For "simpleFoam", using "consistent true;" (SIMPLEC; you will need to change the relaxationFactors, especially p=1) and "bounded" for the schemes can improve the numerical stability. See the web for further info.

And a small warning: OpenFOAM is not free and, to my knowledge, it was never free or was never meant to be. I think it is only another business model, which Google and Facebook use a version of it as well (they are also not "free".) So if you do serious/paid work and are short of time, you may need serious/paid help or training.

Just another thought: What kind of b.c. do you have at the upper and lower wall? Your input stream is full developed there and the b.c. should bi symmtrical or slip. If you have a wall (noslip) bent the wall slightly around the corner.

OpernFoam is free in the sence that is under the Free Software License.
https://openfoam.org/licence/

But it is not a free lunch.

Let me get this straight - you are telling me to replace the inlet BCs that are driving the flow by walls ? What would that achieve ? Make sure that the blame lies with the mesh not the BC ?

It's hard to disagree with that statement, but I guess my reflex recently to refine problematic regions and this turned out to be frustrating and slow. I manage to get a pretty good idea of where the jumbo starts by reducing writeInterval most of the time.
I actually have seen this already ; I tried with a consistent scheme for my case, though it is not required in the NASA case, and it doesn't save me. I hadn't read about the relaxation factor though.
My divSchemes are bounded, my gradSchemes limited.
Guess I learned that lesson the hard way :)
It's all on the graphic and in my previous posts :

• bottom boundary is not specified in my axisymmetric case (I thought wedge did the job automatically)
• top boundary is called far, it has a no slip condition. It has atmospheric pressure and a prescribed velocity that is purely axial with the jet and coincides with coflow maximum

> Let me get this straight - you are telling me to replace the inlet BCs that are driving the flow by walls ?

only the two tanh conditions. We have to be sure that they don't arise a problem.

> top boundary is called far, it has a no slip condition. It has atmospheric pressure and a prescribed velocity that is purely axial

I don't recommend this. If there are tiny differences in the velocity you have nonphaysical b.c. I recommend at least first a slip b.c.

It seems that in the nozzle_stupid case, which is real 2D with a simplistic mesh, the problem arises from a gradual pressure drop at the top right corner.
In my case these are the same. I prefer @HPE suggestion to go back to parabolic profiles then replacing them with walls. I feel like a non-driven cavity would just not move.
I agree that is very sensible. For the nozzle_stupid case, I went back to parabolic BCs and I tried changing the velocity condition at the outlet from zeroGradient to pressureInletOutletVelocity as in the NASA case. The latter broke my case (it now blows up in a few timesteps with negative pressure at the inlet).
I tried that with the simpler mesh I made and posted for the nozzle_stupid case. Logs files are attached, I don't see anything useful to be extracted from these. Both have cells with low determinant, but that doesn't help me much.

Summary
 

To sum things up :

• I went from axisymmetric to 2D
• I simplified massively my mesh
• Moved back to a parabolic velocity profile

Currently my calculation breaks up because of a gradual pressure blow up at the top-right (see attached image).

Ok I have a fix it works fine now.

The secret sauce was a much stronger coflow. with a 1% coflow I was getting instabilities, with a 5% one I get stable stuff, even when dropping continuation, reducing mesh quality, and shrinking the box.

Thanks for your help folks.