[fportela]

Hello!

I have been trying to get results on a turbulent supersonic channel flow test case with rhoCentralFoam. have been trying to get results on a turbulent supersonic channel flow test case with rhoCentralFoam.

The test case is the same as Test A here: http://journals.cambridge.org/action/displayAbstract?fromPage=online&aid=340246

A summary of the (my) test is as follows:

• Domain/mesh
  • Dimensions: 12 x 2 x 4
  • Mesh: 64 x 64 x 64 (and 128 x 128 x 128)
• Boundary conditions
  • Periodic in stream- and spanwise directions
  • At the walls: no-slip, isothermal (T=1), zero gradient on pressure
• Initial conditions
  • Uniform pressure, temperature and density
  • Parabolic velocity profile with random perturbations

The flow is sustained by a source term which is added to the momentum and energy equations.

After running for a while, the flow seems to become laminar (my fluctuations are much lower than those obtained with the DNS) and displays some compression waves travelling through the domain (I attached some contour plots to illustrate this).

Right now I am not using any explicit SGS model, though I have tried Smagorinsky and did not observe much difference. I also tried the same case on a finer mesh and this did not eliminate the wavy pattern.

Has any one experienced something similar? Any suggestions on how to handle this?

Cheers,
Felipe

[cfd.with.openfoam]

Hi Felipe,

I am super interested in your case. So, some questions -

1. Were you able to get rid of the problem?
2. For the incompressible turbulent channel we have the mom. source term in the fvOptions file which we can use to keep the flow going in case we use cyclic BCs (if we use directMapped then we have the setAvg. option). So, you used cyclic BCs - but how did you add source terms to both momentum and energy equations.? More specifically - how to calculate each of them and then how to implement them? Can you please shed some light? Maybe its a stupid question and I should read a little bit more, but maybe its not?

By the way I have been trying to simulate a fully developed compressible channel flow in the RANS framework with rhoCentralFoam using both cyclic and directMapped BCs, without much luck so far. When I use the former (i.e cyclic BC) the average velocity quickly moves towards 0 (which points towards the source term) and in the case of latter (directMapped) my flow gets unstable and starts shedding vortices.

It will be great if you can share your case. If that sounds OK then my id is cfd.with.openfoam@gmail.com

Thanks very much

[fportela]

Hi there,

I did get rid of the problem by stretching the mesh, near wall refinement eliminated the strange compression waves in my original post, however in order to achieve good statistics (I was doing ILES at the time) I had to use a rather fine mesh (256^3) otherwise the flow was barely turbulent.

I am unaware of the momentum source option in the fvOptions dict. In my case I had the source term introduced directly in the solver: I take the momentum source to be given by a constant acceleration (to be determined) such that in the momentum equations we introduce the term , conversely, for the energy equation we must introduce . In my implementation I had a routine which would calculate based on the bulk velocity, however a good estimate for is given by with the density at the wall and the bulk density (h is the half channel height). For the incompressible case I assume the density ratio reduces to 1 so you are left with .

Let me know if something is not clear (if we want me to guide you through the steps above). Regarding implementation, what solver are you planning to use? rhoCentralFoam?

Regarding your last comment, for the cyclic case, if makes sense that your velocity goes to zero, since there is no force driving the flow. I have never used directMapped so I can't help you there, I'm not even sure of what it is supposed to do (maybe you can point me to some documentation?)

Cheers,
Felipe

Originally Posted by cfd.with.openfoam

Hi Felipe,

I am super interested in your case. So, some questions -

1. Were you able to get rid of the problem?
2. For the incompressible turbulent channel we have the mom. source term in the fvOptions file which we can use to keep the flow going in case we use cyclic BCs (if we use directMapped then we have the setAvg. option). So, you used cyclic BCs - but how did you add source terms to both momentum and energy equations.? More specifically - how to calculate each of them and then how to implement them? Can you please shed some light? Maybe its a stupid question and I should read a little bit more, but maybe its not?

