## Blade-Element/Momentum Theory and Implementation

This forum is dedicated to discussions of the redevelopment of AeroDyn.

Moderators: Bonnie.Jonkman, Jason.Jonkman

Marshall.Buhl
Posts: 437
Joined: Fri Oct 21, 2005 10:22 am
Organization: NREL/NWTC
Location: Boulder, CO
Location: Boulder, CO
Contact:

### Blade-Element/Momentum Theory and Implementation

This topic is dedicated to the discussion of BEM theory and implementation.
Mr. Marshall L. Buhl Jr.
NWTC-3811
National Renewable Energy Laboratory
Golden, CO 80401 USA
Marshall.Buhl@nrel.gov
Voice: +1 (303) 384-6914
Cell: +1 (303) 915-6623
Fax: +1 (303) 384-7079

Andreas.Ferber
Posts: 25
Joined: Wed Nov 11, 2009 2:43 pm

### Re: Blade-Element/Momentum Theory and Implementation

Hi my Name is Andreas
and I am pretty new in the wind turbines community. I am currently writing my Bachelor thesis on the topic of twist/bend coupling.

To get a deeper insight what is going on with the BEMT I try to make a simple model of the 5MW Offshore Baseline Turbine with basic BEMT including tip loss, Glauert correction . I am trying to understand the theory better because I want to use programms like Aerodyn or WT_perf for my later work. I now have some question I can not resolve on my own and I can´t find anything in the literature or papers about it.

1.Question:
First I tried to implement a rectangular plate in the BEMT there is my first problem. When I try to go to r-> 0 my algorithm breaks down. R=0 is impossible due to a singularity in the Math I know. But when I go close to zero (what should be possible?) my tangential induction factor begins to get strange values and in the next convergence step the Prandtl tip loss function gets complex numbers. I know that the tangential induction factor has a singularity where it jumps from plus to minus infinity due to its structure which resembles 1/(x-1). Are there ways to go very far to the root r->0 and avoid these problems?

2.Question:
I saw in the WT_perf that there is wind shear included. How is this possible for the BEMT. When I include wind shear every wing sees the same wind speed distribution.I can not simulate a wind turbine where every wing sees a different distribution. Do you use a UNSTEADY BEMT in the WT_perf code or how do you include wind shear?

3.Question:
When I implement the wing you used in the NREL offshore 5MW baseline wind turbine I get convergence problems. The algorithm converges but the tolerance is to high about 2%.
The problem lays at a specific radial position. I think I get the problems when I read out the Airfoildata at this position. The read out data for the lift and drag coefficient at this position always jumps between two entries and so does the aixial induction factor it jumps between 2 values. Did you have similar problems or do you have a solution for that behaviour.

I would be so happy if you could help me. I first want to understand what is happening in the BEMT before I use programmes which use it.

Thank you very much

Andreas

Marshall.Buhl
Posts: 437
Joined: Fri Oct 21, 2005 10:22 am
Organization: NREL/NWTC
Location: Boulder, CO
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Andreas.Ferber wrote:1.Question:
First I tried to implement a rectangular plate in the BEMT there is my first problem. When I try to go to r-> 0 my algorithm breaks down. R=0 is impossible due to a singularity in the Math I know. But when I go close to zero (what should be possible?) my tangential induction factor begins to get strange values and in the next convergence step the Prandtl tip loss function gets complex numbers. I know that the tangential induction factor has a singularity where it jumps from plus to minus infinity due to its structure which resembles 1/(x-1). Are there ways to go very far to the root r->0 and avoid these problems?

I guess I do not understand why you are trying to predict the performance near the center of the hub. The speed due to rotation is extremely low there.

Andreas.Ferber wrote:2.Question:
I saw in the WT_perf that there is wind shear included. How is this possible for the BEMT. When I include wind shear every wing sees the same wind speed distribution.I can not simulate a wind turbine where every wing sees a different distribution. Do you use a UNSTEADY BEMT in the WT_perf code or how do you include wind shear?

When you include shear in WT_Perf, it divides the rotor into four (or more if you ask for it) sectors. Because of the shear, the wind speed is different in each sector depending on how high the analysis node is in that segment. WT_Perf is totally steady. Everything is in equilibrium. It assumes that an infinite amount of time has passed to achieve this state since the last change in state.

