Comparing to QBlade results: No tip loss correction

[hcOnel]

THIS ISSUE IS OBSOLETE. THE CASE WILL BE KEPT FOR POSSIBLE FUTURE NEEDS.

I have written an Octave script which compares turbinesFoam results to QBlade results.

Below is the case details:

• Axial flow, single turbine
• Tip loss correction: None
• Inflow velocity = 11.4 m/s
• TSR = 7
• Multi-Re solution is used

Note that forces are not per density, they are multiplied by rho. Solution has converged.

Local velocity and Reynolds are in very good agreement, however Cl and AoA values are deviant in mid-span region. Also, there is a huge difference in local Axial and Normal force distribution values.

Possible causes: cl cd data are incorrect?

[petebachant]

My guess is that turbinesFoam is computing the forces differently, maybe in a different coordinate system. If I remember correctly, the tangential direction is along the mean chord line of each element, so to put into the rotor coordinate system, the pitch and twist would need to be taken into account. I would need to look at the actuatorLineElement::writePerf code to confirm.

[petebachant]

Can you add the post-processing script to this repository?

[hcOnel]

I have added all pre and post processing scripts under a single directory. The only post-proc script is postProc.m. Please ask anything that is unclear, I will reply ASAP. Scripts are Octave.

[hcOnel]

About your note on calculation of forces in turbinesFoam, as far as I know, normal force component (tangential in QBlade) is used for power and axial (normal in QBlade) is used for thrust calculation in turbinesFoam, so they should be in the rotor plane coordinate system (i.e, normal force is on the rotor plane and perpendicular to blade span, axial force is perpendicular to rotor plane). Your clarification would be really helpful.

[petebachant]

The forces written from each element are written in the global fixed coordinate system. Some additional code would be needed to write them in the rotor coordinate system, or they can be transformed after the fact since the element's global position coordinates are also written.

[hcOnel]

If you are refering to the fx fy fz in the postProcessing elements files, the script converts them into the blade (rotor) local coordinate system by applying a dot product with the unit normal vector of the blade (coincident with the rotor plane and pointing in the rotation direction) at its instantaneous position. So yes, this assumes fx fy fz to be in global CS in the first place and converts it.

[hcOnel]

Element local velocity and Reynolds numbers are in excellent agreement with QBlade, which indicates that chord distribution and local velocity calculations are OK. However, parameters like Cl, Cd, axial and normal forces aside; local angle of attack has a very large deviation, changing its error direction at r~13m. My guess: it may be related to positioning of chord mount point, which affects the angle of attack calculation.

I will now try increasing the number of elements to see the result, because since Lift and Drag equations include the span length and span length is large per element when the number of elements is low, this might also have an effect.

Another question: in turbinesFoam, what is the difference of AoA and AoA_Geom? Since AoA_Geom remains constant, I guess it is a theoretical value which is calculated from BEM method? Also, are forces calculated using AoA or AoA_Geom?

[petebachant]

Forces are calculated using the true angle of attack. The geometric angle of attack is calculated without induction, i.e., assuming the inflow is the free stream velocity. Angle of attack could be different from QBlade due to deflection of the flow by the hub. I don't know if QBlade (AeroDyn?) includes that effect.

[hcOnel]

Results of QBlade here are also from a simple solution of BEMT (Blade Element Momentum Theory) as far as I know, which means that AoA_Geom and QBlade's AoA are both free of induction (assuming uniform freestream, no hub influence), so at least they should be in a very close agreement (red dots and blue line) but they are not.

I should also note that this QBlade solution of NREL5MW is a test case readily comes with QBlade and its results are well-matched with NREL's own results as far as I know, that is why I have chosen it for benchmarking. I am also constantly re-checking my understanding for turbinesFoam and QBlade for possible mistakes.

My primary suspicions are inconsistent positioning of airfoil leading edges (for AoA calculation) and/or an error in cl cd data. Because relative velocity and Reynolds values match perfectly.

