[nuovodna]

I'm working on a NACA0012 in 2D and i must calculate the lift and drag coefficent. Can someone post here a complete liftDrag tool with all files needed to compile and install it? Not only the files under directory src/postProcessing/incompressible/liftDrag but also the files under application/utilities/postProcessing/miscellaneous/liftDrag
.Can someone explain me how the tool calculate the lift and drag coefficents? In the forum there are many thread on this argument but no one is complete.

Thanks a lot

Emanuele

[msrinath80]

Everything that you've asked for (and more) is in this board already. Please use the search utility to the left to find it.

[alberto]

I think the question can be reformulated as "How can I compile liftDrag on OpenFOAM 1.4 and following?" ;-)

Regards,
Alberto

[gtg627e]

Dear Emanuele,

I have been working on a naca0012 as well.
These are the things I found out so far:

1. The solution is very very sensitive to your mesh quality. At the moment, I am using a structured mesh with maximum non-orthogonality of 9 deg, average 4 deg. Maximum aspect ratio 808. Max skewness 0.28. You can find out this information by running checkMesh.

2. You should initialize your case with potential foam to get good initial conditions for velocity. Copy potentialFoam/yourcase/0/U into simpleFoam/yourCase/0

3. At this point you should start your simulation with turbulence off (simpleFoam/yourCase/constant/tubulenceProperties)
Let the pressure and velocity initial residuals drop to say 1e-3, and then turn the turbulence on.
Remember that to be able to have your code see changes as it is running you need to have runTimeModifiable yes in simpleFoam/yourcase/system/controlDict.

4. You should converge your results in simpleFoam with very conservative schemes, i.e.:

file simpleFoam/yourCase/system/fvSchemes

divSchemes
{
default none;
div(phi,U) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,omega) Gauss limitedLinear 1;
div(phi,R) Gauss upwind;
div(R) Gauss linear;
div(phi,nuTilda) Gauss upwind;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

5. After your solution has converged to 1e-6 (residual for all variables), you can go more aggressive with the discretization schemes. for example:

file simpleFoam/yourCase/system/fvSchemes

divSchemes
{
default none;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,omega) Gauss limitedLinear 1;
div(phi,R) Gauss linear;
div(R) Gauss linear;
div(phi,nuTilda) Gauss linear;
div((nuEff*dev(grad(U).T()))) Gauss linear;
}

If you compute your drag at steps 4 and 5, you may notice a difference of 20%. This is due to the better results you get with the more aggressive discretization schemes.

6. Your drag estimate is going to be very sensitive to your boundary conditions at the airfoil also. At the moment I am having very satisfactory results using nutWallFunctions for the airfoil with a y+ of about 7. You can check your y+ during your solution or when the solution is converged issuing checkYplus.

I hope this is a good start.
Keep me posted about your findings,

Alessandro

[nuovodna]

Good morning Alessandro! I have followed your advices but i have a problem with the boundary conditions: where should i set the nutWallFunctions ?

Thanks in advance

Emanuele

[nuovodna]

what turbulence model do you use?

[gtg627e]

Emanuele,

At the moment, I am experimenting with:

1. kOmegaSST (High-Re) meaning you use wall functions.


2. LaunderSharmaKe (Low-Re) meaning you do not use wall functions.

To set wall functions, edit /simpleFoam/yourCase/constant/polymesh/boundary.

For case 1:

frontAndBack
{
type empty;
...
...
}

inlet
{
type patch;
physicalType inlet;
...
...
}

outlet
{
type patch
physicalType outlet;
}

airfoil
{
type wall;
physicalType nutWallFunctions;
(or physicalType wallFunctions);
...
...
}

/simpleFoam/yourCase/0/p

frontAnfback
{
type empty;
}

inlet
{
type zeroGradient;
}

outlet
{
type fixedValue;
value uniform 0;
}

airfoil
{
type zeroGradient;
}

/simpleFoam/yourCase/0/U

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform (50 0 0)
}

outlet
{
type zeroGradient;
}

airfoil
{
type fixedValue;
value uniform (0 0 0)
}

/simpleFoam/yourCase/0/k

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform 0.375;
}

outlet
{
type zeroGradient;
}

airfoil
{
type zeroGradient;
}


/simpleFoam/yourCase/0/omega

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform 12800;
}

outlet
{
type zeroGradient;
}

airfoil
{
type zeroGradient;
}


