airfoil in OpenLB

[balzar29]

Hi,
how can i costruct a 2D airfoil profile instead of cylinder in the cylinder 2D example to run a simulantion over a 2D wing profile?

I mean, for me the example is fine, i just need to remove the cylinder and add a profile and maby extend the domain.

Thx for answers!

[robin.trunk]

Hi balzar29,

right now there is no functor for that. But if your airfoil can be described by a function like here
https://en.wikipedia.org/wiki/NACA_airfoil
you can create your own functor. It should be mainly copy-paste work, just have a look at
src/functors/indicator/indicatorF2D.cpp
src/functors/indicator/indicatorF2D.h
src/functors/indicator/indicatorF2D.hh
You can copy the IndicatorCircle2D, rename it and adapt the equation. Since there are probably more people interested in such simulations it would be great if would consider to send in your code for the airfoil-indicator and thus contribute to future releases of OpenLB.

Best
Robin

[balzar29]

i’ll try! If i’ll be able to make it work i’ll post something!
Thanks for the help!

[balzar29]

I’m sorry, but my c++ is very poor. I don’t understand where to change the function. I think that i’ve found the place where insert the airfoil equation, but if someone can halp me it would be better and faster. I need just some c++ explanation of this code in the file “indicatorF2D.hh”

Code:

template <typename S>
IndicatorCircle2D<S>::IndicatorCircle2D(Vector<S,2> center, S radius)
: _center(center),
_radius2(radius*radius)
{
this->_myMin = _center – radius;
this->_myMax = _center + radius;
}

// returns true if x is inside the circle
template <typename S>
bool IndicatorCircle2D<S>::operator()(bool output[], const S input[])
{
output[0] = ( std::pow(_center[0] – input[0],2) + std::pow(_center[1] – input[1], 2) <= _radius2 );
return output[0];
}

template <typename S>
bool IndicatorCircle2D<S>::distance(S& distance, const Vector<S,2>& origin, const Vector<S,2>& direction, int iC)
{
S a = direction[0]*direction[0] + direction[1]*direction[1];

// returns 0 if point is at the boundary of the sphere
if ( a == _radius2 ) {
distance = S(0);
return true;
}
// norm of direction
a = sqrt(a);

S b = 2.*((origin[0] – _center[0])*direction[0] +
(origin[1] – _center[1])*direction[1])/a;
S c = -_radius2 + (origin[0] – _center[0])*(origin[0] – _center[0])
+ (origin[1] – _center[1])*(origin[1] – _center[1]);

// discriminant
S d = b*b – 4.*c;
if (d < 0) {
return false;
}

S x1 = (- b + sqrt(d)) *0.5;
S x2 = (- b – sqrt(d)) *0.5;

// case if origin is inside the sphere
if ((x1<0.) || (x2<0.)) {
if (x1>0.) {
distance = x1;
return true;
}
if (x2>0.) {
distance = x2;
return true;
}
}
// case if origin is ouside the sphere
else {
distance = std::min(x1,x2);
return true;
}

return false;
}

template <typename S>
bool IndicatorCircle2D<S>::normal(Vector<S,2>& normal, const Vector<S,2>& origin, const Vector<S,2>& direction, int iC)
{
S dist;
if (!(distance(dist, origin, direction, iC)) ) {
return false;
}

Vector<S,2> intresection(origin + dist*direction);

normal = intresection – _center;

return true;
}

Thanks!!

[robin.trunk]

Hi balzar29,

the Indicator is a functor which is called at another point in the library, when the material numbers are set.
In the
IndicatorCircle2D<S>::IndicatorCircle2D
the functor is initialized, here is the place to hand in parameters required for the desired equation.
The indicator return for coordinates, whether this coordinates are inside or outside the object. This is done in
bool IndicatorCircle2D<S>::operator()
This is the place to implement the equation describing the airfoil. “input” is the physical coordinate (input[0] for x-coordinate and input[1] for y-coordinate) which needs to be checked. “output[0]” is a boolean value that is returned. “true” or “1” if the point is inside the airfoil and “false” or “0” if the point is outside the airfoil.

Best
Robin

[balzar29]

Hi robin,

thanks to your help now i can change shape in 2d cylinder example. This is the result for thin wing with no aerodynamic trailing or leading edge, just the shape!

Now my problems are in the boundery conditions. When i start the computation it gives me these errors after defining start parameters:

Code:

[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (0,29) ~ (1.04833,0.751667), in direction 8
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,29) ~ (1.05,0.751667), in direction 1
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,33) ~ (1.05,0.758333), in direction 3
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (0,33) ~ (1.04833,0.758333), in direction 4
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (0,29) ~ (1.04833,0.751667), in direction 7
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,29) ~ (1.05,0.751667), in direction 8
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,29) ~ (1.05167,0.751667), in direction 1
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (0,33) ~ (1.04833,0.758333), in direction 5
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,33) ~ (1.05167,0.758333), in direction 3
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,33) ~ (1.05,0.758333), in direction 4
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,29) ~ (1.05,0.751667), in direction 7
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,29) ~ (1.05167,0.751667), in direction 8
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (3,29) ~ (1.05333,0.751667), in direction 1
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (1,33) ~ (1.05,0.758333), in direction 5
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (3,33) ~ (1.05333,0.758333), in direction 3
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,33) ~ (1.05167,0.758333), in direction 4
[BlockGeometryStatistics2D:6] updated
[sOffLatticeBoundaryCondition2D:6] Cuboid globiC 6 finished reading distances for ZeroVelocity Boundary.
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,29) ~ (1.05167,0.751667), in direction 7
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (3,29) ~ (1.05333,0.751667), in direction 8
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (4,29) ~ (1.055,0.751667), in direction 1
[BoundaryConditionInstantiator2D:3] ERROR: no boundary found at (2,33) ~ (1.05167,0.758333), in direction 5
[BlockGeometryStatistics2D:2] updated
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (240,427) ~ (1.44833,0.708333), in direction 6
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (240,426) ~ (1.44833,0.706667), in direction 7
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (241,426) ~ (1.45,0.706667), in direction 8
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (242,426) ~ (1.45167,0.706667), in direction 1
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (241,426) ~ (1.45,0.706667), in direction 7
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (242,426) ~ (1.45167,0.706667), in direction 8
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (243,426) ~ (1.45333,0.706667), in direction 1
[BoundaryConditionInstantiator2D:2] ERROR: no boundary found at (242,426) ~ (1.45167,0.706667), in direction 7

and so on, then the computation starts and gives results.

Have you some ideas? Thank!

[balzar29]

PS: this is for non symmetrical airfoil, for the symmetrical ones it’s easiest but the boundery condition problem remains.

[mathias]

You should use Bounce Back as boundary condition.

Best
Mathias

[balzar29]

Last questions.
I’ve removed the Poiseuille inflow profile. I need to set periodic boundery condition on top and bottom and a constant inlet velocity and a gauge pressure in outlet. How can i have control on these parameters?
Now with the defined boundery condition when i increase Reynolds number over 200 the solution problem diverge…i can’t find where to change velocity, density and pressure.

[robin.trunk]

Hi balzar29,

for the periodicity:
you can have a look at the thermal2D example, more specific the line
cGeometry.setPeriodicity( true, false );
this sets periodicity for x, but not y direction.

For the velocity:
You can hand in an Analytical functor instead of a Poiseuille profile
std::vector<T> zero( 2,T() );
AnalyticalConst2D<T,T> u( zero );

To control the parameters you want to have a look at the LBconverter at the beginning of main
http://optilb.com/DoxyGen/html/d1/d25/classolb_1_1LBconverter.html

Best
Robin

[balzar29]

Hi,
I’ve set the periodicity boundery layer on top and bottom, than the Poiseuille velocity profile in inlet is removed, but i still have stability problems at high reynolds number (2e6) and also at small reynolds number with a symmetric profile with no angle of attack (it should have only little drag and no lift) at a certain time the profile creates lift and on the top of it separation appears as if there is an angle of attack.

You know why the program runs this way? Thank you.

[robin.trunk]

Hi Balzar,

for high Reynolds numbers you need to apply a turbulence model (e.g. SmagrosinskyBGKdynamics). To do this you have to change the the dynamics, you can have a look at the nozzle3d example, I think about line 349.
For small Reynolds numbers, what do you mean with “symmetric profile”? As far as i understand you are using a contant inflow velocity over the whole inlet? In this case you could check the corners, since there could occur some errors due to the non-smooth transition. Another problem could be the obstacle itself. If there is almost no angle of attac the obstacle maybe is just one cell thick? Such thin spikes can cause numerical errors, maybe increasing the resolution of the grid helps at this point.

Best
Robin

[balzar29]

Hi robin,

i’m back on this post.
I found out some instability problems. I’vs set Smagorinsky model to investigate high Reynolds number, then the boundery condition are set as inlet velocity (left one), pressure outlet (right one) and periodic (top and bottom). If i delate the Poiseuille inlet profile, the problem crash at the top and bottom left corner where periodic boundery meets the inlet. Leaving the Poiseuille gives more stability, but after a certain time the instability strikes out and i really don’t know how to solve this.
I’ve set latticeL=x and Delta_T=x^2 to avoid the compressibility problem. What i don’t understand is why latticeU is defined in the guide like delta_t/delta_x (page 26). I know that delta parameters are adimensional, but velocity is space frac time not time frac space. Anyway like delta_x=latticeL/charL, maby also delta_u=latticeU/charU then if charU=1 the result become delta_u=latticeU, but the problem is still that delta_u=delta_t/delta_x (time/space) and not delta_x/delta_t (space/time).

Just using the guide, if i want delta_t=delta_x^2 it should be latticeU=latticeL because latticeL=delta_x=x, latticeU=delta_u=delta_t/delta_x, but delta_t=delta_x^2=x^2 so delta_u=x^2/x –> delta_u=x that is latticeU=x.

Am I doing something wrong?

Thanks for help!

[balzar29]

*i’ve set walls on top and bottom with Poiseuille inlet not periodic boundery

[balzar29]

PS: how can i prevent the backflow when the perturbation arrives to the outlet?

Sorry for my insistence.

[Author]

Posts