Andreas.Ferber wrote:3.Question:
When I implement the wing you used in the NREL offshore 5MW baseline wind turbine I get convergence problems. The algorithm converges but the tolerance is to high about 2%.
The problem lays at a specific radial position. I think I get the problems when I read out the Airfoildata at this position. The read out data for the lift and drag coefficient at this position always jumps between two entries and so does the aixial induction factor it jumps between 2 values. Did you have similar problems or do you have a solution for that behaviour.

The version of WT_Perf that is available on the web does not converge in some situations. Sometime it does not converge because there may be no solution for a given case and condition. Sometimes it gets hung up on a local minimum/maximum. We are trying to find time to finish an new version with a more-robust iteration algorithm, but we are spread so thin here. I started working on it again last week, but had to drop it due to other priorities. It is at the top of my list of thing to do on Design Codes and I hope to resume work on it next week.
Mr. Marshall L. Buhl Jr.
NWTC-3811
National Renewable Energy Laboratory
Golden, CO 80401 USA
Marshall.Buhl@nrel.gov
Voice: +1 (303) 384-6914
Cell: +1 (303) 915-6623
Fax: +1 (303) 384-7079

Andreas.Ferber
Posts: 25
Joined: Wed Nov 11, 2009 2:43 pm

### Re: Blade-Element/Momentum Theory and Implementation

Hi,
Thanks for the quick response.

[quote="Marshall.Buhl"] Andreas.Ferber wrote:1.Question:
First I tried to implement a rectangular plate in the BEMT there is my first problem. When I try to go to r-> 0 my algorithm breaks down. R=0 is impossible due to a singularity in the Math I know. But when I go close to zero (what should be possible?) my tangential induction factor begins to get strange values and in the next convergence step the Prandtl tip loss function gets complex numbers. I know that the tangential induction factor has a singularity where it jumps from plus to minus infinity due to its structure which resembles 1/(x-1). Are there ways to go very far to the root r->0 and avoid these problems?

I guess I do not understand why you are trying to predict the performance near the center of the hub. The speed due to rotation is extremely low there.[/quote]

I do not want to predict the performance near the root. It was more a first try for me to get in touch with the BEMT. So I wrote a Matlab code and divided the wing in different sections. The first one was very close to the root. But close to the root I could not get values that made any sens to me and mostly the algorithm didn't converge. Then I tried some days to modify my code so that I get reasonable values for the induction factors. But I didn't succeed.
My question was if you have ways in your code to go near the root or if it is impossible. Then I would have to use a cutoff length of let's say 10% of the span. The question was more out of curiosity and interest in the BEMT and not because I wanted to predict the blade performance near the root. Thanks again

Marshall.Buhl
Posts: 437
Joined: Fri Oct 21, 2005 10:22 am
Organization: NREL/NWTC
Location: Boulder, CO
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Andreas.Ferber wrote:I do not want to predict the performance near the root. It was more a first try for me to get in touch with the BEMT. So I wrote a Matlab code and divided the wing in different sections. The first one was very close to the root. But close to the root I could not get values that made any sens to me and mostly the algorithm didn't converge. Then I tried some days to modify my code so that I get reasonable values for the induction factors. But I didn't succeed.
My question was if you have ways in your code to go near the root or if it is impossible. Then I would have to use a cutoff length of let's say 10% of the span. The question was more out of curiosity and interest in the BEMT and not because I wanted to predict the blade performance near the root. Thanks again

To be honest, I've never studied BEM near the root. I don't think I've ever run the code with a non-zero hub radius. I always assumed the impact on the performance that close to the center of rotation would be negligible and the blade always(?) attaches to a hub of non-zero diameter.
Mr. Marshall L. Buhl Jr.
NWTC-3811
National Renewable Energy Laboratory
Golden, CO 80401 USA
Marshall.Buhl@nrel.gov
Voice: +1 (303) 384-6914
Cell: +1 (303) 915-6623
Fax: +1 (303) 384-7079

Marco.Capuzzi
Posts: 23
Joined: Sun May 09, 2010 6:44 am

### Re: Blade-Element/Momentum Theory and Implementation

Hi all,