By the way I have been trying to simulate a fully developed compressible channel flow in the RANS framework with rhoCentralFoam using both cyclic and directMapped BCs, without much luck so far. When I use the former (i.e cyclic BC) the average velocity quickly moves towards 0 (which points towards the source term) and in the case of latter (directMapped) my flow gets unstable and starts shedding vortices.

It will be great if you can share your case. If that sounds OK then my id is cfd.with.openfoam@gmail.com

Thanks very much

[cfd.with.openfoam]

Dear Felipe,

Thank you for your message. Sorry for the delay but I had to suddenly look at another problem. I am still super interested in your case. Since you are happy with your results I will try to follow your method (instead of fvOptions or directMapped) for handling this problem.

Let me first talk about the momentum equation in rhoCentralFoam -

// --- Solve momentum

solve(fvm::ddt(rhoU) + fvc::div(phiUp));

This is a strange form of momentum eqn for me - as U and p are combined for some reason. But if we think about the standard form of mom. eqn (e.g eqn 2 here) - we just need to add -nabla(p) = source term = rho*(rhoW/rho0)*(Utau^2)/h to the above eqn resulting in the following -

solve(fvm::ddt(rhoU) + fvc::div(phiUp) + rho*(rhoW/rho0)*pow(Utau,2)/h);

Does this look OK? I am still not sure about the term - rhoW/rho0? Can you please explain its origin?You also wrote -"In my implementation I had a routine which would calculate based on the bulk velocity"- In this approach I don't see bulk velocity playing any role.

Thank you for your time.

[fportela]

Hi again,

I never said I was happy with my results, as I said: my application was implicit LES on a supersonic channel flow and I found rhoCentralFoam to be a bit too dissipative for such an application in comparison to other solvers (I used the same configuration as Coleman et al). I also used rhoPimpleFoam (I am aware it is not suitable for supersonic applications but due to the absence of shocks I gave it a go) and found it to give better statistics for the velocity components (both mean and RMS) however the temperature field came out a bit messed up... However from the solver readily available, indeed rhoCentralFoam seems to be the best to go with (for now), you can also try sonicFoam but I recall I had some issues with the energy equation.

Regarding the momentum equation in rhoCentralFoam, it is not solving the momentum equation per se as you would run into stability issues using a central scheme: have a look at the schemes here where the solver is based upon.

As for the source term, that seems about right. Sorry for creating confusion with the controller, the issues is the following: if you simply plug in that expression into you equation, your solution (unless you have a massive resolution) will not display the expected Utau, this is because your discretisation itself is introducing dissipation; on the other hand, if you add another term which accounts for the difference between the solution and the desired parameter (bulk velocity in my case) then you will indeed have the desired Utau (or bulk U).

Conclusion: depending on what you want, it may be necessary to increase your source term to account for artificial dissipation (compressible solvers usually have a lot of this), you can add a control variable or just tweak the parameter yourself, on the other hand, if you intend to compare your results with that of others, I recommend you find out what value they use for the source term (through the expression you used) and use that in your case, otherwise you're not really solving the same problem.

I hope things are a bit more clear but of course let me know if not!

Cheers,
Felipe

PS: I can send you the test case if you're interested, though it's not that different from the channel case that comes with OpenFOAM...

Originally Posted by cfd.with.openfoam

Dear Felipe,

Thank you for your message. Sorry for the delay but I had to suddenly look at another problem. I am still super interested in your case. Since you are happy with your results I will try to follow your method (instead of fvOptions or directMapped) for handling this problem.

Let me first talk about the momentum equation in rhoCentralFoam -

// --- Solve momentum

solve(fvm::ddt(rhoU) + fvc::div(phiUp));

This is a strange form of momentum eqn for me - as U and p are combined for some reason. But if we think about the standard form of mom. eqn (e.g eqn 2 here) - we just need to add -nabla(p) = source term = rho*(rhoW/rho0)*(Utau^2)/h to the above eqn resulting in the following -

solve(fvm::ddt(rhoU) + fvc::div(phiUp) + rho*(rhoW/rho0)*pow(Utau,2)/h);