For case 2, I am not so sure about how to enforce no wall functions but I think it is like this:

frontAndBack
{
type empty;
...
...
}

inlet
{
type patch;
physicalType inlet;
...
...
}

outlet
{
type patch
physicalType outlet;
}

airfoil
{
type wall;
physicalType wall;
...
...
}

/simpleFoam/yourCase/0/p

frontAnfback
{
type empty;
}

inlet
{
type zeroGradient;
}

outlet
{
type fixedValue;
value uniform 0;
}

airfoil
{
type zeroGradient;
}

/simpleFoam/yourCase/0/U

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform (50 0 0)
}

outlet
{
type zeroGradient;
}

airfoil
{
type fixedValue;
value uniform (0 0 0)
}

/simpleFoam/yourCase/0/k

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform 0.375;
}

outlet
{
type zeroGradient;
}

airfoil
{
type fixedValue;
value uniform 1.0e-10;
}


/simpleFoam/yourCase/0/omega or epsilon (since I use LaunderSharmaKE)

frontAnfback
{
type empty;
}

inlet
{
type fixedValue;
value uniform 12800; Make sure this value is consistent with your analysis
}

outlet
{
type zeroGradient;
}

airfoil
{
type fixedValue;
value uniform 1.0e-10;
}

I hope this helps,

Alessandro

[ariorus]

Hello,

I thought physicalType was only used in FoamX and that it didn't have any effect on solvers.

Wall functions should be used only if an high Reynolds turbulent model is employed otherwise they are not used, regardless of what physicalType" one is using.

Please correct me if I'm wrong or something has been changed in openfoam about this.

Ciao.

[nuovodna]

this is the output of checkYPlus:

Create mesh for time = 1000

Time = 1000
Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model kOmegaSST
Patch 3 named ala y+ : min: 0.0658652 max: 102.546 average: 65.7859

End

[gtg627e]

Emanuele,

The inlet velocity seems very low. What Re are you trying to obtain?
What is your kinematic viscosity?

Also, remember that in liftDrag a value of Aref is computed, which is equivalent to the wetted area of your wing, but your are interested in your
sectional drag coefficient I assume. This is what I do:

1. Set Aref=1 in src/postProcessing/incompressible/liftDrag/liftDrag.C
recomplile your liftdrag tool by issuing ./Allwmake.
Then divide the drag coefficient that you obtain this way by your depth in the z-direction.

Make sure that in applications/postProcessing/miscellaneous/liftDrag/liftdrag.C
turbDragCoeffcient is called, instead of the default dragCoefficient. If you change the function call I mentioned, recompile by issuing wclean && wmake.

As an after thought, does your velocity field look reasonable?
What about your pressure field?

From "Theory of wing sections", you get the maximum velocity over a naca0012 at 0-deg aoa as v_local/Vinf = 1.188. Does your velocity field agree with this?


Let me know,

Alessandro

[nuovodna]

Hi Alessandro,
in transportProperties file there is this line:

nu nu [0 2 -1 0 0 0 0] 1.4e-06;

I'm trying to obtain Re = 1e+06


If my depth in the z-direction is 1, can i unchange the liftDrag.C file?

The velocity field looks reasonable in paraView also the pressure field.

[gtg627e]

Emanuele,

What is the chord length of your airfoil?
To get Re=1e6, given your nu, it looks like the chord is ~ 0.72.

So, the aspect ratio of your faces in the z-direction is very high. If you search the postings, you will find a posting where Dr. Jasak mentioned the procedure that is used to compute the drag coefficient; he also mentioned that while making your z-depth = 1 is appealing, since you don't have do recompute your Cd (by normalizing about z-direction, if Aref = 1), this leads to significant errors.

I think you should reduce your z-dimension to something like 1%-10% of your chord, to be safe.

Furthermore, I do not know the details of the Navier-stokes formulation in OpenFOAM, but I remember reading that solving the Navier-Stokes equations at really low speed may be difficult (I think there are special formulations for this, but look this up I am not a CFD expert). So I would suggest changing the inlet velocity to something reasonable like 50, and playing with nu to get your desired Re.

Alessandro

[gtg627e]

Emanuele,

Sorry about the chord length;
This is what I mean Re=chord*U/nu; I messed that up in my previous posting.
I see your chord must be 1, right?

Alessandro

[nuovodna]

