[dancfd]

Hello all,

Has anyone run a simulation that has produced good results (e.g. that compare well to published data) on a 2D airfoil? After running a series of simulations with M=0.3 and Re=3E6, I cannot produce useable results.

I would appreciate any assistance anyone could offer.

Thanks,

Dan

[FelixL]

Good evening, Dan,


could you provide some details about what exactly you mean with "good results" and "useable results"?

When it comes to airfoils, it shouldn't pose a problem to get lift and pressure distribution in very good agreement with experimental data using OpenFOAM. Drag is a bit trickier, as this quantity is much more sensitive to the grid, choice of turbulence modeling, spatial descretization and so on.

What exactly is causing you trouble?

Greetings,
Felix.

[dancfd]

Hello Durwin and Felix,

What I mean by good/useable results are data that are in agreement with published results. I am attempting to reproduce McAllister 1982 Experimental Study of Dynamic Stall Vol 2 Page 46. Conditions are alpha = 13.5 deg, NACA0012 airfoil, Re 3.9E6, M 0.301. This yields Cl = 1.40, with a Cp distribution as in the attached image.

My results are very different. I have M and Re similarity, using sonicDyMFoam as a solver, and my case is attached. Unfortunately, I could not attach the polyMesh folder because it is too large. I essentially used the mesh from the wingMotion2D tutorial, with a NACA0012 airfoil. My results are a Cl of 0.07 (using the forces subdict of controlDict) and a Cp plot as attached.

Thanks for taking the time to respond.

Dan
Reference.xcf.zip

case01.zip

mycpplot.xcf.zip

[FelixL]

Hi, Dan,


thanks for posting your case details. Well that's a rather complex setup you're having there. You're using sonicDyMFoam but your mesh isn't moving or changing at the moment, does it? I don't have any experiences with that solver so if the problem results from sonicDyMFoam, I can't really help you there.

But before we reach that point I would suggest you a bottom-up approach to find the possible source. For example I'd start with a steady state, incompressible solver (simpleFoam) to rule out transient effects. If the resulting C_l is reasonable, I'd switch to rhoSimpleFoam and see what happens under the influence of compressibility and so on.

I've seen in your fvSchemes file that you're using Gauss linear on every divergence term of your equations. Didn't you have any stability problems using it? It's an unbounded scheme and requires good care with the mesh quality. Does changing to a limited scheme or first order influence the results significantly?

A look at your checkMesh would also improve the overview of your setup.


Greetings,
Felix.

[dancfd]

Hi Felix,

You're right about sonicDyMFoam. Eventually I will use a scenario with M=0.8 and a moving airfoil, and I wanted to conduct all of the verification and validation on the same solver to ensure consistency. Right now all I need is an incompressible, steady and static solver - but I have had trouble doing this for a number of reasons:

1) For incompressible cases, the units on P are different (incompressible - m2s-2, compressible - kgm-1s-2) not a huge deal, but additional complication. I get the same "is the mesh actually moving?" error as described below with simpleFoam.

2) For static cases, I attempted to run sonicFoam on the same case, but received the following error (log attached). This happens despite deleting the pointDisplacement and dynamicMeshDict files, and commenting out the cellDisplacement subdictionary in fvSolution.

--> FOAM FATAL ERROR: 
mesh flux field does not exists, is the mesh actually moving?

    From function fvMesh::phi()
    in file fvMesh/fvMeshGeometry.C at line 392.

FOAM exiting

3) For steady cases, I cannot decide between rhoSimpleFoam, rhoSonicFoam and rhopSonicFoam - what is the difference? It is of little consequence, since none of these will run anyway. I changed the thermoType in thermoPhysicalProperties to the appropriate value, and added the SIMPLE subdict from the airfoil2D example to the fvSolution dict. I get an error asking for rhoMax, though, and it seems to want a word followed by a dimension set. I cannot determine which word to insert, and there are no tutorials for rhoSimpleFoam.