Does this look OK? I am still not sure about the term - rhoW/rho0? Can you please explain its origin?You also wrote -"In my implementation I had a routine which would calculate based on the bulk velocity"- In this approach I don't see bulk velocity playing any role.

Thank you for your time.

[cfd.with.openfoam]

Hi Felipe,

I think I should have used "OK" instead of "Happy" to describe your experience with rhoCentralFoam. Sorry for that.

After reading your message I am starting to lean towards the directMapped method again as it is very clean - simply because we can maintain the Ubulk automatically with the setAverage option. So, along that line -

In my first message, I mentioned that with directMapped my channel starts shedding vortices. I managed to get rid of the issue by increasing the domain length but with the longer channel -

I start off really well (i.e the velocity looks quite OK - picture is attached) but as the simulation progresses my BL almost disappears (see attached picture). Can you guess the reason for such a behaviour? I think it might have something to do with compressibility - so my next step would be to try a 10 times lesser velocity with pimpleFoam.

Any comments? Thank you

[cfd.with.openfoam]

Sorry - forgot to clarify that I am only doing RANS for now and not LES.

[fportela]

What RANS model are you using? I have little to none experience with RANS but as I recall some are more appropriate for wall bounded flows than others, this associated with the over diffusive behaviour I mentioned above may be influencing your results, I don't see how compressibility could affect the velocity field (if you read the paper of Coleman I posted above this is one of their comments).

What boundary conditions do you use for temperature? The velocity profile seems ok (not laminar) but indeed k seems quite low...

Can you also give some details on the mesh?

Cheers,
Felipe

[cfd.with.openfoam]

Hi Felipe,

Answers to some of your questions -

1. k-Omega SST
2. Total Temperature BC at the inlet (value for To is 340 K)
3. Pressure is zeroGradient since it is subsonic
4. Velocity as I said is directMapped with an average value of (213 0 0). As you can see the avg. value is maintained throughout but the profile is messed up.

Regarding your comments on diffusion - to me more diffusivity would have meant a thicker BL and not the other way around. This seems to be a less diffusive behaviour Or maybe I have missed your point completely.

I have to admit that my mesh might not be that well designed (by the way a picture is attached - it is 50*40 in the streamwise and wall normal directions). I think y+ should be less than 5, but not sure . So, I will certainly look into that - maybe I should check this out before going for the incompressible case with pimple/piso-Foam

Any comments?

[fportela]

But what is the BC at the wall for temperature? If you use constant temperature I believe the core of the flow should be hotter and not the other way around.

Again, I am not sure about the directMapped thing: by including the source terms in the equations as I had done you could also account for the energy introduced by the source term (by adding it appropriately to the energy equation) however since you are rescaling the velocity profile (I assume this is what directMapped does) shouldn't you be re-scaling temperature as well?

