[mkhm]

Dears,



I am an OpenFOAM user. I use this software to simulate 2 dimensional (planar) nozzles. The geometry is obtained by Method of characteristics. As you know, the input of such method is the gamma (=cp/cv), the desired Mach number at outlet and the number of characteristic lines. My question is about the validation of my simulations.

• I use the totalPressure and totalTemperature as boundary condition at inlet.
• I use waveTransmissive boudnary condition at outlet.

The simulations converge well and I do the mesh sensitivity analysis. The higher is the number of characteristic lines and the more refined is the mesh, you get a higher mach number at outlet. The mach number obtained by simulations are very close to the desired value of Mach which was inserted to get the geometry based on the method of characteristics. However, if I use the area ratio between the outlet and the throat and I consider a 1D isentropic model (isentropic equations), I see that I should get a higher Mach number.



With the method of characteristic, the flow is uniform at the outlet. So, by the conservation of the mass, I should get the same mach number in 2D than 1D at the outlet.



So my question is why I am not getting the same results ? Why Mach number with 2D simulations is lower than the 1D (like instead of 7, I get 6) ?



I checked the area at the outlet and the throat as I thought that I might have loosed some precision passing from gmsh to OpenFOAM. But with that area ratio, I should get higher mach number.



Someone sees why it's not the case ? Is there a problem with the solver ? With the boundary condition ? with my hypothesis that I should get the same results than 1D ?



Please help me.

[LuckyTran]

Honestly I hear nothing alarming. The solver probably isn't broken and there is a whole laundry list of reasons why they could be different.


Is this a C-D nozzle or just the diverging part of the nozzle? But just to give an example: 1D theory says the critical pressure ratio for a C-D nozzle to be choked is ~1.8 but in reality, and if you run 3D CFD, it chokes much earlier (as low as 1.2). 1D theory always thinks the flow is uniform and that turns out to be quite a difference.

[mkhm]

Originally Posted by LuckyTran

Honestly I hear nothing alarming. The solver probably isn't broken and there is a whole laundry list of reasons why they could be different.


Is this a C-D nozzle or just the diverging part of the nozzle? But just to give an example: 1D theory says the critical pressure ratio for a C-D nozzle to be choked is ~1.8 but in reality, and if you run 3D CFD, it chokes much earlier (as low as 1.2). 1D theory always thinks the flow is uniform and that turns out to be quite a difference.

It is a converging-diverging nozzle. My supervisor told me that I should get nearly the same Mach number than 1D study, otherwise the simulations are not valid because the conservation of mass is not respected.



I use ParaView to get the mass flow rate at inlet and outlet, the results are similar. Ideally, I should check the mass flow rate for the intermediate cross sections (between the inlet and outlet). But for that I should write some post processing tools. As the fluid is not uniform between the inlet and outlet, by considering the integral of rho U dS (what ParaView uses) we have a different mass flow rate than inlet/outlet. That comes from the fact that this integral takes into account only the normal contribution of the velocity. So some time is needed to write another postProcessing tool.



You are right. The list of things that can go wrong is long. However, I want to make sure that my 2D results are correct. Do you know a reference or a book where we say that at the outlet of a well designed (where the method of characteristic is giving the geometry) 2D case, we should not necessarily have the 1D Mach number ?



Otherwise do you know any methods to check the validity of a 2D case ? The focus of my project is not on the method of characteristics. Maybe the matlab code that I am using to have the geometry is not a good one. But I achieve the Mach number that I introduce to this program as input ...

[LuckyTran]

Originally Posted by mkhm

You are right. The list of things that can go wrong is long. However, I want to make sure that my 2D results are correct. Do you know a reference or a book where we say that at the outlet of a well designed (where the method of characteristic is giving the geometry) 2D case, we should not necessarily have the 1D Mach number ?

Forget that it was designed using method of characteristics. You have a nozzle contour somehow? It doesn't matter where it came from. You are asking why the 1D and real flow is different? Also the method of characteristics is normally used only to design the diverging part of the nozzle. What about the converging part? How did you design it?