Also, if you check the postProc.m script, I calculate elementary Fa and Fn values (using data of turbinesFoam only) by calculationg Lift and Drag forces manually (0.5 rho V^2 span chord * coeff) and plotting it on top of Fa Fn data available in postProcessing elements files, which also has some inconsistency in Fn curve. I am also trying to figure that out.

My ultimate goal here is to find out why I get a very large difference in power coefficient for tip loss on/off cases in turbinesFoam, and Cp passes the Betz limit for tip loss off case even though the local velocity velocity magnitue is consistent.

[petebachant]

In BEMT the induction is computed from momentum theory, so its reported angle of attack should be compared with the true, not geometric angle of attack.

One thing I noticed is that the QBlade force results seem to be in units of N/m, which would explain why the turbinesFoam results show larger discrepancies at larger span values.

[hcOnel]

Of course... So I should divide the reported Fa and Fn of turbinesFoam by span length. Doing it right away, thanks sir.

[hcOnel]

The results become: Now I need to find out why Cp is 0.71 in turbinesFoam and 0.48 in QBlade while the force distributions along the blade are almost identical (although there are some elements where the forces per length are slightly higher).

[petebachant]

That I'm not sure of. It could be a bug in how the moments are calculated/normalized. If you have a chance, you can peer through the turbinesFoam source code and double check for errors.

[hcOnel]

Actually, after zooming in to the normal force curve, difference is significant: I guess the problem may boil down to AoA calculation after all.

[hcOnel]

In the script, I have calculated power coeff (cp) and thrust coeff (ct) manually and compared them to turbinesFoam results. (Needs verification: in turbinesFoam, ct and cd are torque and drag coefficients respectively, and cd corresponds to thrust coefficient in that sense.)

Power: multiply normal force per element span by element span by radial position for each element and sum them up. Thrust: multiply axial force per element span by element span for each element and sum them up. Since I have results for one blade only, I multiplied power and thrust calculated for one blade by number of blades to find turbine power and turbine thrust respectively. I know this would cause an error but I expect it not to be significant in a sense of seeing the overall behavior. Coefficients: Divide turbine power and turbine thrust by wind power and wind force respectively.

Below is the result I have obtained:

The valleys are caused by the presence of tower, but overall, curves have a large deviation. I have explained my calculations above deliberately to see if I have a different method than turbinesFoam.

[petebachant]

That all seems consistent, yet there is a significant discrepancy. This definitely warrants double checking the turbinesFoam code.

[hcOnel]

UPDATE: Those latest results are incorrect, since I have forgotten to multiply force components by density (for the sake of comparison with QBlade) whereas wind power calculation included density; so there were a mistake. Manually calculated cp is closer turbinesFoam's calculated cp (but still not equal). Basically, axial and normal force graphs are scaled by the density value (1.225 in that case).

However, it bugs me that before this "error" of mine, axial force matches perfectly, normal force graph is closer and power coefficient is really in good arrangement with QBlade and some publications. Might there be a bug involving density specification of power and forces? p.s: I have used to multiply with rho as I said in my first post but it somehow went away during my editings I guess.

[petebachant]

It looks like the AoA calculation is fine, but the Cl values are significantly higher. Are you using the same foil data for both?

[hcOnel]

Yes. The polar data I am using in turbinesFoam is generated by QBlade's Multi-threaded batch analysis tool. I have generated them directly from the QBlade's NREL_5MW.wpa project; so same airfoils. I have used free transition (top and bottom forced transition location set to 1). Overall, the same polar is used for both. You said AoA calculation is fine, but if you zoom into 0-10 degrees range, you will see that there is a 2-4 degrees error in midspan region. And cl cd curves are very sensitive to small AoA changes, a 1 degree change in AoA may result in a significant cl change. Also, any ideas on the density multiplication I have mentioned? Did you have the time to check the code? Best regards.

[petebachant]

