Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Axisymmetric source viscous terms #1095

Closed
wants to merge 23 commits into from
Closed

Conversation

FlorianDm
Copy link
Member

@FlorianDm FlorianDm commented Oct 20, 2020

I added the viscous terms to the axisymmetric source but I am not sure it is all correct. I get very reasonable, well converged results with my cases but dissipation seems a little excessive, however I cannot say for sure because there are many other factors in my own cases. There is no test case I believe. Should I try to create one with air like this one https://turbmodels.larc.nasa.gov/axibump_val.html ?

How come in the incompressible version the turbulent heat transfer is ignored? I added it here with the Reynolds analogy as eddyvisc * cp/PrTurb * dT/dy and it seems significant in my cases (with PrTurb = 0.9).

Also, I created a generalised version of the axisymmetric source in a separate class which goes with the generalised Roe and AvgGradFlow solvers with any eos. I use it for my cases with a table based fluid model (not included here) but I am not sure it is needed for the existing fluid models because they all use a constant gamma anyway. Perhaps it is because "Pressure_i = (Gamma-1.0)U_i[0](U_i[nDim+1]/U_i[0]-0.5*sq_vel)" may not be true for the "polytropic Van der Waals" and "polytropic Peng Robinson" models. If that is the case the Driver should instantiate the general one for anything other than incompressible, ideal gas and standard air.

@WallyMaier
Copy link
Contributor

Hi Florian, I was also working on the same thing (feature_axi). I was unable to find a suitable case to test, and it seems we went about it in different ways. Im happy you got something working! Could you pass along the literature you used to create these terms?

@bigfooted
Copy link
Contributor

Hi Florian, great work, and important to have. The transonic bump case seems to have enough data available for proper validation.
What are you unsure about? The actual derived source terms? If you have a document with your derivation of the source terms that would help in checking if it is correct.

Copy link
Member

@pcarruscag pcarruscag left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Florian, looks tidy, just a few comments.
Do you have a comparison between 3D and axisymmetric?
If you have time, add a regression test, it would also serve as an example for when we get questions on CFD online.

Comment on lines 1868 to 1869
if (config->GetKind_FluidModel() == 0)
numerics[iMGlevel][FLOW_SOL][source_first_term] = new CSourceGeneralAxisymmetric_Flow(nDim, nVar_Flow, config);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

elseif, otherwise it may leak memory.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does the general (non ideal?) implementation only work for standard air? (I did not think our standard air was non ideal)
Don't compare to 0, use the appropriate enum values.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh yes the 0 was a placeholder but a bad idea. I made it default to the general version now for anything but ideal gas and standard air. I think that is how it should be since the non-ideal models use the generalised solvers for fluxes too

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is completely general and works with any eos such as my table interpolation. Which btw works very well but before I can propose to integrate that I want to create a user interface in the config file to provide the data in ascii files. Standard air is of course ideal gas.

SU2_CFD/src/solvers/CEulerSolver.cpp Outdated Show resolved Hide resolved
Comment on lines 228 to 236
residual[0] -= 0.0;
residual[1] -= yinv*Volume*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][1]
- TWO3*PrimVar_Grad_i[2][0]);
residual[2] -= yinv*Volume*total_viscosity_i*(2*(PrimVar_Grad_i[2][1]-Velocity2_i*yinv)
- TWO3*PrimVar_Grad_i[2][1]-FOUR3*Velocity2_i*yinv);
residual[3] -= yinv*Volume*(total_viscosity_i*(Velocity1_i*(PrimVar_Grad_i[1][1]
+ ONE3*PrimVar_Grad_i[2][0])
- FOUR3*Velocity2_i*PrimVar_Grad_i[1][0])
+ total_conductivity_i*PrimVar_Grad_i[0][1]);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this part looks the same as for the ideal gas source, please avoid this duplication.
For example, implement this section as a method of CSourceAxisymmetric_Flow, then make the general class derived from it so it can re use the code.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay I will see how i can do that

@FlorianDm
Copy link
Member Author

FlorianDm commented Oct 21, 2020

Hi everyone, glad this subject is of interest.

I will try to make that test case.

I edit my earlier comment to avoid a mess. So I changed the terms. Still not sure but the derivation is simple.

source term viscous = (0, tau_xy, tau_yy - tau_thetatheta, u* tau_yx + v* tau_yy - q)/y, right?

then from Bird:

IMG_20201027_115003

bulk viscosity = 0, any derivative wrt theta = 0

For the generalised inviscid part I am pretty sure it is all correct including the jacobian. You can compare with very similar terms in any generalised flux jacobian like in Glaister's paper https://www.sciencedirect.com/science/article/pii/002199918890174X