Regarding diffusivity, indeed I agree with you, the point is: it should be thicker with respect to what? I assume your initial condition is not that of a fully developed turbulent flow (maybe I'm wrong ) so comparing just those two profiles doesn't say much. I guess what I mean is that you still have a turbulent profile (with a thicker BL than it should but that we have no means to check) and k is reduced due to this dissipation.

W.r.t. the mesh, I guess if you are going down to y+<5 you should be fine without any special wall treatment. So I would say the issue is not mesh-related.

Have you compared the solution with and without the turbulence model on? what happens to UMean?

Cheers,
Felipe

Originally Posted by cfd.with.openfoam

Hi Felipe,

Answers to some of your questions -

1. k-Omega SST
2. Total Temperature BC at the inlet (value for To is 340 K)
3. Pressure is zeroGradient since it is subsonic
4. Velocity as I said is directMapped with an average value of (213 0 0). As you can see the avg. value is maintained throughout but the profile is messed up.

Regarding your comments on diffusion - to me more diffusivity would have meant a thicker BL and not the other way around. This seems to be a less diffusive behaviour Or maybe I have missed your point completely.

I have to admit that my mesh might not be that well designed (by the way a picture is attached - it is 50*40 in the streamwise and wall normal directions). I think y+ should be less than 5, but not sure . So, I will certainly look into that - maybe I should check this out before going for the incompressible case with pimple/piso-Foam

Any comments?

[cfd.with.openfoam]

Hi Felipe,

Regarding Temperature - Upper and Lower walls have T=340K (at the walls static and stagnation are the same as there is no velocity). The inlet has To=340 K. As we know To has two components T (static) and the velocity component U^2/2Cp. Since in the middle of the channel we have a velocity of 213 m/s which makes the second component = around 25 K - which basically means that T (static) = To-25 K = 340 - 25 = 315 K. So, I don't think we have any problem there.

So, if I use To I don't think that I need to use any rescaling on temperature as the TotalTemperature BC will pick the velocity up to give me a corresponding T (static). OR maybe I am wrong and I should use directmapped on T as well.

Yeah I am also worried about the source term in the energy equation which is why I would like to test an incompressible case.

Diffusivity - I see your point but in a fully developed channel we should expect the two BLs to meet at the center but here the BLs disappear far away from the channel center. So one might say that it is not diffusivity problem.

Without Turbulence Model - Well I haven't thought about it but maybe I will give it a shot as well.

Thank you for your comments and time.

[fportela]

Well if you have the walls set to constant temperature then there is something wrong with the temperature profile, so this suggests the total temperature approach may not be the best option (though I am unable to give a better alternative than the inclusion of the source term).

Here is what I have for rhoCentralFoam with the source term included (I didn't send before as I didn't have the files with me):

solve(fvm::ddt(rhoU) + fvc::div(phiUp);

becomes

solve(fvm::ddt(rhoU) + fvc::div(phiUp) - source);

and

 solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU)
        );

becomes

solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU) - (U & source)
        );

Indeed I agree with your point on the BCs, the U profile does seem too flat at the centreline which leads me to think of two possibilities: maybe the simulation didn't run for long enough (I'm not sure how to estimate this in RANS but for my LES I just waited until the running mean reached a "steady state") or the scaling in the plot is not adequate, I don't see how the flow could develop to a state where the BL on both sides don't meet.

The reason why I keep insisting on a dissipation related issue is that for my implicit LES I simply simulated the flow on a "coarse" mesh: I got reasonable (but still too diffusive) results for meshes finer than 128^3.

Let me know how the simulation without model goes, compare the results with that of you k-omega and maybe it will give hints on how to proceed.

Cheers,
Felipe

Originally Posted by cfd.with.openfoam

Hi Felipe,

Regarding Temperature - Upper and Lower walls have T=340K (at the walls static and stagnation are the same as there is no velocity). The inlet has To=340 K. As we know To has two components T (static) and the velocity component U^2/2Cp. Since in the middle of the channel we have a velocity of 213 m/s which makes the second component = around 25 K - which basically means that T (static) = To-25 K = 340 - 25 = 315 K. So, I don't think we have any problem there.

So, if I use To I don't think that I need to use any rescaling on temperature as the TotalTemperature BC will pick the velocity up to give me a corresponding T (static). OR maybe I am wrong and I should use directmapped on T as well.

Yeah I am also worried about the source term in the energy equation which is why I would like to test an incompressible case.

Diffusivity - I see your point but in a fully developed channel we should expect the two BLs to meet at the center but here the BLs disappear far away from the channel center. So one might say that it is not diffusivity problem.

Without Turbulence Model - Well I haven't thought about it but maybe I will give it a shot as well.

Thank you for your comments and time.

[cfd.with.openfoam]

Hi again Felipe,

Some updates -

1. directMapped works very well for an incompressible case. Works very well means that the case reaches steady state and the velocity profile looks quite OK and it doesn't change in the streamwise direction.
2. The thin BL issue - I made a very stupid mistake as I used zeroGradient both at the inlet and the outlet. So, the pressure was not anchored and it just kept on increasing throughout the domain leading to very thin BLs
3. directMapped doesn't seem to work for a compressible case, simply because if we fix Mass Flow Rate (which is what directMapped does) then we get a pressure gradient in the results but since the flow is compressible it means T, rho and U change in the streamwise direction as well. Please see attached velocity contour. Any comments on this theory?