3a) I have also tried rhoSonicFoam, but I receive an error that I do not understand (log1 attached).

I ran checkMesh, the results are attached (log2). No problems running the simulation, I am trying to run a second-order implicit simulation to allow me to use large timesteps and ensure high accuracy, but the settings I have chosen still restrict me to a Courant number < 1.

As you can see, I have been running into problems at every turn. Do you know how to solve any of them?

Thank you for taking the time to review my case. I really appreciate your help.

Regards,
Dan

[FelixL]

Hey, Dan,

I am sorry about the very late reply. The holidays were keeping me busy.
Here's what I've got to say:


1) Yeah, you're right, the normalized pressure for incompressible cases makes the transition to compressible calculations a bit more difficult. There was a script somewhere around these forums which automatically converted incompressible cases to compressible cases, if I recall correctly. But I don't know if it still works and where exactly to find it. Shouldn't be a problem to write one by yourself, though

2) That error might occur, because you still have BCs for moving walls (if you're using the same BCs as in the case01.zip-file you provided). Change movingWallVelocity in the 0/U-file to fixedValue and see if it works then. Same goes for the simpleFoam-case.

3) rhoSimpleFoam is a steady state solver for compressible flows. It's the only compressible solver using the SIMPLE algorithm! I'd recommend this solver for steady state compressible simulations.

rhoSonicFoam and rhopSonicFoam are transient solvers for compressible flows, the latter one uses pressure based equations. I can't give you any details on the precise differences between these two, but these solvers aren't recommended for steady state simulations. Of course you can use them to march to a steady state with an implicit time scheme, for example, but it will take much more time than with rhoSimpleFoam.

The SIMPLE subdict makes only sense for rhoSimpleFoam. You have to define rhoMax and rhoMin inside this subdictionairy. Something like this:

SIMPLE
{
nNonOrthogonalCorrectors 0;
pMin pMin [ 1 -1 -2 0 0 0 0 ] 0.100;
rhoMax rhoMax [1 -3 0 0 0 0 0] 100;
rhoMin rhoMin [1 -3 0 0 0 0 0] 1e-3;
}

3a) That error should be resolved, when you change the waveTransmissive definition of phi in your pressure BCs to phiv.
You might also wanna read this article: http://openfoamwiki.net/index.php/Ho...dary_condition


Your checkMesh results look good despite of the one poor cell. (Face Skewness > 7!). You can view the face in paraview by using the command

foamToVTK -faceSet skewFaces

and loading the resulting *.vtk-file into paraview. It's better to get rid of that cell, because it might be causing some numerical trouble.

If you still encounter problems, feel free to report here.


Greetings,
Felix.

[dancfd]

Hi Felix,

Thank you for yet another detailed response.

1) As far as I can tell the only change that needs to be made to use a compressible case with an incompressible solver is to the dimensions of the /0/p file. The value of p itself is not required for the calculations since only the gradient is determined, however the value of p/rho vice p must be used in an incompressible case if forces are to be output correctly.

2-3) Thanks, I have now got my case working in simpleFoam and rhoSimpleFoam.

3a) Thanks for the foamToVTK tip; I found the discrepant cell. Unfortunately this brings up another issue whereby I cannot get OF to produce the mesh motion I am prescribing without errors. In this case, it was a simple rotation of 13.5 degrees (angularOscillatingDisplacement BC, inverseDistance diffusivity). I have found a way around this for the meantime (create the mesh around the already-rotated geometry) but for my dynamic simulations it looks like I will have to code a new diffusivity/mesh motion solver. The RBF solver in 1.5-dev works, however it is too computationally expensive.

Unfortunately my results remain incorrect (Cl ~ 0.56 vice 1.4). I believe the source of the issue now lies with turbulence modeling. I am uncertain how to select values for k and omega in the komegaSST model, since Menter's 1994 paper gives a range for both k infinity and omega infinity.

I think the turbulence and mesh motion are the last nuts to crack. I cannot thank you enough for your assistance.