@FlorianDm
Copy link
Member Author

I will do my case in 3D to compare

Comment on lines +134 to +136
void CSourceAxisymmetric_Flow::ResidualDiffusion(){

su2double laminar_viscosity_i = V_i[nDim+5];
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nicely done 👍

@FlorianDm
Copy link
Member Author

In this current form there is too much dissipation in most situations in my RANS simulations with k-omega-SST (in some situations it appears there is not enough, probably related to a low magnitude of du/dy). I don't have time to try any pure Navier-Stokes cases.

The source terms (convection and diffusion) for k and omega should also be implemented I suppose which should be very much analogous so I think I could try to do that.

Otherwise I don't see what could be wrong with the flow source terms. Please let me know if you do.

I have a couple of questions and I would really appreciate if anybody could give any input.

  1. What happens when a user specifies several sources? Will only the last one be retained since they are always assigned to the same numerics[iMGlevel][FLOW_SOL][source_first_term]? This matters most for the turbulent source since there is always already another source.

  2. Why is the axisymmetric source term for the incompressible variable density version what it is (the publication related to this does not give any derivation)? If it is correct, something about my derivation must be wrong. My derivation :

    • take the equations in cylindrical coordinates (convection term as well as the tau definition, see "Bird - Transport Processes" Annex B if needed)
    • set any d/dtheta and velocity_theta (ignoring swirl) to zero
    • rearrange into two parts, one which resembles the cartesian 2D form (r=y, z=x) and the remaining part is the source term to be added.
  • Is the above derivation wrong?

@pcarruscag
Copy link
Member

pcarruscag commented Nov 3, 2020

Hi Florian,
You can have a look in CDriver.cpp around line 1850, the logic is mutually exclusive, so if multiple sources are specified the order does not matter, the way things are coded determines which source is used. (I hope somewhere in CConfig we check for conflicting sources).
Note that the first term / second term is per solver, so there is actually a lot of room for different sources.

Regarding 2 there are more knowledgeable people about those topics out there than me. Maybe you can join the developers meeting tomorrow at 3pm CET (https://meet.jit.si/SU2_DevMeeting) and you can ask those questions to e.g. @economon or @TobiKattmann

@FlorianDm
Copy link
Member Author

Okay thank you Pedro (I meant the order in the code not in the config file). So I would just assign any k-omega axi source to numerics[iMGlevel][TURB_SOL][source_second_term].

See you tomorrow for the meeting

@WallyMaier
Copy link
Contributor

@FlorianDm When I was working on the feature_axi to attempt what you have done, I unsed this reference, which should match the incompressible version exactly.
AxiSymm_1.pdf
AxiSymm_2.pdf

@jtneedels Has ran/is running a laminar viscous wedge with your formulation of axisymmetric and 3D to do some verification. We have also tried to run the NASA TMR axisymmetric Validation cases, but they are proving to be quite difficult. Hope this is useful to the discussion.

@FlorianDm
Copy link
Member Author

@WallyMaier oh great I think this book was cited in the publication but I couldn't get access to it

@pcarruscag
Copy link
Member

@FlorianDm sorry I got my timezones mixed up, the meeting usual time is 4pm CET (3pm where I am)

@TobiKattmann
Copy link
Contributor

Unfortunately I am not too much of a help in this, but I hope @WallyMaier 's reference is helpful. I will follow this PR though to see if I can provide s.th. in the process :)

@vdweide
Copy link
Contributor

vdweide commented Nov 5, 2020

Maybe this helps a bit

2D_Axi_Part1

2D_Axi_Part2

@vdweide
Copy link
Contributor

vdweide commented Nov 5, 2020

Although, this seems to be exactly the same reference as @WallyMaier uploaded already....

@FlorianDm
Copy link
Member Author

Hi, thank you all a lot for your help. Here is my derivation :
SourceAxi.pdf

Mistakes have been corrected and it should be the same as your book now apart from the fact that I am neglecting the gradient of viscosity. To take it into account it needs to be evaluated and I suppose that is what those AuxVar variables are for (something like u *mu or u *mu/y I guess). This is new to me and unfortunately I do not have time to look into it atm (unless someone points me in the right direction and it's quickly done)

@WallyMaier I believe that you had taken that into account already and I just came along and cut out your work. So I should close this PR and once you open a new one I can add the generalisation for the convection terms. Were you planning on implementing the k and omega sources as well?

What book is that btw?

@WallyMaier
Copy link
Contributor

@FlorianDm, I had no plans to implement the k-omega terms. I dont think you need to close this PR, I am very happy you took this on! I can add some of the changes I have to this branch. The branch feature_axi has been giving us spurious results and needs some debugging.

As for the book, maybe @vdweide or @economon may know....I only have those two pages haha.

@bigfooted
Copy link
Contributor

@FlorianDm @WallyMaier

I have a different version, but it looks like this is from the book of Hoffmann and Chiang, Computational Fluid Dynamics, volumeII

@FlorianDm
Copy link
Member Author

Aah thanks so it was cited but without the volume so I looked in volume I not knowing there was a volume II

@economon
Copy link
Member

I'm late to the party here, but just a note to say that the original implementation for the incompressible source terms are indeed from the text that @WallyMaier / @vdweide shared. It was added as part of the work in this paper (https://economon.github.io/docs/AIAA-2018-3111.pdf), but I did not test it much beyond a simple laminar channel case or really attempt to treat turbulence at the time. Thanks for putting in more effort on these terms!

@WallyMaier
Copy link
Contributor

WallyMaier commented Nov 15, 2020

@FlorianDm @economon and everyone else

I have merged Florian's derivation with my earlier derivation. I have added the gradients with respect to viscosity. @jtneedels is running some tests on it currently. I can push those changes to this branch so everyone can look over them...what do you all think?

@FlorianDm
Copy link
Member Author

@WallyMaier fantastic! Yes please do push them to this branch

@FlorianDm
Copy link
Member Author

FlorianDm commented Nov 16, 2020

Actually I had a look at your branch and the way it is right now is not correct I believe because I had taken the other variables v and 1/y out of the derivative using the chain rule and combined them with the other derivatives to end up with the terms as they are now so only d(mu)/dy was missing. The AxiAuxVarGrad you are using is apparently d(v*mu/y)/dy so the other terms have to be different. I will change them.

Why not just simply compute the viscosity gradient? Is there any reason not to pull the other variables out?

Is there not already an AuxVar being just v*mu or something?

Anyway I guess it will work the same either way

@FlorianDm
Copy link
Member Author

FlorianDm commented Nov 16, 2020

To do it with your AxiAuxVar i would suggest this:

residual[0] -= 0.0;
residual[1] -= Volume*(yinv*total_viscosity_i*(PrimVar_Grad_i[1][1]+PrimVar_Grad_i[2][0]) 
                      - TWO3*AxiAuxVar_Grad_i[0][0]);
residual[2] -= Volume*(yinv*total_viscosity_i*2*(PrimVar_Grad_i[2][1]-v*yinv)
                      - TWO3*AxiAuxVar_Grad_i[0][1]);
residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1])
                                              + v*(FOUR3*PrimVar_Grad_i[1][1]+TWO3*PrimVar_Grad_i[1][0])
                                              - TWO3*v*v*yinv)
                            - total_conductivity_i*PrimVar_Grad_i[0][1])
                      - TWO3*(AxiAuxVar_Grad_i[1][1]+AxiAuxVar_Grad_i[2][1]));

