ABL in OF2.0.x
 

Hi,
I'm testing the ABL in OF 2.0.X.
I build a flat "2D" mesh with grading in Z.
In the picture attached you can see the velocity probed at 5 m from the inlet.
I found a very strange jump between the log function and the free-stream.
Can anyone help me to understand it? (I attached also the case)
Have anyone found problem with this BC?
Or maybe is there a bug?

Tahnks

Andrea

Hello Andrea,

I have also struggled with this problem for quite some time.

I can probably guess what the problem is, but I am not sure. I can't seem to open your case.

The problem I think it is, is that in OF-2.0 the code uses 2 different equations for ABL. If you look at the code for ABL you can see that it uses all the parameters dictated in the ABLcondition file except Ustar. It uses all these parameters to calculate this Ustar with an equation Ustar=U_ref*kappa/log((Z_ref+Z0_)/min(Z0,0.001) So it will take the minimum of Z0 or 0.001.

A little bit further down the code it calculates the ABL inlet wind speed with equation Ustar=U_ref*kappa/log((Z_ref+Z0_)/(Z0). So there will be a shift in your Ustar and this can be the cause for your problem.

This difference in equation is new in OF-2.0, in OF-1.7.x there was no difference. I don't know the reason yet. You can solve this in two ways, you can adjust the code of you can take into account the difference when you choose the parameters.

I hope this helps.

I also have some problems with the nutkroughwallfunction for the turbineSiting case where this ABL is used. If someone reads this tread, maybe they can join the discussion. I am trying to simulate a flow over a flat plane and try to have always the same wind profile using the ABL condition and the wallfunction.

Regards,

Jochem

Hello,

I'm using ABL to calculate the flow over a terrain. Now I'm doing the mesh with snappyHexMesh options, but still don't have a good surface inside the box yet.

At this moment I have a stupid question. I read but still it isn't so clear for me. So, here is the question:
How can you calculate
ABLConditions: Ustar, Uref, Href, z0, zGround,turbulentKE; and
initialConditions: flowVelocity, pressure, turbulentKE, turbulentEpsilon?

Best Regards,
Ovidiu

Could you replace the min(z0, 0.001) by max(z0,0.001)? This approach is traditionally used to prevent divide by zero problem. I think that it could like be a mistake.

Also the parenthesis do not match on the expression: Ustar=U_ref*kappa/log((Z_ref+Z0_)/(Z0).

Hi Ovidi,

zGround is the lowest part in your domain, this is the easy one. The Abl condition and initial conditions can be calculated from winddata, I suspect that is the right way.

I use the simplified loglaw to calculate a surface roughness with 2 datapoints (same point of time, different hight). With this data I calculate Ustar, and fill in the surface roughness with z0. H_ref is the place where the free stream velocity U_ref is reached. H_ref is in my case the maximum height of the domain, because I suspect that the freestream velocity isn't reached within the height of my case. I use the formula to calculate U_ref.

Flow velocity and pressure can remain the same is in the tutorial. For turbulentKE and turbulent epsilon I have found a document on the forum, I don't know where so I've attached the document. You can see the formulas to calculate the constants for the initial condition.

Julien,

The way you've suggest is also a good way to prevent this mistake. And you are correct, I've made a mistake with the parenthesis.

If anyone wants to make additions or corrections to my explanation, please do so.

Regards,

Jochem

Hi Jochem,

Thank you for useful information.
Also, can you attache again the document which calculates TurbulentKE and turbulent epsilon. Thks.

Best Regards,
Ovidiu,

Dear All,
I solved it using the 17x version where the ABL bc is ok.
For the zGround I think is better to use the lowest z of the inlet patch (Z_low_inlet).
If you have a model with a global lowest Z (Z_low_global) below the lowest Z inlet (Z_low_inlet) and you choice the first Z as zGround = Z_low_global, you make a mistake in your inlet log profile.
I think another problem is when your Z_low_inlet changes in Y direction.
Did anyone found this problem? Is in OpenFOAM something to take in account this?
I think the fast way to overcome this, it's to separate the inlet patch in different patches @ different Z_low_inlet.

Regards

Hi,

I'm studding the wind over a complex terrain. I toke like example the simpleWindFoam tutorial. I start the simulations and I'm receiving this error after a few iterations:

#0 Foam::error::printStack(Foam::Ostream&) in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#1 Foam::sigFpe::sigFpeHandler(int) in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libOpenFOAM.so"
#4 void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#5 Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::operator/<Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<doub le, Foam::fvPatchField, Foam::volMesh> > const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#6 Foam::incompressible::RASModels::kEpsilon::correct () in "/home/openfoam/OpenFOAM/OpenFOAM-1.7.0/lib/linuxGccDPOpt/libincompressibleRASModels.so"
#7 main in "/home/openfoam/OpenFOAM/openfoam-1.7.0/applications/bin/linuxGccDPOpt/simpleFoam"
#8 __libc_start_main in "/lib/libc.so.6"
#9 _start at /usr/src/packages/BUILD/glibc-2.10.1/csu/../sysdeps/i386/elf/start.S:122
Floating point exception


If you receive something similar can you explain why I'm receiving this kind of error?
Thank you in advance,

Best Regards,
Ovidiu

If I don't remember bad, int the tutorial you have in the file k:

value uniform 0.0;

So you divide per 0 and got floating point.
Change it as:

value $internalField;



Regards

Hi,

Thank you guys for the information. It's working.

Hello,

I'm running a simulation with a wind passing on a terrain. The solution converge and I want to plot the velocity profile. I manage to do this but in excel (see the picture) and i want to plot in ParaView. How can I do it?

Thank you

Hello,

You have to use the filter->Data Analysis -> plot over line. Type in the coordinates you want to use.

Regards,

Jochem

Hi to everyone, I'm simulating the flow throw a delta wing in open sky, I need to generate the ABL but it doesn't work! I'm using the OF2.1.1 maybe that's the problem.

Please Ovidiu can you shre your 0 and patch files for taking a look at them?

thanks!

Hi there I am trying to run a different terrain in turbineSiting tutorial.But I ran into problems with that so I went back to the tutorial case and tried to edit the wind direction.I think my mesh is fine since it passes the checkmesh.The strange thing is when I set the ABL conditions.

Uref 8.16;
Href 70;
z0 uniform 0.1;
turbulentKE 1.3;
windDirection (1 0 0);
zDirection (0 0 1);
zGround uniform 0.0;

The problem is simple like in the Turbinesiting tutorial.From the west there is an inlet from east there is an outlet the terrain is defined like a wall and top and sides are no slip

With these initial Condition the solution converges in 66 iterations.So when I change the wind direction to (1 1 0) the solution crashes.But if I use (5 2 0) the solution still converges.

Create time

Create mesh for time = 1

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting RAS turbulence model kEpsilon
bounding k, min: 0 max: 3.42645150226 average: 1.30628839732
kEpsilonCoeffs
{
Cmu 0.09;
C1 1.44;
C2 1.92;
C3 -0.33;
sigmak 1;
sigmaEps 1.11;
Prt 1;
}

Creating field source list from sourcesProperties


SIMPLE: convergence criteria
field p tolerance 0.01
field U tolerance 0.001
field "(k|epsilon)" tolerance 0.001


Starting time loop

Time = 2

smoothSolver: Solving for Ux, Initial residual = 0.473678546445, Final residual = 0.0253675867531, No Iterations 3
smoothSolver: Solving for Uy, Initial residual = 0.545508812572, Final residual = 0.0410704803125, No Iterations 3
smoothSolver: Solving for Uz, Initial residual = 0.490732100575, Final residual = 0.0369494371199, No Iterations 3
GAMG: Solving for p, Initial residual = 0.235186667396, Final residual = 0.0115677539951, No Iterations 4
time step continuity errors : sum local = 0.000123942210288, global = -1.07094754161e-05, cumulative = -1.07094754161e-05
smoothSolver: Solving for epsilon, Initial residual = 0.0142418185171, Final residual = 0.00129445182949, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.443069945282, Final residual = 0.0234141579111, No Iterations 3
ExecutionTime = 1.13 s ClockTime = 1 s



Time = 3

smoothSolver: Solving for Ux, Initial residual = 0.46653783383, Final residual = 0.0263371603412, No Iterations 3
smoothSolver: Solving for Uy, Initial residual = 0.309793353435, Final residual = 0.0170232186791, No Iterations 3
smoothSolver: Solving for Uz, Initial residual = 0.205484283561, Final residual = 0.0129615975667, No Iterations 3
GAMG: Solving for p, Initial residual = 0.16854903988, Final residual = 0.0155268143013, No Iterations 3
time step continuity errors : sum local = 0.000210952976415, global = 1.76511520803e-05, cumulative = 6.94167666417e-06
smoothSolver: Solving for epsilon, Initial residual = 0.00969959500694, Final residual = 0.000656979811585, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.344065799334, Final residual = 0.021449332508, No Iterations 3
ExecutionTime = 1.78 s ClockTime = 1 s

Time = 4

smoothSolver: Solving for Ux, Initial residual = 0.149131904737, Final residual = 0.00828018673637, No Iterations 3
smoothSolver: Solving for Uy, Initial residual = 0.186625072917, Final residual = 0.0100424579985, No Iterations 3
smoothSolver: Solving for Uz, Initial residual = 0.139523403237, Final residual = 0.00937017636232, No Iterations 3
GAMG: Solving for p, Initial residual = 0.20638215344, Final residual = 0.0181687040399, No Iterations 5
time step continuity errors : sum local = 0.000109995029081, global = -1.4898546506e-05, cumulative = -7.95686984185e-06
smoothSolver: Solving for epsilon, Initial residual = 0.00758607525273, Final residual = 0.000583556877572, No Iterations 2
smoothSolver: Solving for k, Initial residual = 0.215939251023, Final residual = 0.0141821403433, No Iterations 3
bounding k, min: 0 max: 3.32130962682 average: 1.32598953792
ExecutionTime = 2.47 s ClockTime = 2 s

Time = 5

smoothSolver: Solving for Ux, Initial residual = 0.231364300243, Final residual = 0.0131708756433, No Iterations 3
smoothSolver: Solving for Uy, Initial residual = 0.103479618454, Final residual = 0.00572050956008, No Iterations 3
smoothSolver: Solving for Uz, Initial residual = 0.278417864979, Final residual = 0.0194834887617, No Iterations 3
GAMG: Solving for p, Initial residual = 0.269803139772, Final residual = 0.0259705324425, No Iterations 2
time step continuity errors : sum local = 0.000175945244605, global = 5.86352385936e-05, cumulative = 5.06783687518e-05
#0 Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#2 Uninterpreted:
#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"
#4 void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"
#5 at kEpsilon.C:0
#6 Foam::incompressible::RASModels::kEpsilon::correct () in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"
#7
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/simpleFoam"
#8 __libc_start_main in "/lib/i386-linux-gnu/libc.so.6"
#9
in "/opt/openfoam211/platforms/linuxGccDPOpt/bin/simpleFoam"
Floating point exception (core dumped)


As you can see there is a bounding at k just before it crashes.The last iteration I persume says there is a problem at Kepsilon model.Floating point expection usually means division by zero I figured that out but simply all the parameters are same and the solution just crashes.

I am using OpenFoam version 2.1 and I am using exactly the same wall conditions and solvers as proposed by the tutorial.

for Epsilon : AtmBoundaryLayerInletEpsilon and EpsilonWallFunction

for k: kqrWallFunction

for nut: nutkAtmRoughWallFunction

is used.

I also checked the previous errors where you had to change min(Zo,0.001) to max(Zo,001) but that doesnt exist.I also tried to change the relaxation values;no use.

I won't post my solution and other sheets since they are exactly the same with the tutorial ,since I didn't change anything.The terrain is also AskervereinHill.

If somebody could give me some insight that would be terrific

My Regards
Burak

Well I did some more digging and operation about my problem;

#2 Uninterpreted:


#3 Foam::divide(Foam::Field<double>&, Foam::UList<double> const&, Foam::UList<double> const&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libOpenFOAM.so"


#4 void Foam::divide<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"


#5 at kEpsilon.C:0


#6 Foam::incompressible::RASModels::kEpsilon::correct () in "/opt/openfoam211/platforms/linuxGccDPOpt/lib/libincompressibleRASModels.so"


Floating point exception (core dumped)



The highlighted parts indeed meant that in the model K or epsilon value was getting the value of "zero" and division was causing a floating point.




For my case for the modified "turbine siting tutorial" I changed the values of k at case/0/k where the highlighted parts were changed from 0.0 to 1.3 which is the internal field value that $turbulentKE calls from "initial conditions" .But I guess something other than zero would work since considering what I read it gets revised and updated no matter what anyway.


At the end it looks like this;


#include "include/initialConditions"

dimensions [0 2 -2 0 0 0 0];

internalField uniform $turbulentKE;

boundaryField
{
#include "include/ABLConditions"

outlet
{
type inletOutlet;
inletValue uniform 1.3;
value $internalField;
}
inlet
{
type uniformFixedValue;
value $turbulentKE;
uniformValue constant $turbulentKE;
}
"terrain_.*"
{
type kqRWallFunction;
value uniform 1.3;
}

ground
{
type zeroGradient;
}

#include "include/sideAndTopPatches"
}



This worked for me and the solution converged.It even worked for a vertical flow like wind direction ( -5 -5 0) Northeast-->SouthWest wind which was causing an immidiate crash for me.

I hope somebody finds this usefull;I was getting annoyed on the posts where the guy explained his situation and when he got the answer he wrote like "Okay I got the solution thanks" and didn't explain the outcome :).

Cheers
Burak