Due to the 3rd point above I am trying to get a case working with the sources being added directly to rhoCentralFoam - but I have a problem. First let me describe the way I have added source terms (based on your suggestions) -

This is my momentum equation with the Source term in blue -

solve(fvm::ddt(rhoU) + fvc::div(phiUp) - Source*sourceMag);

Right now I am just adding a no. (i.e no Utau calc.), with correct dimensions and size, extracted from the directMapped simulation as the source (see code below) -

volVectorField Source=rhoU;
vector forSource(1,0,0);
forAll(mesh.C(),i)
{
   Source[i]=forSource;
} 
dimensionedScalar sourceMag
(
    "sourceMag",
    dimensionSet(0, 0, -1, 0, 0, 0 ,0),
    15000
);

So, the no. is 15000 and is not unrealistic. Below is my energy eqn with Source term in blue -

        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
          - fvc::div(sigmaDotU) - (U & Source*sourceMag)
        );

All this is really simple and straightforward (after your suggestions) but it is just not working. What happens is that the Temp. and Pressure in my domain just keep on increasing. Velocity seems to stablize but it is very high. Any comments?

If possible can you please send me a dummy case where you get a compressible turbulent channel case working in the RANS framework where all the properties are not changing in the streamwise direction. Maybe I am asking too much.

Thank you very much

[cfd.with.openfoam]

Hi Felipe,

I am also attaching the contours of P,T and rho from the simulation mentioned in the third point (i.e Compressible Turbulent channel with directMapped where all the things change in the streamwise direction as well) in my last message.

Also I must say that I had to use the waveTransmissive BC to get this case to stabilize. I know very little about this BC as of now but I think I might have used it correctly.

Hope it will give you a better idea.
Thanks

[cfd.with.openfoam]

Hi Felipe,

Any comments? I am really looking forward to your response.

Thanks

[fportela]

Hi there,

Sorry for the late reply, haven't had a chance to look at this for a while.

At a first glance I agree with your implementation of the source term. Would you share a few plots of your results for the source term case?

I attached a test case I had for rhoCentralFoam with source term. It is the Mach=1.5 case of Coleman. All variables are non-dimensionalised much like he does in his paper.

With this mesh I managed to get something "close" to the DNS, without ANY turbulence model. Unfortunately I no longer have the plots nor solver available so I can't be more precise with how close it was...