i have change my z-depth to 1% of chord (my chord's lenght is 1 ---> z-depth 0.01) and i've adjusted nu to obtain Re = 1e6 with a velocity inlet of (50*cos(18°), 50*sin(18°), 0) supposing an aoa = 18 degree but the lift coefficent is too low.

[gtg627e]

Emanuele,

With aoa = 18 deg, you have massive separation. SimpleFoam, being a steady solver, is not applicable to this case. You should use turbFoam.

For a naca0012, I think you are safe using simpleFoam up to 8-9 deg maximum.

Alessandro

[nuovodna]

Alessandro,

with simpleFoam now i obtain a reasonable value of lift coefficent from 0 to 6 deg. I'm traying to use turbFoam for aoa > 6 deg but i have always the floating point error. In what way can i solve this error??
Thanks a lot

Emanuele

[gtg627e]

Emanuele,

Could you be more specific about your turbFoam case:

1. turbulence model
2. checkMesh output and y+
3. Boundary conditions
4. Initial conditions
5. deltaT and Courant number you observe during the simulation.
6. fvSchemes and fvSolution specifications.

Alessandro

[gtg627e]

Emanuele,

I don't know whether or not the spalartAllmaras uses wall functions. If it does, you should set zeroGradient for all turbulence variables at the boundary patch ala like you have done. If spalartAllmaras does not use wall functions, you should set all turbulence variables to 1e-10 at the patch ala. I think there is something wrong since your y+ is 0 everywhere.

How do you specify nuTilda?
I think you need nuTilda for the spalartAllmaras turbulence model. In order to find out which variables you need to specify for the SpalartAllmaras turbulence model edit /src/turbulenceModels/incompressible/SpalartAllmaras/SpalartAllmaras.C

At line 125 you can see that this turbulence model definitely wants to read nuTilda, so you need to specify it. for external flows, (please double check what I am saying) nuTilda = beta * nu, where 1 <= beta <= 10.
For more details see:

http://www.cfd-online.com/OpenFOAM_D...es/1/3476.html

fvSchemes and fvSolutions look good. Later on when you get your model working you might what to tighten up convergence tolerances to at least 1e-6.

I am no expert in CFD and much less in turbulence modeling, but are you sure you can use spalartAllmaras for separated flows? (I assume you are still trying to analyze your naca0012 at 18-deg aoa).

If you are still trying to simulate a naca0012 at 18-deg aoa, I think you would be better off with kOmegaSST.

I assume you farfield1 patch is the inlet correct?

I hope this helps,

Alessandro

[ariorus]

Hello,

yPlus are equal to

yPlus_* =Cmu^0.25 sqrt(k) yp / nu

so for Spalart Allmaras the field is always zero since in that model k is not used and its value is set to zero.

On the other hand I think it could be useful to have the "other" yPlus (without *):

yPlus= u_* yp / nu

yPlus_* and yPlus should be comparable in equilibrium turbulent boundary layers.

I'm attaching here an application which computes both yPlus* and yPlus and writes them on files so that the fields can be visualized by paraFoam. I hope it can be useful to anybody, and/or that anybody will correct it if I made some mistakes ;-) .

computeYplus.tgz

Just one more thing:
in the liftDrag utility the computation of the turbulent force is mesleading it is better to compute laminar and turbulent forces at the same time:

laminarTurbForce += gSum ( -turbulence->nuEff()().boundaryField()[patchI]* U.boundaryField()[patchI].snGrad()* mesh.magSf().boundaryField()[patchI] );

This explains why in the automotive test case thread I mentioned strange values obtained for turbulent forces.


Ciao.

[gtg627e]

Rosario,

Thank you for your elucidations on the Spalart-Allmaras model, and the associated method to compute y+.

Now I have a question for you. I have been going crazy with trying to understand wallFunctions. I now understand that there are Hi-Re models that do use wallFunctions such as kOmega, kEpsilon on so forth. I experimented with LaunderSharma (no wallFunctions) but I get a really bad looking velocity field at the leading edge of my naca0012.

I went through all threads discussing wallFunctions, and the only thing I saw was people enforcing the use of them with physicalType ....

As you have pointed out though, physicalType specifications are only used by FoamX. Do hi-Re turbulence models automatically use wallFunctions, no matter what you specify?

what about low-Re turbulence models? Do you set you wall patch as type wall?

Thank you,

Alessandro