My example comes from venturi's (which are well designed c-d nozzles, but not designed for M>>1 at the exit): You can find a reference. Btw this reference is a sort of bible for building critical flow venturi's / sonic nozzles. Just look at figures at the end if you are too lazy to read the full text. Check out Fig. 10 in particular. In short, 1D theory is nice and describes a lot of things, but it also misses at a lot of things. But notice the year of publication! 1961! Isn't it amazing how this was common knowledge in 1961 yet unbelievable now?

Now I don't know why you get a number like Mach 6 instead of 7 and whether this is a small or big difference. But to expect that a 2D/3D simulation reproduces exactly the same result as 1D theory sounds quite... academic. Maybe in a pipe of constant cross section with no flow.

[mkhm]

Originally Posted by LuckyTran

Forget that it was designed using method of characteristics. You have a nozzle contour somehow? It doesn't matter where it came from. You are asking why the 1D and real flow is different? Also the method of characteristics is normally used only to design the diverging part of the nozzle. What about the converging part? How did you design it?

My example comes from venturi's (which are well designed c-d nozzles, but not designed for M>>1 at the exit): You can find a reference. Btw this reference is a sort of bible for building critical flow venturi's / sonic nozzles. Just look at figures at the end if you are too lazy to read the full text. Check out Fig. 10 in particular. In short, 1D theory is nice and describes a lot of things, but it also misses at a lot of things. But notice the year of publication! 1961! Isn't it amazing how this was common knowledge in 1961 yet unbelievable now?

Now I don't know why you get a number like Mach 6 instead of 7 and whether this is a small or big difference. But to expect that a 2D/3D simulation reproduces exactly the same result as 1D theory sounds quite... academic. Maybe in a pipe of constant cross section with no flow.

The converging part indeed is not coming from the method of characteristics. For the simplicity, I did a symmetrical geometry or in other terms, the converging part is the mirror of the diverging part obtained by MOC.


Thanks for your explanation. I will read this document. I just need as many documents as possible to show to my supervisor why I should not get the same results than 1D.



Assume that we should not get the same results than 1D, how do we can so validate the 2D results ? The boundary conditions at outlet somehow dictates the behavior of the fluid within the nozzle. The sensitivity analysis ,convergence of the simulations and checking the conservation of the mass at the inlet and outlet are enough ?

[mkhm]

I forgot to mention that I consider an inviscid flow and there is no boundary layer as the boundary condition for the velocity at the wall has been fixed to slip.

[mkhm]

Should the 2D results be isentropic ? because if I calculate the variation of entropy between the inlet and outlet, it is not equal to zero. And that's explain why I don't get the same mach number than 1D. But the question is am I supposed to have some isentropic behavior ? The only thing that the solver does is to save Euler equation for conservative variables. Not to fulfill the isentropic formulation.



The whole day today, I searched for kind of figures where we compare eventually 1D and 2D. But I didn't succeed.



I would be very grateful if someone can help me with this problem.

[LuckyTran]

First you claimed that you have a beautifully designed nozzle, but now you tell us it is mirrored. People that make venturi nozzles would laugh at you. Mockery aside, your laundry list is now growing. =)

So what equations are you solving? Just setting slip walls does not mean it is inviscid. You still have velocity gradients and viscous effects if you are solving N-S. Are you solving N-S or Euler? Again I refer to Fig. 10. Viscous effects account for a small effect near the wall, centrifugal effects (i.e. non-uniform flow effects) are much greater.

Is the velocity you get at the outlet purely axial? And is it uniform? If no to either, then the 2D result will deviate with 1D theory. Check out this paper on averaging non-uniform flow.

There are volumes upon volumes of books dealing with numerical validation and verification. OpenFOAM has been used for over a decade by many people with innumerable papers published. The problem is you need to convince yourself whether you trust the solver or not. And only you can convince yourself. I won't give any recommendations here. Even if you were convinced that solver is accurate, I wouldn't ever trust that there is no user error, since you haven't provided any pictures reassuring anyone that you have even a converged solution (not that I am requesting them).