I have two questions:
1) By looking at the document "AeroDyn Theory Manual", I do not undestand expression 22 at page 9.
From my understanding it seams that the thrust coefficient is equal to the second term present in
the right hand of the expression shown, therefore the one should be removed. Am I wrong?
2) Post processing the results obtained by WT_Perf in a certain condition on a certain blade element, if I
compute the thrust by using Blade Element and Actuator Disc theories I obtain two different values
(quite different). So, if I use the induced speed result, it seams this value does not give equal thrust values
in the two theories, as it should be. I'm doing the check unabling the tangential induced speed functionality
and the hub loss factor.
How is this possible? Are expression 3 and 5 (always from "AeroDyn Theory Manual")
the ones used in the program to compute these two thrusts? (in expression 5 I'm adding the tip loss factor because
I enabled this in WT_Perf)

Regards

Marco

Marco.Capuzzi
Posts: 23
Joined: Sun May 09, 2010 6:44 am

### Re: Blade-Element/Momentum Theory and Implementation

Hi all,

I need to point out my last post. I found out where I was wrong in calculating the thrust by using
blade element theory. I forgot to put the air density! Sorry!
Now everything is fine and I can completely understand and reproduce WT_Perf results.
For what concerns the expression 22 at page 9 of the "AeroDyn Theory Manual", I can still
not understand why there is that number "one" in the right hand.

Regards

Marco

Jason.Jonkman
Posts: 5082
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Dear Marco,

I'm not sure what problem you have with Eq. 22 from page 9 of the AeroDyn Theory Manual. Eq. 22 is simply found as explained below.

The equation for thrust coefficient on an annular section is:

CT = dT/(0.5*rho*Uinf^2*dA),

where dA = 2*pi*r*dr is the differential area of the annular section. The differential thrust in an annular section from Eq. 3 of the AeroDyn Theory Manual is:

dT = B*0.5*rho*Vtotal^2*[ Cl*cos(phi) + Cd*sin(phi) ]*c*dr.

The total velocity seen by the blade element is:

Vtotal = Uinf*(1-a)/sin(phi).

Also, the local tip speed ratio is:

sigmaprime = B*c /( 2*pi*r ).

Substituting all of these equations into the equation for CT above yields Eq. 22:

CT = sigmaprime*(1-a)^2*[ Cl*cos(phi) + Cd*sin(phi) ] / sin(phi)^2

I hope that helps.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Marco.Capuzzi
Posts: 23
Joined: Sun May 09, 2010 6:44 am

### Re: Blade-Element/Momentum Theory and Implementation

Hi Jason,

thanks for the answer. The expression you wrote in the last post is the one I consider correct.
Nevertheless, the expression reported in "AeroDyn Theory Manual" (2005) is different, because
there is a "number one" which is added to the correct expression. If you want you can check it.
I'm writing about this just to know whether it is a typing mistake or I'm wrong.

Regards

Marco

Jason.Jonkman
Posts: 5082
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Dear Marco,

I did some research and I think I found the problem. There was a correction to the AeroDyn Theory Manual, dated December 2005. The corrected version is available from: http://wind.nrel.gov/designcodes/simulators/aerodyn/. The version of the AeroDyn Theory Manual that I found in the NREL publications database is the old version, dated January 2005. In the corrected version (December 2005), Eq. (22) is stated correctly. In the old version it is not, as you rightly pointed out. I'll see if we can update the NREL publications database to have the corrected version uploaded.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Abdulaziz.Abutunis
Posts: 3
Joined: Mon Sep 28, 2015 8:52 pm
Organization: Missouri University of Science and Technology
Location: Missouri, USA

### Re: Blade-Element/Momentum Theory and Implementation

Hi all,

Please I have question regarding the airfoil data file that BEM reads from.
When establishing the cl and cd vs alfa curves, what Reynolds number (Re No.) should we considered. For example, if I am calculating at different flow speeds or tip speed ratios, should I use the incoming flow velocity or the relative velocity (which changes with section location, RPM, and flow velocity) to calculate Re No.

If I am considering stall delay at a section which is function of radius/chord (r/c) and rotational speed (w), that means for each section we need data curve but establishing this curves is time consuming.

Can I calculate Re numbers based on the incoming flow and then cl, cd vs alpha are created at these Re Numbers. Then use Viterna method to extrapolate. But what about the 3D effect on stall delay, if we consider it, then the cl cd vs alpha curves needed to be extrapolated to higher angle of attack for each of these r/c and w.