Best regards,
Dan

[FelixL]

Hi, Dan,


glad I could be of assistance.

I am wondering why the lift coefficient is still over 50% off. Lift over an airfoil is mainly pressure driven and turbulence (or friction in general) plays a minor role, as long as the BL remains attached. I was expecting that at least the simpleFoam computation even with 1st order discretization leads to more reasonable results, but as it seems this is not the case.

Well anyways, I am sure you will obtain better results and be successful in the end. If there still are any questions coming up (e.g. regarding kOmegaSST), feel free to ask.

Good Luck!

Greetings,
Felix.

[jms]

Hi Felix!

I am running a computation to get the characteristics of a NACA0015 at RE=2x10^6 using OpenFOAM. The results I get are pretty good for the lift coefficient (error of 10% at most) and very bad for the drag(I get a drag 4 times bigger than the one I should). This is why I want to explain the BC applied to see if you can help me.
I am using the Spalart Allmaras turbulent model with a value for nuTilda and nut of 0.006 for the initial value on the inlet, outlet and the internal field, coming from the equation ((U*I*l)*sqrt(3/2)). For nut I am using the "nutSpalartAllmarasWallfunction" with a value of 1e-10 (based on what I found on the forum). I read that Spalart Allmaras is a hybrid model so this should work fine for different meshes. I tried for 1st cell height of 10^-6, 10^-5, 10^-3 and I get wrong results for all of them, even though y+ looks ok. However, the drag is better for the case of 10^-3 (but still 3-4 times the value I should get), getting a worse result for the lift coefficient on this case...

Thank you very much for your attention.
I hope you can give some ideas because I am struggling a lot with it.

José

[FelixL]

Hello, José,


I need some more details before I can help you out.

First of all: does checkMesh tell that everything is okay with your mesh?
What is the order of accuracy of the discretization for the divergence terms? The numerical dissipation of first order accurate schemes usually leads to much higher drag coefficients.

I am assuming you are simulatin external air flow with a kinematic viscosity of roundabout 1.5e-5 m²/s? What's the freestream turbulence characteristics like? Low turbulence assuming fully turbulent boundary layers on the airfoil?

What angle of attack are you simulating?


Greetings,
Felix.

[jms]

Hi Felix!

First of all thank you for your quick asnwer!
checkMesh complains about the AR. It says that I have a maximum aspect ratio of 2863.93 and it considers I have 3581 cells with high AR (a 4% of the total amount of cells). But I think this is not my problem since for the case I run for 1st cell height of 10^-3, it does not complain at all.
I am using the following divergence schemes:
divSchemes
{
default none;
div(phi,U) Gauss linearUpwind Gauss linear;
div(phi,nuTilda) Gauss linearUpwind Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}
should I change them ?

I am simulating for a range from 0 to 16 degrees. But the one I am doing more trials on is for AOA=8degrees. The kinematic viscosity I use is 5e-7, so for a velocity of 1m/s I get a RE=2milions.
For the turbulence...I was expecting to have, for those values of nut and nutilda I wrote above, viscous and turbulent boundary layer (where I used I=5% and length scale l=0.07*L=0.07, for a chord length of 1).

Thank you very much for your attention.
I really appreciate all the suggestions you will give to me.

José

[FelixL]

Hello, José,


the AR values shouldn't do no harm, according to my experience with SpalartAllmaras. Your div-schemes are fine.

The alpha=8° case is the one you described in your previous post? Or do you encounter much too high values of drag over the whole range of alpha?

For the inlet values of nuTilda I'd recommend to you to read [1]. The SpalartAllmaras turbulence model implemented in OpenFOAM doesn't account for transition between laminar and turbulent boundary layers.
Finding exact values of nuTilda at the inlet isn't easy. For a start I would use a value of , as suggested in [1].

This should result in fully turbulent boundary layers everywhere and hopefully improves your drag results.

Greetings,
Felix.