Originally Posted by mkhm

The boundary conditions at outlet somehow dictates the behavior of the fluid within the nozzle.

Wait, what? If you have M>1 at the outlet then the BC at the outlet doesn't matter. Something doesn't make sense.

[mkhm]

Originally Posted by LuckyTran

First you claimed that you have a beautifully designed nozzle, but now you tell us it is mirrored. People that make venturi nozzles would laugh at you. Mockery aside, your laundry list is now growing. =)

If I am asking a question on this forum, that does mean that I take the risk that people mock at me. However, in answer I would like them to pinpoint the real problem to be more constructive. I said earlier :

I checked the area at the outlet and the throat as I thought that I might have loosed some precision passing from gmsh to OpenFOAM. But with that area ratio, I should get higher mach number.

So already, I have a doubt about the diverging part of the nozzle. I do not understand why the same area ratio should give higher mach number at outlet when the 1D equation is used. However, I get with my 2D simulation, the same mach number that I introduce to the program. Isn't it strange ? after we can discuss that getting a mach number of 6 instead of 7, is not a big deal and as the sonic line is not a real line in 2D case and it is a curve, we can accept this difference. I don't know honestly. I asked my question on the main CFD forum hopping that people encountered the same issue and having an answer for that.

for the converging part, is there any rule/method to be able to claim that is well designed ? Should not flow adapt itself to have the sonic condition at the throat ?

Please see the shape of the nozzle below.

Originally Posted by LuckyTran

So what equations are you solving? Just setting slip walls does not mean it is inviscid. You still have velocity gradients and viscous effects if you are solving N-S. Are you solving N-S or Euler? Again I refer to Fig. 10. Viscous effects account for a small effect near the wall, centrifugal effects (i.e. non-uniform flow effects) are much greater.

I am solving Euler equation. There is no viscosity terms. I knew that there is this non-uniformity between the throat and the outlet. What sincerely I don't know is how big their influence could be.

Is the velocity you get at the outlet purely axial? And is it uniform? If no to either, then the 2D result will deviate with 1D theory. Check out this paper on averaging non-uniform flow.

We can say that the flow is uniform at outlet and the velocity has only one component. Why could I say that ? Please have a look at the figure below for the Mach variation within the geometry. If I do the calculation for the mass flow rate with "Surface Flow" option of ParaView, it considers integral of rho U dS. If the velocity has two components, this integral would take only the axial velocity. So, if the mass flow rate that you obtain is different from the inlet, you can say the flow is not uniform. But, at the outlet I get the same flow rate than inlet. So, I assume that non uniformity is very very low. Otherwise, I would not get the same inlet mass flow rate.

There are volumes upon volumes of books dealing with numerical validation and verification. OpenFOAM has been used for over a decade by many people with innumerable papers published. The problem is you need to convince yourself whether you trust the solver or not. And only you can convince yourself. I won't give any recommendations here. Even if you were convinced that solver is accurate, I wouldn't ever trust that there is no user error, since you haven't provided any pictures reassuring anyone that you have even a converged solution (not that I am requesting them).

I have more doubts about the errors that could happen by myself by setting something wrong than having doubts about the openfoam solver. That's why I tried to explain everything in details.

Wait, what? If you have M>1 at the outlet then the BC at the outlet doesn't matter. Something doesn't make sense.

In my very first message, I said that for the output I use the "waveTransmissive" boundary condition for the pressure. Without that, I would get normal shocks. Euler equations are not linear. So having a normal shock could be one of the solution of the Euler equation. I don't fixe the pressure at outlet to the 1D supersonic pressure solution, because I would see the solver crashes. So, I try to fixe the pressure in a far field (openfoam will extend the geometry virtually and at the the end of this far field extension, you have the fixed pressure).

To sum up, there are two fundamental questions : Are we supposed to get a 1D Mach number at the outlet when the flow is uniform ? If not, how far it could be our result ?

[LuckyTran]

Originally Posted by mkhm

for the converging part, is there any rule/method to be able to claim that is well designed ? Should not flow adapt itself to have the sonic condition at the throat ?

I already gave you someone's thesis for how to design the throat. Btw those recommendations are now part of manufacturing specs for sonic nozzles. Hence why it is a bible. The converging part is much easier to design than the diverging part. In the diverging part you use method of characteristics. In the converging part, simple radius filleting with some rules of thumb are given. The challenge in the converging part is minimizing boundary layer growth (which is not an issue for your case) and trying your best to deliver a uniform flow to the throat. Unfortunately, you mirrored the diverging part and don't have a nice smooth throat. Design rules for rocket nozzles are similar, but the throat is extended (more straight section is added) for combustion to take place.

The sonic condition does not happen exactly at the throat but (as you already correctly understand) does not occur at a single axial location but as a 3D disc. The axial extent of this disc is on the order of the throat curvature (i.e. the throat diameter, because that's what determines the level of centrifugal forces).

Originally Posted by mkhm

If I do the calculation for the mass flow rate with "Surface Flow" option of ParaView, it considers integral of rho U dS. If the velocity has two components, this integral would take only the axial velocity. So, if the mass flow rate that you obtain is different from the inlet, you can say the flow is not uniform. But, at the outlet I get the same flow rate than inlet. So, I assume that non uniformity is very very low. Otherwise, I would not get the same inlet mass flow rate.

This surface integral always correctly returns the mass flow rate so you won't be able to detect non-uniformity this way. So you should still check. But looking at your Mach contours, they look pretty nice. Can't tell what is the +/- Mach number variation at the outlet alone but it is a starting point.

Originally Posted by mkhm

In my very first message, I said that for the output I use the "waveTransmissive" boundary condition for the pressure. Without that, I would get normal shocks. Euler equations are not linear. So having a normal shock could be one of the solution of the Euler equation. I don't fixe the pressure at outlet to the 1D supersonic pressure solution, because I would see the solver crashes. So, I try to fixe the pressure in a far field (openfoam will extend the geometry virtually and at the the end of this far field extension, you have the fixed pressure).

If your outlet is supersonic (and the Mach number contour you provide convincingly says that it is since M>>1) then you don't need the wave transmissive BC. So I would rule that out. I find it odd that you said the outlet BC affect anything, because it shouldn't. Not really important, but the wavetransmissive BC doesn't virtually extend the geometry. If you care to learn more about it, read about it, if not, don't.

Originally Posted by mkhm

To sum up, there are two fundamental questions : Are we supposed to get a 1D Mach number at the outlet when the flow is uniform ? If not, how far it could be our result ?

In general, no.

You can see in your Mach number contour at the upstream part of the diverging section that the Mach number is not uniform (I can see green, yellow, and red at the same axial location). If you want your OF to match the 1D theory, it should always be uniform in the diverging part.


Btw, can't you just run a case where you simulate only the diverging part of the nozzle without any contracting part?

[mkhm]

Thanks LuckyTran for your answer.

Originally Posted by LuckyTran

In the converging part, simple radius filleting with some rules of thumb are given.

I tried that for the converging part but after I removed it as I had difficulties to find the real surface area at the throat. And after as I was not sure if modifying the real curvature (by designing like that the curvature of the converging part, you are somehow obliged to change the coordinate of first two points of the diverging part) was the reason of not getting the 1D mach number at outlet, I started keeping the exact same curvature given by the method of characteristics.

Originally Posted by LuckyTran

Btw, can't you just run a case where you simulate only the diverging part of the nozzle without any contracting part?

Good idea. I will try that.

However, I would like to check for other codes giving the supersonic curvature by MOC. Do you know if MOC is still reliable for hypersonic cases ? My experience show that with much lower mach number, between 1-4, the outlet mach number is much closer to 1D one. However, when we go for higher mach number, the difference is big. I am wondering if this method does not lose its accuracy for higher mach number cases. Do you know any reliable contouring codes ?

Thanks a lot for your time