The value I used for acceleration (as mentioned a few posts above) was 0.004005 (don't forget the case is non-dimensionalised).

To initialise the case (don't forget to run blockMesh) use funkySetFields (if you don't have swak4foam yet, it's worth looking into it). The initialisation is also the same as that of Coleman: parabolic profile with random perturbations on the velocity field, I worked out the remaining variables from his data.

Regarding wave transmissive BC, from what I recall it is a 'modified' (as in extended) version of the advective BC. The point is that you have waves travelling in and out of the domain and you'd like to ensure they leave unreflected. However, as far as I recall, OpenFOAM's implementation merely 'simulates' this behaviour, I gave a shot at non-reflective BCs and it was a nightmare... If you're curious on this have a look at the paper of Poinsot and Lele as I believe it's the seminal paper on non-reflective conditions.

Keep me updated.

Cheers,
Felipe

Originally Posted by cfd.with.openfoam

Hi Felipe,

Any comments? I am really looking forward to your response.

Thanks

[cfd.with.openfoam]

Hi Felipe,

Thanks for your detailed response. I am sorry for the delay - but a no. of things kept me away.

Plots from the source term case -

1. Velocity - Looks OK and there is no change in the streamwise direction contrary to the directMapped plots in my last message.
2. Pressure - Doesn't look good. Did you get a similar plot. This kind of contour generally indicates that there is something wrong with the conservation eqns. Another problem here is that I started with 6000 Pa and it continues to increase, which again points towards the source term. With a proper source term the pressure should not increase. Any comments?
3. Temperature - also starts increasing like pressure. I started with 340 K. Also note that the upper and lower walls have zeroGradient.

Your Case - I would first try your case in 2D RANS framework. However I still feel that my case was OK but my source terms are messed up. So, I would expect something similar (i.e pressure and temp. continuously increasing) in your case as well, with my source term implementation. Any comments?

Poinsot and Lele's paper - I will have to take a look at it sooner than later.

Thank you very much for your time.

[fportela]

Hey,

OK, my first question is regarding temperature: did you add a source term to the energy equation? If so, you can easily see why temperature increases: adiabatic walls mean you don't "lose" energy through the wall, but if you're adding energy through the source term your temperature will keep on rising. For this you have (as far as I can imagine) three options:
1 - Set constant temperature on the wall such that the flow develops to some equilibrium state
2 - Set the temperature gradient at the wall to a fixed value (non-zero) which should balance the energy inputed by the source term
3 - Eliminate the source term from the energy equation: in the paper of Coleman they present results for something similar, I believe it's called "case AX"

Regarding pressure, I don't see why it should remain constant: this would be the case if you would initialise your flow with the "exact" solution and carry on from there, however the flow is initialised from some non-equilibrium state so we can't say much about what happens in between. From my results, I observed that the pressure rose from the initial prescribed value, I also observed it rose more for the coarser grids: the profile of pMean is also not exactly constant across the channel, I have attached results I obtained on the case I sent you with rhoPimpleFoam (for these results I set the source term to the constant I mentioned before plus a term proportional to the difference between the target bulk velocity and the actual bulk velocity).

A note on the attached plots: as I mentioned before, my experience with rhoPimpleFoam (for this case) was that it overestimated the temperature profile, with rhoCentralFoam and the other solvers I tried the gap w.r.t. the DNS was much smaller. The plots are for a statistically stationary solution of the case I sent you before on two different meshes (the one I sent was 256^3).

Keep me updated on the results you get on the case I sent, I would expect them to be reasonably close to the DNS

Cheers,
Felipe

Originally Posted by cfd.with.openfoam

Hi Felipe,

Thanks for your detailed response. I am sorry for the delay - but a no. of things kept me away.

Plots from the source term case -

1. Velocity - Looks OK and there is no change in the streamwise direction contrary to the directMapped plots in my last message.
2. Pressure - Doesn't look good. Did you get a similar plot. This kind of contour generally indicates that there is something wrong with the conservation eqns. Another problem here is that I started with 6000 Pa and it continues to increase, which again points towards the source term. With a proper source term the pressure should not increase. Any comments?
3. Temperature - also starts increasing like pressure. I started with 340 K. Also note that the upper and lower walls have zeroGradient.

Your Case - I would first try your case in 2D RANS framework. However I still feel that my case was OK but my source terms are messed up. So, I would expect something similar (i.e pressure and temp. continuously increasing) in your case as well, with my source term implementation. Any comments?

Poinsot and Lele's paper - I will have to take a look at it sooner than later.

Thank you very much for your time.

[anshul bansal]

hello all,
i doing my project in rhoCentralfoam..

i have some question...

i) in surfaceScalarField phiv_neg(U_neg & mesh.Sf());
what mean by this '&' operator??
and mesh.Sf() define only direction or it is also have magnitude??

ii) in momentum equation
fvm::ddt(rho, U) - fvc::ddt(rho, U)

why we use fvm and fvc for same quantity ??

thanks
Anshul

[fportela]

Hi anshul,

Regarding your first question:

The '&' operator is used for inner product as you can see here. If I am not mistaken mesh.Sf() contains face normal vectors whose magnitude is the area of the face.

As for the second question, I suggest you go through this thread.

Cheers,
Felipe

Originally Posted by anshul bansal

hello all,
i doing my project in rhoCentralfoam..

i have some question...

i) in surfaceScalarField phiv_neg(U_neg & mesh.Sf());
what mean by this '&' operator??
and mesh.Sf() define only direction or it is also have magnitude??

ii) in momentum equation
fvm::ddt(rho, U) - fvc::ddt(rho, U)

why we use fvm and fvc for same quantity ??

thanks
Anshul