Did you do a grid size dependence check? I'm wondering if the proper amount of force is not being imparted on the flow field.

[hcOnel]

Although I did not note the exact grid sizes and results, yes, I did do informally. I obtained solutions for 3 different mesh resolutions, starting from the coarsest and doubling the number of cells in each directions at each step. My measurement was the final c_p and c_t values, and they remained the same. The axial flow tutorial provided in the original code also returns high c_p values without tip loss correction. Can the amount of force imparted on the flow field has that much effect on turbine power anyways?

[petebachant]

If the force it too low, it won't slow the flow enough, which would lower the local speed ratio and therefore raise the angle of attack.

Following SOWFA, turbinesFoam sets the momentum source as the force per unit density (tapered by the Gaussian function), but OpenFOAM's actuator disk appears to add force in the cell proportional to the fraction of the total disk volume the cell takes up. This makes sense since it's computing the total thrust from the entire disk rather than from a blade section, so that may not be relevant.

template<class RhoFieldType>
void Foam::fv::actuationDiskSource::addActuationDiskAxialInertialResistance
(
    vectorField& Usource,
    const labelList& cells,
    const scalarField& Vcells,
    const RhoFieldType& rho,
    const vectorField& U
) const
{
    scalar a = 1.0 - Cp_/Ct_;
    vector uniDiskDir = diskDir_/mag(diskDir_);
    tensor E(Zero);
    E.xx() = uniDiskDir.x();
    E.yy() = uniDiskDir.y();
    E.zz() = uniDiskDir.z();

    vector upU = vector(VGREAT, VGREAT, VGREAT);
    scalar upRho = VGREAT;
    if (upstreamCellId_ != -1)
    {
        upU =  U[upstreamCellId_];
        upRho = rho[upstreamCellId_];
    }
    reduce(upU, minOp<vector>());
    reduce(upRho, minOp<scalar>());

    scalar T = 2.0*upRho*diskArea_*mag(upU)*a*(1 - a);

    forAll(cells, i)
    {
        Usource[cells[i]] += ((Vcells[cells[i]]/V())*T*E) & upU;
    }
}

It may be worth comparing with the SOWFA NREL 5 MW example case. I wonder if that gets Cp and/or force distributions correct.

[hcOnel]

Can you suggest;

• Number of cells along the blade (rotor radius) for mesh size,
• Rotation per time step (i.e. radians/timeStep) for time step size?

[petebachant]

I'd say a good start would be about 200 time steps per rev, and 32 cells per meter in the rotor region.

[hcOnel]

Ok, I have reduced my time step from 0.2s to 0.05s (which means rotor revolution per time step has dropped from 14.4deg/dt to 3.6deg/dt, or time steps per revolution has increased from 25 to 100) before reading your post. Result: cp dropped from 0.54 to 0.5075. There is definitely a convergence problem I am having, and I'll try 200 time steps per rev and observe the convergence. About the mesh size: I used to think of it in a cell per rotor radius basis but I will also try per meter approach (which will dramatically increase my number of cells per diameter value from 50 to 4032, since rotor diameter is 126m).

[petebachant]

From that perspective I ran my simulations around 32 cells per rotor diameter, since it was a one meter rotor. It sounds like you may have a decent spatial resolution, but the timesteps may be too large--the wakes from the blades don't spend enough time in a certain area to affect the flow (and cause the induction).

[hcOnel]

I guess, since the tip is the node which moves the largest distance per time step, size of time step should depend on this value. Say, the tip must travel at most 0.5m per time step (I am making it up for the sake of argument), and then a number of cells per this distance (0.5m) can be defined that would allow the resolution of flow within this distance between time steps. This makes me think that an unstructured mesh may come in handy within the rotor plane, which provides some homogeneity of "cells per meter in the angular direction at each radial position of blade" in a sense of computational resources. Does that sound logical? Thanks for your feedback.

[petebachant]

That does sound logical, but may be a bit overkill. If you decide to give it a shot, let me know.