There should be clear steps of how to do this yet I could not find it in the literature.
Note: I am using xfoil to generate 2D airfoil data at different Re's.

Thanks,
Abdulaziz Abutunis

Jason.Jonkman
Posts: 5082
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Dear Abdulaziz,

When calculating the Reynolds number for airfoil data, you should use the velocity of the wind relative to the blade, not the undisturbed wind-inflow velocity. This relative velocity will depend on the undisturbed wind-inflow velocity, rotational speed, local radius, and induced velocities, but the induced velocities can usually be ignored.

Yes, stall delay depends on local radius, chord, and rotational speed. We typically only use the rated speed, but for a variable-speed rotor, you could assess the effects of different speeds. Stall delay is most important inboard, so, we typically use several airfoil data boards at least across the inboard stations of the blade. Ideally, you would use a separate airfoil table for each blade station, but I understand that this time consuming and the additional accuracy might not warrant the additional time. You'll have to make your own choices between the accuracy you need and the time you want to spend. BEM is only as accurate as the input data used, but without validation using data from an operational rotor, it is not always easy to know how accurate the input data is.

Our AirfoilPrep tool can be used to generate the airfoil data needed by AeroDyn: https://nwtc.nrel.gov/AirFoilPrep. There are basic instructions in the "ReadMe" worksheet.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Dayuan.Ju
Posts: 9
Joined: Wed Mar 06, 2013 10:25 am
Organization: university of Calgary

### Re: Blade-Element/Momentum Theory and Implementation

Hello all,

Did anyone try to consider blade velocity in BEMT implementation by using equation (2) in "AeroDyn Theory Manual"? I always get convergence problem if I use equation (2) for BEMT implementation. If I use equation (1) in "AeroDyn Theory Manual", it works well. I think blade velocity is important for wind turbine aeroelastic analysis and it should be included in inflow angle calculation. If anyone have met this problem, can you explain to me how you solved it? Thanks.

Best regards
Dayuan Ju

Jason.Jonkman
Posts: 5082
Joined: Thu Nov 03, 2005 4:38 pm
Location: Boulder, CO
Contact:

### Re: Blade-Element/Momentum Theory and Implementation

Dear Dayuan,

AeroDyn v12-v14 does use Eq. (2) without problems. In AeroDyn v15, the induction is also applied to the structural motion v_e-op and v_e-ip, also without problems.

Best regards,
Jason Jonkman, Ph.D.
Senior Engineer | National Wind Technology Center (NWTC)

National Renewable Energy Laboratory (NREL)
15013 Denver West Parkway | Golden, CO 80401
+1 (303) 384 – 7026 | Fax: +1 (303) 384 – 6901
nwtc.nrel.gov

Guilherme.Terceiro
Posts: 0
Joined: Thu Nov 10, 2016 6:23 pm
Organization: Universidade Federal do Ceará
Location: Brazil

### Re: Blade-Element/Momentum Theory and Implementation

First of all I would like to thank you for being part of this community. Congratulations to the Team for the work, helping students and professionals in the field.

I am a master's student and I am studying the control of structural vibrations in wind turbines. I'm building a wind turbine model and I use Fast to validate it. But despite all my effort in developing the model it still responds very differently from Fast and to solve this error I need the help of experts in implementing the BEM method.

The same conditions were used to simulate the displacement of the structure in my model and Fast. I used NREL's 5MW reference turbine. The wind was simulated in steady conditions varying exponentially with height with a wind speed of 12m/s in hub height. The BEM method with steady airfoil aerodynamics model was used in the simulation. The structural movement was added the equations of the BEM method, to include aerodynamic damping.

I checked the inflow angle and other parameters of the BEM method, I noticed a significant difference in the values ​​between the model and the Fast, but I can not identify what caused this difference in the displacement in the plane of rotation. The amplitude of the gravitational force is greater than the amplitude of the aerodynamic force, making the total force contrary to the rotation. How do I correct this?

I am attaching pictures to illustrate the difference between the answers.

I am very grateful for your help in this implementation.

Best Regards
Guilherme Terceiro
Attachments
tower tip displacement
tower displacement.jpg (113.85 KiB) Viewed 9857 times
blade 1 tip displacement
blade1 displacement.jpg (131.81 KiB) Viewed 9857 times