Could you put that into your branch and then push to this branch? I double checked the signs.. should be fine. If I commited it here I now would break the compilation because these things are not defined in this branch

And for the setters maybe combine the if guards so it's cleaner

      /*--- Set primitive variables for viscous terms and/or generalised source ---*/
      if (!ideal_gas || viscous) numerics->SetPrimitive(nodes->GetPrimitive(iPoint), nodes->GetPrimitive(iPoint));

      /*--- Set secondary variables for generalised source ---*/
      if (!ideal_gas) numerics->SetSecondary(nodes->GetSecondary(iPoint), nodes->GetSecondary(iPoint));

      if (viscous) {
        
        /*--- Set gradient of primitive variables ---*/
        numerics->SetPrimVarGradient(nodes->GetGradient_Primitive(iPoint), nodes->GetGradient_Primitive(iPoint));
      
        /*--- Set gradient of auxillary variables ---*/
        numerics->SetAxiAuxVarGrad(nodes->GetAxiAuxVarGradient(iPoint), nullptr);
      }

@FlorianDm
Copy link
Member Author

In fact i just created a new branch from yours: feature_axi_general_viscous where I made the changes so I will do a new PR from there ok?

@WallyMaier
Copy link
Contributor

@FlorianDm Thanks for doing all that! Unfortunately time zone differences prevented my reply haha.

@FlorianDm
Copy link
Member Author

FlorianDm commented Nov 17, 2020

Yes I was a bit impatient there :)

@WallyMaier WallyMaier deleted the axisym_source_viscous branch December 16, 2020 02:46
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants