OpenFAST / openfast

Main repository for the NREL-supported OpenFAST whole-turbine and FAST.Farm wind farm simulation codes.
http://openfast.readthedocs.io
Apache License 2.0
663 stars 452 forks source link

AeroDyn 15 mesh refinement #1344

Open andrew-platt opened 1 year ago

andrew-platt commented 1 year ago

Feature overview For coupling AeroDyn to CFD codes, it would be really useful to have a method built directly into AeroDyn 15 to refine the aero mesh. In issue #1232, there was a mismatch in position between the aero position mesh and the force mesh passed to CFD codes such as NALU-Wind, AMR-Wind, and SOWFA. This issue had come about as a result of including the structural mesh in the mesh interpolation within the OpenFOAM module (OpFM). While the structural mesh was removed with PR #1324, the interpolation of the aero mesh to a finer mesh that the CFD requires is still handled within the OpFM module. This refinement might be giving rise to discontinuities in the aero loads passed back to the CFD as shown in the following figure.

Figure: Comparison of CFD couplings with blade resolved (BR) CFD. Note the discontinuities in the normal force in the BEM, Vortex-Method, and ALM curves . These discontinuities are likely due to a combination of changes between airfoils in AD15 and mesh mapping from the course aero mesh to the finer CFD mesh near the airfoil change. MicrosoftTeams-image (8)

In addition to possible discontinuities from the mesh mapping, it would be very convenient to include a simple mesh refinement in AD15 directly. This would greatly simplify using OpenFAST with a CFD code.

Describe the solution you'd like We would like to include a mesh refinement parameter in the AD15 input file. This would trigger AD15 to interpolate additional nodes along the blade span, and would require some interpolation of airfoils.

