Hi guys,
I implemented a simple subroutine in order to estimate the induced flow distribution on the rotor. I considered 9 elements for the blade of a three-bladed rotor (44m of blade lenght), starting from 8m of radial position (therefore, every element is 4 meters wide spanwise). I’m computing the normal component of the induced flow (tangential induced flow neglected). The subroutine consists of an iterative cycle which finds out the value of inuced flow on every blade element. The iteration is performed in order to find the normal induced flow (on an annular area of the rotor) which equalizes the thrust values obtained by actuator disc model and blade element model. The problem is that for some wind speeds, just above cut-in, the software does not converge on several elements of the blade. From my checks it seams there is not a solution at all. So, it seams that for this conditions, with fixed wind speed and speed of rotation, there is not a value of induced flow which equalizes the two thrust expressions obtained by the two aerodynamic models. May someone let me know whether this makes sense or I’m doing something wrong? I want to remark that I’m considering a simple static analysis, with a constant speed of wind distribution along the rotor and without any engineering correction.
regards
Marco
Marco,
Depending on your lift curve, you can get 0, 1, 2, or more solutions, so I am not surprised you are not converging in some cases. This is generally not discussed in textbooks on rotor aerodynamics.
I noticed this with WT_Perf a few years ago and we’ve been working on modifications to make the iteration solution more robust.
I’d like to recommend an exercise for you to do. Usually, we make an initial guess for the induction, run it through the equations to get a new guess for the induction. We use the new value as the next initial guess. Instead of doing that, put a loop around the iteration algorithm where you vary the initial guess from -1 to 2 in steps of, say, 0.1. At the end of each iteration and right before the END DO, print out the initial guess and the difference between the initial guess and the iterated value. That is, print Aold and (Anew-Aold). Plot it. This is your solution space and wherever the curve crosses the x-axis is a solution. If you curve looks something like a quadratic that never crosses the x-axis, you have no solution. Do this for a variety of conditions–especially the ones you are having trouble with.
BEM theory is developed with the assumption that the induction is constant across the rotor. Because multiple solutions are possible, your program really should find all solutions for each annulus and pick a combination that gives the lowest deviation across the rotor. We are debating doing that for WT_Perf, but it will really slow down the code. As far as I know, no one does that.
You may want to compare your program to WT_Perf. The version that is currently available does not have the new, more-robust algorithm, but you may find it informative to compare the two programs.
Dear Marshall,
thanks a lot for the answer. It was very helpful.
Before your reply I had already verified that the non-converging cases are not related to numeric issues of the code.
I did this simple check: I plotted, for different values of dimensionless induced speed, the AD (actuator disc) and BE (blade element) thrust in corrispondence of fixed wind speed, rotational speed and radial posiotion.
In the non-converging cases the two curves do not intersect at all. In few cases they intersect, but only for dimensionless induced speeds higher than 0.5. (It is my understanding that this result does not correspont to a realistic condition).
If you do not mind I would like to ask you the following questions:
-Does it make any sense to take a dimensionless induced speed higher than 0.5?
-What value should I consider when the two curves do not intersect? (I need these values in order to compute the complete power curve)
-What is the physical meaning of such behaviour?
In the end I want to remark that the paremeters domain of non-convergence becomes wider by adding a tip loss factor. Indeed, in this way the BE thrust is reduced further on.
looking forward to hearing from you,
best regards
Marco
Marco,
I apply a modified version of the Glauert Empirical Correction to cases when the induction exceeds 0.4. I wrote a short paper on a simple way to tweak Glauert’s correction to work with tip losses. It is purely a numerical trick–it is not based on physics. You can see it here:
[url]http://www.nrel.gov/docs/fy05osti/36834.pdf[/url]
I’m not sure what to use. Are all the cases when they don’t intersect ones where it looks like induction is > 0.5?
I think the physical meaning is that BEM really does not apply for that condition, but I am not really an aerodynamicist–I just play one on TV. 
I can try to talk one of ours into providing better insight if you like.
Marshall
Dear Marshall,
thanks for the answer. I want to specify that in some cases the two curves intersect for values higher than 0.5. On the contrary, in some others they do not intersect at all. Indeed, in these cases the BE thrust is always above (greater) the AD one. Typically, these non-converging cases correspond to low wind speeds (between 4 and 8 m/s) and high radius position (the turbine radius is 44 meters and the convergence problems start nearby 38m). If you can ask to someone inside NREL I would really appreciate it. I want just to know which is the best value that I can take for the induced flow for such non-converging conditions.
I will look at the document you attached to see if I can use also this trick.
Thanks again for your help.
Regards
Marco
Dear Marshall,
I read the document you attached. Thanks a lot, you made me aware of the turbolent windmill state existence, because I wasn’t. So, resuming I can use the Glauert correction and modify it (as shown in the document) when tip loss is included. In this way I should not obtain conditions where the BE thrust is above the AD one for all the “a” values between zero and one. If you don’t mind, in order to refine my modelling I want to ask you something more:
-Are there in literature values of the induced flow distribution at the different wind speeds along a large WT blade. These would allow me to check if my induced flow estimation is reliable or not.
-Do you think that adding also a tangential component of the induced flow I can strongly improve my power curve estimation accuracy. Or it should have a small effect?
-For very small wind speeds ( at cut-in and just a little bit above it) I obtain a negative induced flow at the most inner section that I have considered (at 8m of radial position on a total lenght of 44m); does this result seam wrong to you?
Looking forward to hearing from you,
best regards
Marco
Marco,
Please keep in mind that I am not really an aerodynamicist, but I have learned a few things in my years of playing with WT_Perf. I did ask our aerodynamicists to respond, but it looks like they have not.
I know of nothing in literature that tabulates induction along a large blade. I recommend you compare your results to those from other codes such as Bladed, PROP-ID, or WT_Perf if you want to do a sanity check.
There was debate in the old days about adding tangential induction to the codes. Many thought it wouldn’t make significant difference. We added it to our codes and found that it sometimes did matter. I can’t remember if it was at low or high wind speeds. If you test your code against WT_Perf, you can try WT_Perf with and without it to see how much difference it makes.
At very low wind speeds, the rotor can enter the vortex ring state, where the wind flows around the rotor and then turns back upwind to flow through the rotor. Eggleston’s and Stoddard’s Wind Turbine Engineering Design book has a classic picture of the various rotor states on pages 32 and 33. Or, you can see a similar image on page 41 of Lock’s classic paper here:
[url]National Wind Technology Center's Information Portal | Wind Research | NREL
Have fun!
Marshall
Dear Marshall,
Thanks a lot, I will use WT_Perf for all the checks I need.
I hope that some aerodynamicist could join our discussion about the physical sense of dimensionless induced flow higher than 0.5.
Best regards
Marco
Marshall and Marco,
Reading through your posts, I agree with everything you have said. The state with the induction factor is between 0.5 and 1 is known as the turbulent wake state, where so much momentum is extracted by the rotor that there is some reversed flow in the wake and perhaps some through the rotor - Figs. 4d in Lock’s paper. As you increase the thrust further at induction factor = 1.0, the vortex ring state is reached - Lock’s Figs. 4e.
As far as having negative induced velocities - I have seen this in some calculations where this happens near the root and near the tip, which physically makes sense to me near the tip in the turbulent/vortex wake state, but less so near the root. I’ve seen this is somewhat dependent on grid resolution, so maybe try to add a few more points near the root and see if the calculation still go negative.
cheers,
Pat
Talk of the vortex-ring state brings me back to school days when helicopter flight theory was being discussed. Opened up old Shapiro’s book to refresh my memory and though the bias of the subject matter is for helicopter rotors, there is a lot of discussion in practical terms of the transitions to the other states. Perhaps this would be a helpful line of comparison in your investigation. For instance a quick check that “1/f < 4” would quickly alert you when approaching the turbulent state. I should put this into my own models. Maybe convergence problems that I’ve been having are the same thing, and I didn’t realize it!
PS it’s Principles of Helicopter Engineering, by Jacob Shapiro, 1956, McGraw-Hill, which I find much more readable than Wayne Johnson’s book, in case you are interested.
Hi guys,
thanks a lot to everyone for joining the discussion. I have a new question about the induced flow estimation in BEM.
I Implemented a subroutine for the axial induced flow estimation, including Buhl correction and tip and hub loss factors.
Consideirng many operative conditions (speeds of wind and rotation fixed in every of these), I obtain a dimensionless axial induced flow which is
always (10-20)% higher than the one obtained by WT_Perf. I’m performing the comparison using the same inputs for the aerodynamic
coefficients and the same blade discretization. I’m trying to understand whether this difference is due to a mistake present in the code
or simply to the aerodynamics which is modelled in WT_Perf and is not in mine code (tangential induced flow or other).
Do you have any suggestions?
Because looking at the tangentail induced flow values in WT_Perf output, they seam too little in order to justify this big percentage
different in results. Usually their dimensionless value is lower than 1%.
Is there something else which makes WT_Perf aerodynamics modelling more refined than mine?
For example in my software I have costant aero coefficients with the Reynolds number and in WT_Perf
I input these coefficients also in correspondence of only one Reynolds.
Does the code consider by itself a change of the coefficient with such parameter?
Regards
Marco
Hi,
sorry just one correction to my previuos post: the induced flow values I obtain are
10-20% lower (and not higher) than WT_Perf ones.
Marco
Marco,
I think you answered your own question. Yes, aerodynamic coefficients, particularly airfoil drag, vary with Reynold’s number. Also consider the effect of Re on the slope of the lift-curve. This will have an effect on the induced flow velocity.
Marco,
WT_Perf has a lot of switches for enabling options for the induction algorithm. I would recommend that you review the parameters in the “Algorithm Configuration” section of the input file and try to build a WT_Perf model that uses the same algorithm you are using with your code before you try to compare the codes.
Marshall
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,