[1] Effective Inflow Conditions for Turbulence Models in Aerodynamic Calculations, P. Spalart and C.L. Rumsey, AIAA Journal, vol. 45, issue 10, pp. 2544-2553

[fletc900]

Originally Posted by FelixL

Hello, José,


the AR values shouldn't do no harm, according to my experience with SpalartAllmaras. Your div-schemes are fine.

The alpha=8° case is the one you described in your previous post? Or do you encounter much too high values of drag over the whole range of alpha?

For the inlet values of nuTilda I'd recommend to you to read [1]. The SpalartAllmaras turbulence model implemented in OpenFOAM doesn't account for transition between laminar and turbulent boundary layers.
Finding exact values of nuTilda at the inlet isn't easy. For a start I would use a value of , as suggested in [1].

This should result in fully turbulent boundary layers everywhere and hopefully improves your drag results.

Greetings,
Felix.


[1] Effective Inflow Conditions for Turbulence Models in Aerodynamic Calculations, P. Spalart and C.L. Rumsey, AIAA Journal, vol. 45, issue 10, pp. 2544-2553

Felix,

Did you mean nuTilda = 5*nu = 5*1.4607e-5 = 7.3035e-5 m^2/s rather than the 2.5e-6 m^2/s value? I'm assuming you guys are talking about air here (sea level)...

Also is it common convention of setting nut = nuTilda? Looking at the airfoil2D example for the SimpleFOAM I see nut and nuTilda being set with the same values for both initial and boundary conditions...

Jose,

The 1e-10 value for nut seems to be very low. The turbulent viscosity ratio (nut/nu) for external flow is normally in the range between 1 and 10. For air at sea level, that would put nut in the range of 1.4607e-5 to 1.4607e-4.

Regards,
Bruno

[FelixL]

Hello, Bruno,


José said in his post above, that he's using a value of 5e-7 m²/s for the kinematic viscosity.

Setting nuTilda=nut is only reasonable, if , which is not the case here, yeah.
You can calculate nut using the formula:

(see http://turbmodels.larc.nasa.gov/spalart.html )

Setting the boundary of nut to calculated at the inlet makes OpenFOAM doing the conversion automatically. The initial values of nut will be corrected before the first iteration step so those can be set arbitrarily, as far as I know.

And regarding the value of 1e-10 for nut: José meant the near-wall value. Essentially nut is zero on the surface of a no-slip wall. This is because of the no-slip condition, turbulent fluctutations cannot occur and so the turbulent viscosity is damped in the viscous sublayer. Setting nut to 0 might cause some numerical troubles (division by zero), so it's often set to a very small, nonzero value at the wall.
But since José is using wallFunctions anyways, the nut-value at the wall will be ignored either way, so this value is just a dummy value and doesn't affect the solution at all.


Greetings,
Felix.

[jms]

Hi!
Thank you again for all you suggestions! They really work!

But, there is 1 point where I am not very confident yet though...
What about the internalField for both nut and Nutilda? which value do I set to it? For the outlet I just set zeroGradient and for the inlet and the wall what you just told me.

And I would like to ask you also about this "freestream" condition I use (because I saw it on the tutorial), what does it exactly do? I set it for both the inlet and the outlet. The other option would be setting a fixedValue at the inlet and a zeroGradient at the outlet, right? which is the difference between these 2 options?

Again...thanks! this really helps me!

Greetings,

José

[FelixL]

Hello, José,


choose the inlet value as your initial internal values of nuTilda. This is common practice when there's no other data available. As for nut, those initial values can be set arbitrarily. Choose the inlet value of nut, if you have one. If nut is set to calculated at the inlet, choose 1 or anything else. the internal field of nut will be automatically corrected before the first iteration step.

I don't know for sure, what the freestream BC does, because I never used it before. Search the forums or have a look at the source code, that might give you an idea. Personally, I prefer to stick to standard BCs whenever I can, i.e. using fixedValue for p and zeroGradient for the other quantities at the outlet and vice versa at the inlet.


Greetings,
Felix.

[fletc900]

Felix,

Thanks for the quick response and for the clarifications. Your explanation and the link that you have sent for the SA model are going to be very useful to me...

Best regards,
Bruno

[jms]

Hello Felix,

Thank you for your recommendations. They worked better than before!

Now I am doing the same simulation on this airfoil but with the k-Omega SST model.
Before to continue I want to make sure some aspects related with this 2 models (k-omega SST and Spalart Allmaras being low-RE models, High-RE models or hybrid models...). It has already been explained on other threads but I have read so many things that I prefer asking it straight!

When you implement wall functions they work as hybrid models both of them, don´t they? (nutSpallartAllmarasWallfunction for Spalart Allamaras and KqrWallFunction/nutWallFunction/omegaWallFucntion for K-OmegaSST). For which range of y+ they work well?
The meaning of this constant you define just after this wall function what is it? I read on one thread that it is one approach for the result you will get, is it like that?
If it is, and correct me if I am wrong please, I would put all these values for the different wallFunctions to 0 (or 1e-10 to avoid dividing by 0) but the value of omegaWallFunction, which, based on what I found on http://turbmodels.larc.nasa.gov/sst.html is something based on this omega_wall), right.

I am really looking forward to reading your answer because it may make things clear due to different answers found, as I said above.

Thanks,

Greetings!

José

[FelixL]

Good morning, José,


if I recall correctly, with OpenFOAM 1.7 nutWallFunction and nutSpallartAllmarasWallfunction are the same now. Both wall functions use Spalding's law, which is suitable for any range of y+, from the viscous sublayer up to the log layer. Same goes for KqrWallFunction, which actually is simply imposing a zeroGradient BC for k at the wall.

omegaWallFunction uses a blending method between viscous sublayer and log layer values to obtain values for omega in the buffer layer, so it's also suitable for the whole inner region of a turbulent boundary layer.

I didn't use any of these wall functions (except for omegaWallFunction) for meshes with y+ < 30, because I like to have control over the near wall values when using a mesh, where the viscous sublayer is resolved. So my experience with these continuous wall functions is quite limited.

The meaning of the constant after the wallFunction definition is nothing else than a dummy value. I think ParaView needs this value to initialize fields, but this dummy value doesn't affect the simulation at all.


If you wanna use kOmegaSST on lowRe meshes (which of course is possible), I'd recommend to you to use these BCs at the walls:

• fixedValue 0 or 1e-10 for nut and k
• omegaWallFunction for omega

It's better to use omegaWallFunction, because omega_wall depends upon the first cell height and omegaWallFunction automatically uses the equation for the viscous sublayer, if your y+ is low. That saves you some time.

You've got to be aware, though, that the original implementation of kOmegaSST in OF 1.7 leads to drag values, which are slightly too high. See this thread to see how to correct this error: http://www.cfd-online.com/Forums/ope...nce-model.html


Greetings,
Felix.

[jms]

Dear Felix,

I am actually using OpenFOAM v1.6 but I hope your explanations given for OpenFOAM v1.7 will be useful.
I am running the simulations using high-Re meshes and what I want to tell you is that this constant after the wallFunction definition is not a "dummy value", at least for omegaWallFunction, since modifying it, changes my results. I used the equation for this value

10*6*nu/(beta1*(h1)^2),

where h1 is the 1st cell height and beta1 is one constant as defined on page http://turbmodels.larc.nasa.gov/sst.html on "w_wall". I am wondering if this value for "w_wall" is referring to low-Re meshes ...because I don´t get the desired results.

For all the parameters defined on the folder "0/" the only values of the initial conditions that affect the results are the ones at the inlet, aren´t they? the ones of the internalField do not affect the results at all, do they? (of course giving reasonable values to start the simulation...)

If you can give me any help, again...it will be very welcome!

Thanks,

Greetings,

José