From @mchurchf (comment https://github.com/OpenFAST/openfast/issues/1232#issuecomment-1317119119 from #1232)

Tony and I just had a conversation about interpolation of Cl/Cd polars between airfoils. There is a method implemented in this old NREL tool AirfoilPrep that does a Cl/Cd polar interpolation that makes most sense that you may want to consider. I use it when I prepare actuator line or Aerodyn input files to have more resolution than whatever I am given to start with. I will try to find a reference, but the main idea is that it transforms the alpha in these tables to what I'd call alpha'. In this transformed space, peak Cl occurs at the same alpha' no matter the airfoil (so, for each airfoil, you need to know alpha for Cl,max), then it does linear interpolation based on thickness of airfoil. Once the interpolation is complete, the inverse alpha transform is applied and all is done. The effect is that as you blend between two airfoils with different max Cl alphas, the peak in Cl will smoothly move left or right. If you don't do this, you get double Cl peaks for blended polars, which is completely undesirable. I'm happy to explain this process in more detail. It is a feature we've really wanted for a long time directly in OpenFAST, instead a preprocessing step.

To summarize, AD15 would get the following changes:

In addition, the following changes would be made to the OpFM module:

And some API changes for communication with CFD codes:

Describe alternatives you've considered The current best practices method of refining the AD15 mesh is to preprocess the AD15 blade file and polars to create a more refined blade file with modified polars. This method still leaves some room for error as the CFD mesh is linearly spaced along the blade span or tower, but the AD15 meshes are not uniformly spaced. As show in the figure above, this may be leading to artifacts.

Those who may want to weigh in on this proposal @tonyinme @mchurchf @michaelasprague @gantech @jjonkman @ebranlard

ebranlard commented 1 year ago

This feature would indeed be great.

As a note, on top of the careful interpolation of polars (to avoid double stall), we should implement 3D polar corrections into AeroDyn, because 3D polar corrections typically depend on the radial location (c/r).

The Snel and Du-Selig implementation can be found here together with a fairly standard blending procedure (to blend the corrected polar with the uncorrected polar). This will require:

It's likely fine to first implement the interpolation of polars, and then implement the 3D correction. The error introduced by omitting the 3D correction are likely small.

jjonkman commented 1 year ago

I agree that this feature would be great to have. In addition to the comments above, I would add:

Best regards,

jjonkman commented 1 year ago

This issue is related to https://github.com/OpenFAST/openfast/issues/374.

ietqlw commented 1 year ago

Dear all,

Subroutine “ReadBladeMeshFileAD” in ElastoDyn.f90 reads in the AeroDyn v14.00.00 input file to get the blade discretization used in the structural dynamics module. When employing AD15, how the program get the blade discretization?I do not find the codes that get the blade discretization from AD15 input file.

jjonkman commented 1 year ago

Dear @ietqlw,

When AeroDyn v15 is enabled, the aerodynamic blade discretization used by AeroDyn is defined in the AeroDyn blade input file and the structural blade discretization used by ElastoDyn is defined by ElastoDyn input BldNodes. (That is, the aerodynamic and structural discretizations are independent.) Basically, when AeroDyn v15 is enabled, the ElastoDyn blade is discretized into BldNodes equal-length elements with the blade structural analysis nodes located at the centers of these elements.

Best regards,

ietqlw commented 1 year ago

Dear @jjonkman,

Thank you for your kind answer.

I set _AFTabMod_as 3 in AD15 primary input file for enabling 2D interpolation on AoA and UserProp. The following AD15 blade definition input file was selected from the reference NREL-5MW wind turbine. Eight airfoil-data tables (BlAFID=1~8) were originally used. Now I choose a element to apply UserProp and set BlAFID as 9.

The following doubt confuses me a lot: (1) The table defines 19 blade span properties. what about the properties between two blade nodes? If BlSpn = 9m, which number is the corresponding BlAFID, 2 or 3? (2) I used to guess that the blade airfoil kept the same within each "DRNodes" as indicated in the following figure. But I am not sure whether the gusess is correct. (3) I change BlAFID to 9 for enabling UserProp at BlSpn = 47.15 m. With this setting, I am confused about the span location range which enables UserProp. From 43.05 m to 47.15 m ? Or the "DRNodes" containg BlSpn = 47.15 m? (4) As you mentioned above, the aerodynamic and structural discretizations are independent. If I want to employ a trailing edge flap at span location 70-90%, I first change the corresponding aerodynamic discretizations to contain 70-90% span location in AD15 blade files. But the structural discretizations is defined by ElastoDyn input BldNodes. I am not sure whether the span location where ElastoDyn uses the aerodynamic forces with trail edge flaps is 70-90%.

------- AERODYN v15.00.* BLADE DEFINITION INPUT FILE ------------------------------------- NREL 5.0 MW offshore baseline aerodynamic blade input properties; note that we need to add the aerodynamic center to this file ====== Blade Properties ================================================================= 19 NumBlNds - Number of blade nodes used in the analysis (-) BlSpn BlCrvAC BlSwpAC BlCrvAng BlTwist BlChord BlAFID BlCb BlCenBn BlCenBt (m) (m) (m) (deg) (deg) (m) (-) (-) (m) (m) 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 1.3308000E+01 3.5420000E+00 1 0.0 0.0 0.0 1.3667000E+00 -8.1531745E-04 -3.4468858E-03 0.0000000E+00 1.3308000E+01 3.5420000E+00 1 0.0 0.0 0.0 4.1000000E+00 -2.4839790E-02 -1.0501421E-01 0.0000000E+00 1.3308000E+01 3.8540000E+00 1 0.0 0.0 0.0 6.8333000E+00 -5.9469375E-02 -2.5141635E-01 0.0000000E+00 1.3308000E+01 4.1670000E+00 2 0.0 0.0 0.0 1.0250000E+01 -1.0909141E-01 -4.6120149E-01 0.0000000E+00 1.3308000E+01 4.5570000E+00 3 0.0 0.0 0.0 1.4350000E+01 -1.1573354E-01 -5.6986665E-01 0.0000000E+00 1.1480000E+01 4.6520000E+00 4 0.0 0.0 0.0 1.8450000E+01 -9.8316709E-02 -5.4850833E-01 0.0000000E+00 1.0162000E+01 4.4580000E+00 4 0.0 0.0 0.0 2.2550000E+01 -8.3186967E-02 -5.2457001E-01 0.0000000E+00 9.0110000E+00 4.2490000E+00 5 0.0 0.0 0.0 2.6650000E+01 -6.7933232E-02 -4.9624675E-01 0.0000000E+00 7.7950000E+00 4.0070000E+00 6 0.0 0.0 0.0 3.0750000E+01 -5.3393159E-02 -4.6544755E-01 0.0000000E+00 6.5440000E+00 3.7480000E+00 6 0.0 0.0 0.0 3.4850000E+01 -4.0899260E-02 -4.3583519E-01 0.0000000E+00 5.3610000E+00 3.5020000E+00 7 0.0 0.0 0.0 3.8950000E+01 -2.9722933E-02 -4.0591323E-01 0.0000000E+00 4.1880000E+00 3.2560000E+00 7 0.0 0.0 0.0 4.3050000E+01 -2.0511081E-02 -3.7569051E-01 0.0000000E+00 3.1250000E+00 3.0100000E+00 8 0.0 0.0 0.0 4.7150000E+01 -1.3980013E-02 -3.4521705E-01 0.0000000E+00 2.3190000E+00 2.7640000E+00 9 0.0 0.0 0.0 5.1250000E+01 -8.3819737E-03 -3.1463837E-01 0.0000000E+00 1.5260000E+00 2.5180000E+00 8 0.0 0.0 0.0 5.4666700E+01 -4.3546914E-03 -2.8909220E-01 0.0000000E+00 8.6300000E-01 2.3130000E+00 8 0.0 0.0 0.0 5.7400000E+01 -1.6838383E-03 -2.6074456E-01 0.0000000E+00 3.7000000E-01 2.0860000E+00 8 0.0 0.0 0.0 6.0133300E+01 -3.2815226E-04 -1.7737470E-01 0.0000000E+00 1.0600000E-01 1.4190000E+00 8 0.0 0.0 0.0 6.1499900E+01 -3.2815226E-04 -1.7737470E-01 0.0000000E+00 1.0600000E-01 1.4190000E+00 8 0.0 0.0 0.0

image

jjonkman commented 1 year ago

Dear @ietqlw,

Regarding your questions 1-3, AeroDyn will only compute aerodynamic loads at the aerodynamic analysis nodes. In AeroDyn v14 and earlier, the aerodynamic analysis nodes were located at the centers of blade elements and the aerodynamic load was assumed to be constant across an element. In AeroDyn v15, the aerodynamic analysis nodes are located at the end points of blade elements and the aerodynamic loads are assumed to be linearly varying along the element between its two end points.

Regarding you last question, with AeroDyn v14 and earlier, the aerodynamic discretization and structural discretization in ElastoDyn where identical and the load transfer was one-to-one. With AeroDyn v15, the aerodynamic and structural discretization are independent and spatial mesh-mapping is used to transfer the aerodynamic loads along the AeroDyn blade analysis nodes to those applied at the structural analysis nodes.

Best regards,

ietqlw commented 1 year ago

Dear @jjonkman,

Thank you for your clear answer.

Inspired by your comments, I can conclude that there are some differences between AD14 and AD15 regarding selecting the span location range of the UserProp properties. I would like to share our views about the differences. Please determine whether my views are correct.

(1) For AD14, RNodes is the distance relative to the hub center. The aerodynamic analysis nodes were located at the centers of blade elements and the aerodynamic load was assumed to be constant across an element. We change the NFoils of three Nodes from 8 to 9 for employ UserProp properties. Hence, the UserProp properties are allowed from 48.65-4.1/2 m to 56.1667+2.7333/2 m relative to the hub center.

--------- AeroDyn v14.04.* INPUT FILE ------------------------------------------------------------------------- NREL 5.0 MW offshore baseline aerodynamic input properties. "BEDDOES" StallMod - Dynamic stall included [BEDDOES or STEADY] (unquoted string) "USE_CM" UseCm - Use aerodynamic pitching moment model? [USE_CM or NO_CM] (unquoted string) "EQUIL" InfModel - Inflow model [DYNIN or EQUIL] (unquoted string) "SWIRL" IndModel - Induction-factor model [NONE or WAKE or SWIRL] (unquoted string) 0.005 AToler - Induction-factor tolerance (convergence criteria) (-) "PRANDtl" TLModel - Tip-loss model (EQUIL only) [PRANDtl, GTECH, or NONE] (unquoted string) "PRANDtl" HLModel - Hub-loss model (EQUIL only) [PRANdtl or NONE] (unquoted string) "NEWTOWER" TwrShad - INSTEAD OF: 0.0 TwrShad - Tower-shadow velocity deficit (-) True TwrPotent - Calculate tower potential flow (flag) INSTEAD OF 9999.9 ShadHWid - Tower-shadow half width (m) False TwrShadow - Calculate tower shadow (flag) INSTEAD OF 9999.9 T_Shad_Refpt- Tower-shadow reference point (m) "NRELOffshrBsline5MW_Onshore_AeroDyn_Tower.dat" TwrFile - Tower drag file name (quoted string) True CalcTwrAero - Calculate aerodynamic drag of the tower at the ElastoDyn nodes. TwrPotent must be true. 1.225 AirDens - Air density (kg/m^3) 1.464E-05 KinVisc - Kinematic air viscosity [CURRENTLY IGNORED] (m^2/sec) "default" DTAero - Time interval for aerodynamic calculations (sec) 8 NumFoil - Number of airfoil files (-) "AeroData/Cylinder1.dat" FoilNm - Names of the airfoil files [NumFoil lines] (quoted strings) "AeroData/Cylinder2.dat" "AeroData/DU40_A17.dat" "AeroData/DU35_A17.dat" "AeroData/DU30_A17.dat" "AeroData/DU25_A17.dat" "AeroData/DU21_A17.dat" "AeroData/NACA64_A17.dat" 17 BldNodes - Number of blade nodes used for analysis (-) RNodes AeroTwst DRNodes Chord NFoil PrnElm 2.8667000E+00 1.3308000E+01 2.7333000E+00 3.5420000E+00 1 NOPRINT 5.6000000E+00 1.3308000E+01 2.7333000E+00 3.8540000E+00 1 NOPRINT 8.3333000E+00 1.3308000E+01 2.7333000E+00 4.1670000E+00 2 NOPRINT 1.1750000E+01 1.3308000E+01 4.1000000E+00 4.5570000E+00 3 NOPRINT 1.5850000E+01 1.1480000E+01 4.1000000E+00 4.6520000E+00 4 NOPRINT 1.9950000E+01 1.0162000E+01 4.1000000E+00 4.4580000E+00 4 NOPRINT 2.4050000E+01 9.0110000E+00 4.1000000E+00 4.2490000E+00 5 NOPRINT 2.8150000E+01 7.7950000E+00 4.1000000E+00 4.0070000E+00 6 NOPRINT 3.2250000E+01 6.5440000E+00 4.1000000E+00 3.7480000E+00 6 NOPRINT 3.6350000E+01 5.3610000E+00 4.1000000E+00 3.5020000E+00 7 NOPRINT 4.0450000E+01 4.1880000E+00 4.1000000E+00 3.2560000E+00 7 NOPRINT 4.4550000E+01 3.1250000E+00 4.1000000E+00 3.0100000E+00 8 NOPRINT 4.8650000E+01 2.3190000E+00 4.1000000E+00 2.7640000E+00 9 NOPRINT ! we change NFoil to 9 5.2750000E+01 1.5260000E+00 4.1000000E+00 2.5180000E+00 9 NOPRINT ! we change NFoil to 9 5.6166700E+01 8.6300000E-01 2.7333000E+00 2.3130000E+00 9 NOPRINT ! we change NFoil to 9 5.8900000E+01 3.7000000E-01 2.7333000E+00 2.0860000E+00 8 NOPRINT 6.1633300E+01 1.0600000E-01 2.7333000E+00 1.4190000E+00 8 NOPRINT User MultiTab - multiple airfoil table option ("ReNum","User","Single")

(2) For AD15, BlSpn is the distance relative to the blade root. The aerodynamic analysis nodes are located at the end points of blade elements and the aerodynamic loads are assumed to be linearly varying along the element between its two end points. If I want to employ UserProp properties from 47.15 m to 57.4 m relative to the blade root, I try to add two extra nodes (BlSpn = 47.1499 m, and BlSpn = 57.4001 m) respectively. Then I change NumBlNds from 19 to 21. I wonder whether my modification of the blade properties is correct. Is there any negetive effect on load analysis because the element from BlSpn =47.1499 m to 47.15 m is too short and the aerodynamic loads varies too largely due to different BlAFID?

    ------- AERODYN v15.00.* BLADE DEFINITION INPUT FILE -------------------------------------

NREL 5.0 MW offshore baseline aerodynamic blade input properties; note that we need to add the aerodynamic center to this file ====== Blade Properties ================================================================= 21 NumBlNds - Number of blade nodes used in the analysis (-) !add two nodes close to the airfoil sections with UserProp BlSpn BlCrvAC BlSwpAC BlCrvAng BlTwist BlChord BlAFID BlCb BlCenBn BlCenBt (m) (m) (m) (deg) (deg) (m) (-) (-) (m) (m) 0.0000000E+00 0.0000000E+00 0.0000000E+00 0.0000000E+00 1.3308000E+01 3.5420000E+00 1 0.0 0.0 0.0 1.3667000E+00 -8.1531745E-04 -3.4468858E-03 0.0000000E+00 1.3308000E+01 3.5420000E+00 1 0.0 0.0 0.0 4.1000000E+00 -2.4839790E-02 -1.0501421E-01 0.0000000E+00 1.3308000E+01 3.8540000E+00 1 0.0 0.0 0.0 6.8333000E+00 -5.9469375E-02 -2.5141635E-01 0.0000000E+00 1.3308000E+01 4.1670000E+00 2 0.0 0.0 0.0 1.0250000E+01 -1.0909141E-01 -4.6120149E-01 0.0000000E+00 1.3308000E+01 4.5570000E+00 3 0.0 0.0 0.0 1.4350000E+01 -1.1573354E-01 -5.6986665E-01 0.0000000E+00 1.1480000E+01 4.6520000E+00 4 0.0 0.0 0.0 1.8450000E+01 -9.8316709E-02 -5.4850833E-01 0.0000000E+00 1.0162000E+01 4.4580000E+00 4 0.0 0.0 0.0 2.2550000E+01 -8.3186967E-02 -5.2457001E-01 0.0000000E+00 9.0110000E+00 4.2490000E+00 5 0.0 0.0 0.0 2.6650000E+01 -6.7933232E-02 -4.9624675E-01 0.0000000E+00 7.7950000E+00 4.0070000E+00 6 0.0 0.0 0.0 3.0750000E+01 -5.3393159E-02 -4.6544755E-01 0.0000000E+00 6.5440000E+00 3.7480000E+00 6 0.0 0.0 0.0 3.4850000E+01 -4.0899260E-02 -4.3583519E-01 0.0000000E+00 5.3610000E+00 3.5020000E+00 7 0.0 0.0 0.0 3.8950000E+01 -2.9722933E-02 -4.0591323E-01 0.0000000E+00 4.1880000E+00 3.2560000E+00 7 0.0 0.0 0.0 4.3050000E+01 -2.0511081E-02 -3.7569051E-01 0.0000000E+00 3.1250000E+00 3.0100000E+00 8 0.0 0.0 0.0 4.7149900E+01 -1.3980013E-02 -3.4521705E-01 0.0000000E+00 2.3190000E+00 2.7640000E+00 8 0.0 0.0 0.0 !add a node close to the airfoil sections with UserProp 4.7150000E+01 -1.3980013E-02 -3.4521705E-01 0.0000000E+00 2.3190000E+00 2.7640000E+00 9 0.0 0.0 0.0 !BlAFID to 9 5.1250000E+01 -8.3819737E-03 -3.1463837E-01 0.0000000E+00 1.5260000E+00 2.5180000E+00 9 0.0 0.0 0.0 !BlAFID to 9 5.4666700E+01 -4.3546914E-03 -2.8909220E-01 0.0000000E+00 8.6300000E-01 2.3130000E+00 9 0.0 0.0 0.0 !BlAFID to 9 5.7400000E+01 -1.6838383E-03 -2.6074456E-01 0.0000000E+00 3.7000000E-01 2.0860000E+00 9 0.0 0.0 0.0 !BlAFID to 9 5.7400100E+01 -1.6838383E-03 -2.6074456E-01 0.0000000E+00 3.7000000E-01 2.0860000E+00 8 0.0 0.0 0.0
!add a node close to the airfoil sections with UserProp 6.0133300E+01 -3.2815226E-04 -1.7737470E-01 0.0000000E+00 1.0600000E-01 1.4190000E+00 8 0.0 0.0 0.0 6.1499900E+01 -3.2815226E-04 -1.7737470E-01 0.0000000E+00 1.0600000E-01 1.4190000E+00 8 0.0 0.0 0.0

!bjj: because of precision in the BD-AD coupling, 61.5m didn't work, so I changed it to 61.4999m 6.1500000E+01 -3.2815226E-04 -1.7737470E-01 0.0000000E+00 1.0600000E-01 1.4190000E+00 8 0.0 0.0 0.0

Sincerely.

jjonkman commented 1 year ago

Dear @ietqlw,

I agree with your approaches.

In AeroDyn v15 you can have nodes close together as you indicated. My only concern would be if the nodes are so close together that they would be within the internal tolerances used by the spatial mesh mapping.

Best regards,

ietqlw commented 1 year ago

Dear @jjonkman,

Thank you for your patient answer.

As you mentioned above, error was reported because the nodes were too close. I changed 47.1499 m to 47.1490 m and 57.4001 m to 57.4010 m. The error was then addressed.

Sincerely