This topic is dedicated to the discussion of BEM 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
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.
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.
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.
Hi,
Thanks for the quick response.
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.
Hi all,
I have two questions:
- 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? - 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
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
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.5rhoUinf^2*dA),
where dA = 2pir*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 = B0.5rhoVtotal^2[ Clcos(phi) + Cdsin(phi) ]cdr.
The total velocity seen by the blade element is:
Vtotal = Uinf*(1-a)/sin(phi).
Also, the local tip speed ratio is:
sigmaprime = Bc /( 2pi*r ).
Substituting all of these equations into the equation for CT above yields Eq. 22:
CT = sigmaprime*(1-a)^2*[ Clcos(phi) + Cdsin(phi) ] / sin(phi)^2
I hope that helps.
Best regards,
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
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: 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,
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
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: nwtc.nrel.gov/AirFoilPrep. There are basic instructions in the “ReadMe” worksheet.
Best regards,
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
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,
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
Dear Guilherme,
I’m sorry, but I’m not sure I understand what your question is enough to answer it. Also, the images you uploaded are a bit too coarse to interpret e.g. I’m not sure which curve is FAST and which is your model. Please clarify your question.
Best regards,
Dear Jason,
I’m sorry for the confusing question. First of all by clarifying the image, the FAST response is in blue, while the response taken from the model is black. Second, in the image it can be seen that the out-of-plane displacement are very close in both simulations (Fast and Model), but the in-plane displacement is very different.
For example, the tip blade in-plane displacement on FAST oscillates around -0.5 whereas in the model the mean is zero. And in the tower, the model’s displacement is very small compared to the FAST’s displacement.
My question is wider and I’m sorry for that, but I checked the model and I did not find something that would make that difference. I would like to know, with your experience, what could cause this difference only in displacement in the plane?
Thank you for your patience and cooperation.
Best Regards
Guilherme Terceiro
Dear Guilherme,
Does your model account for structural pretwist, as FAST does? As discussed in the following forum topic: Coupled blade modes in FAST, because the blade flapwise stiffness is typically quite a bit less than the edgwise stiffness, the influence of the edgwise bending on the flapwise tip deflection (due to pretwist) will be much less than the influence of flapwise bending on the edgewise tip deflection. Perhaps this is why the results between FAST and your model match for out-of-plane deflection, but not for in-plane deflection.
Best regards,
Dear Jason,
Thanks for your answer.
Yes, my model account for structural pretwist. But after reading the topic you recommended, I do not know if I’m doing it right.
In the model I used the second area product of inertia obtained through the principal second area moments to take into account the coupling. I used the mode shape of the first edgewise mode as in-plane mode shape and the first flapwise mode as out-of-plane mode shape.
Can this big difference in results be due to these considerations? The equations that I used are attached.
Again, thank you very much for the help.
Best Regards
Guilherme Terceiro
Dear Guilherme,
Your approach for accounting for the structural pretwist is definitely different than the approach used within FAST. Without going into the details, your approach seems to be missing the influence of the twist on the beam curvatures. Many years ago back in graduate school, I wrote a paper about finite elements of pretwisted beams; please find this paper attached and compare equation (21) to your formulation of the strain energy of the beam. Hopefully a review of this paper will help you resolve the problem.
Best regards,
PretwistedBeamReport.pdf (391 KB)