Skip to content

Commit 45796c8

Browse files
authored
Merge pull request #2945 from ebranlard/f/srcPanels
[OLAF] Implementation of a source panel method
2 parents da4502a + 66acdf2 commit 45796c8

File tree

19 files changed

+2573
-454
lines changed

19 files changed

+2573
-454
lines changed

docs/source/acknowledgements.rst

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,8 +11,37 @@ from U.S. Department of Energy Wind Energy Technology Office through the
1111
NREL gratefully acknowledges development contributions from the following
1212
organizations:
1313

14-
- Envision Energy USA, Ltd
14+
- Accelerate Wind
15+
- AirLoom Energy
16+
- ARPA-E
1517
- Brigham Young University
18+
- CENER
19+
- Delft University of Technology
20+
- DOE Water Power Technologies Office (WPTO)
21+
- Equinor
22+
- Envision Energy USA, Ltd
23+
- Institute for Energy Technology
24+
- GoogleX
25+
- Hayman Consulting
26+
- Norwegian University of Science and Technology
27+
- NOWRDC
28+
- Octue
29+
- Orcina
30+
- Oregon State University
31+
- Principle Power
32+
- RRD Engineering
33+
- Sandia National Laboratory
34+
- SENER
35+
- Shell
36+
- SOWENTO
37+
- Technical University of Denmark
38+
- TODA Corporation
39+
- TotalEnergies
40+
- University of Strathclyde
41+
- University of Stuttgart
42+
- University of Texas-Austin
43+
- Vestas
44+
- Windward Engineering
1645

1746
NREL gratefully acknowledges additional development support through designation
1847
as an `Intel® Parallel Computing Center (IPCC) <https://software.intel.com/en-us/ipcc>`_.

docs/source/user/aerodyn-olaf/ExampleFiles/ExampleFile--OLAF.dat

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,4 +43,8 @@ default VTK_fps - Frame rate for VTK output (frames per second) {"all
4343
0 nGridOut - Number of grid outputs
4444
GridName GridType TStart TEnd DTGrid XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ
4545
(-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-)
46-
------------------------------------------------------------------------------------------------
46+
===============================================================================================
47+
--------------------------- ADVANCED OPTIONS --------------------------------------------------
48+
===============================================================================================
49+
! Advanced options may be placed here in arbitrary order, with the regular format:
50+
! Value Key - Comment

docs/source/user/aerodyn-olaf/InputFiles.rst

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -312,6 +312,118 @@ of a box of shape 5x20x30 and dimension 1200x300x295. The grid contains both th
312312
The two other grids are vertical and horizontal planes containing only the velocity.
313313

314314

315+
Advanced Options
316+
~~~~~~~~~~~~~~~~
317+
318+
319+
Advanced options (typically used for developers or beta features) can be placed at the end of the OLAF input file:
320+
321+
322+
- These options use the regular format: `Value Key - Comment`.
323+
- They can be placed in arbitrary order.
324+
- They can be commented out with a `!` character at the beginning of the line.
325+
- Blank lines or unsupported options are ignored (display to screen).
326+
- If not provided, a default value will be used, but the "DEFAULT" keyword is not supported for the advanced options.
327+
- These are **beta features** and should mostly be used by developers.
328+
329+
330+
The end of the OLAF input file would look as follows:
331+
332+
.. code::
333+
334+
[...]
335+
GridName GridType TStart TEnd DTGrid XStart XEnd nX YStart YEnd nY ZStart ZEnd nZ
336+
(-) (-) (s) (s) (s) (m) (m) (-) (m) (m) (-) (m) (m) (-)
337+
===============================================================================================
338+
--------------------------- ADVANCED OPTIONS --------------------------------------------------
339+
===============================================================================================
340+
! Advanced options may be placed here in arbitrary order, with the regular format:
341+
! Value1 Key1 - Comment1
342+
! Value2 Key2 - Comment2
343+
! Lines starting with `!` are ignored, empty lines are ignored
344+
[...] etc
345+
346+
Currently, the advanced options supported are as follows:
347+
348+
.. code::
349+
350+
===============================================================================================
351+
--------------------------- ADVANCED OPTIONS --------------------------------------------------
352+
===============================================================================================
353+
"Panels.vtk" SrcPnlFile - Name of VTK file containing source panels {default: ""}
354+
1 nSrcPnlUpdate - How often do src panel updates (in time steps of OLAF), {default: 1}
355+
True Induction - Compute induced velocities from wake to blade and wake to wake, {default: True}
356+
True InductionAtCP - Compute induced velocities at nodes or CP, {default: True}
357+
True WakeAtTE - Start the wake at the trailing edge, or at the LL, {default: True}
358+
False DStallOnWake - Dynamic stall has influence on wake, {default: False}
359+
0.75 kFrozenNWStart - Fraction of wake induced velocity at start of frozen wake, {default: 0.75}
360+
0.5 kFrozenNWEnd - Fraction of wake induced velocity at end of frozen wake, {default: 0.5}
361+
0.0 zGround - Ground height, used to enforce that no vortices go into the ground {default: 0.0}
362+
0.1 zGroundPush - Ground push back, vortices that are lower than zGround are placed back at zGroundPush {default: 0.1}
363+
364+
365+
366+
**SrcPanlFile** [string] specifies the name of the VTK file used for the source panel method.
367+
The VTK file should be a legacy ASCII VTK file, with DATASET POLYDATA POINTS and POLYGONS.
368+
A sample VTK file is provided below for two panels forming a regular grid in the XY plane.
369+
The connectivity of the polygons needs to be such that the normal points away from the body and
370+
into the fluid.
371+
In the example below, the polygons are defined clockwise when viewed from above, which results in the normals internally defined by OLAF pointing in the `+z` direction. This configuration is typical of a bottom wall with fluid above it: the OLAF normals point from the wall toward the fluid. However, when applying the right-hand-rule convention, the resulting normals point from the fluid toward the interior of the body.
372+
When using **WrVTK**, OLAF will write separate VTK files containing various information related to the source panels, such as the pressure coefficient, area, force per area, normals. Looking at the orientation of the normals is extremely important (In Paraview, select 3D Glyph, Arrows, and select the `Normals` output, scaled by the `Normals`).
373+
The BodyID CELL_DATA can be used to separate different patches, which can help the post processing. This is optional, OLAF doesn't use it currently, but it is written back in the output files of OLAF.
374+
375+
Curious readers may look at the unittest `Test_SrcPnl_Sphere` of OLAF that test the pressure coefficient
376+
about a sphere.
377+
378+
.. code::
379+
380+
# vtk DataFile Version 2.0
381+
Comment
382+
ASCII
383+
DATASET POLYDATA
384+
POINTS 6 double
385+
0.0 0.0 0.0
386+
0.0 25.0 0.0
387+
0.0 50.0 0.0
388+
10.0 0.0 0.0
389+
10.0 25.0 0.0
390+
10.0 50.0 0.0
391+
392+
POLYGONS 2 10
393+
4 0 1 4 3
394+
4 1 2 5 4
395+
396+
CELL_DATA 2
397+
SCALARS BodyID int
398+
LOOKUP_TABLE default
399+
0
400+
0
401+
402+
403+
404+
**nSrcPnlUpdate** [int] Define how often the source panel updates (in time steps of OLAF), the default is `1`, at each time step.
405+
406+
**Induction** [switch] Compute induced velocities, otherwise all induced velocity are 0 (no wake, no panels, etc) ! Default is `True`.
407+
408+
**InductionAtCP** [switch] When performing the lifting line calculations, compute induced velocities at nodes (False) or at control points (CP, True). Default is `True`.
409+
410+
**WakeAtTE** [switch] Start the wake at the trailing edge (True), or directly at the LL (False, no chordwise panel). Default is `True`.
411+
412+
**DStallOnWake** [switch] Include the influence of the dynamic stall on the wake (True), i.e. use the dynamic Cl to update the circulation of the lifting line and wake panel. Default is `False`.
413+
414+
**kFrozenNWStart** [float] Fraction of wake induced velocity at start of frozen wake, default is `0.75`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`.
415+
416+
**kFrozenNWEnd** [float] Fraction of wake induced velocity at end of frozen wake, default is `0.5`. See OLAF theory related to Frozen NW, :numref:`sec:vortconvfrozen`.
417+
418+
**zGround** [float] Height in meters below which vortex points are not allowed to be and if any are present they will be pushed back above the ground at the height defined by **zGroundPush** (in meters). Default is `0.0`.
419+
For MHK, the sea bed location, based on the water depth is added to **zGround**.
420+
421+
**zGroundPush** [float] Ground push back, see **zGround**. Default is `0.1`.
422+
423+
424+
425+
426+
315427
AeroDyn Input File
316428
--------------------
317429
Input file modifications

modules/aerodyn/src/AeroDyn.f90

Lines changed: 26 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -412,7 +412,12 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
412412

413413

414414
! set the rest of the parameters
415-
p%Skew_Mod = InputFileData%Skew_Mod
415+
! NOTE: some parameters need to be set even if no rotors
416+
p%Skew_Mod = InputFileData%Skew_Mod
417+
p%UA_Flag = InputFileData%UA_Init%UAMod > UA_None
418+
p%CompAeroMaps = InitInp%CompAeroMaps
419+
p%DT = InputFileData%DTAero
420+
p%Wake_Mod = InputFileData%Wake_Mod
416421
do iR = 1, nRotors
417422
p%rotors(iR)%AeroProjMod = AeroProjMod(iR)
418423
call WrScr(' AeroDyn: projMod: '//trim(num2lstr(p%rotors(iR)%AeroProjMod)))
@@ -1472,9 +1477,8 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err
14721477

14731478
p%MHK = InitInp%MHK
14741479

1475-
p_AD%DT = InputFileData%DTAero
1476-
p_AD%Wake_Mod = InputFileData%Wake_Mod
14771480
p%DBEMT_Mod = InputFileData%DBEMT_Mod
1481+
14781482
p%TwrPotent = InputFileData%TwrPotent
14791483
p%TwrShadow = InputFileData%TwrShadow
14801484
p%TwrAero = InputFileData%TwrAero
@@ -1765,9 +1769,11 @@ subroutine AD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg )
17651769
if (p%Wake_Mod == WakeMod_FVW ) then
17661770

17671771
if ( p%UA_Flag ) then
1768-
do iW=1,p%FVW%nWings
1769-
call UA_End(m%FVW%W(iW)%p_UA)
1770-
enddo
1772+
if (allocated(m%FVW%W)) then
1773+
do iW=1,p%FVW%nWings
1774+
call UA_End(m%FVW%W(iW)%p_UA)
1775+
enddo
1776+
endif
17711777
end if
17721778

17731779
call FVW_End( m%FVW_u, p%FVW, x%FVW, xd%FVW, z%FVW, OtherState%FVW, m%FVW_y, m%FVW, ErrStat, ErrMsg )
@@ -5158,11 +5164,25 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta
51585164
allocate(InitInp%W(nWings) , STAT = ErrStat2); ErrMsg2='Allocate W'; if(Failed()) return
51595165
allocate(InitInp%WingsMesh(nWings), STAT = ErrStat2); ErrMsg2='Allocate Wings Mesh'; if(Failed()) return
51605166

5167+
! --- Default inputs (not per retor)
5168+
! UA and inputs that are not per rotor
5169+
InitInp%UA_Flag = p%UA_Flag
5170+
! Important inputs if no rotors are present!
5171+
if (size(p%rotors)==0) then
5172+
InitInp%numBladeNodes = 0
5173+
InitInp%AirDens = InputFileData%AirDens
5174+
InitInp%KinVisc = InputFileData%KinVisc
5175+
InitInp%MHK = 0 ! TODO this is an initinp of AeroDyn
5176+
InitInp%WtrDpth = 0.0_ReKi ! TODO this is an initinp of AeroDyn
5177+
InitInp%RootName = p%RootName(1:len_trim(p%RootName)-2) ! Removing "AD"
5178+
endif
5179+
51615180
! --- Inputs per wings/blades
51625181
iW_incr=0
51635182
do iR=1, size(p%rotors)
51645183

51655184
InitInp%numBladeNodes = p%rotors(iR)%numBlNds ! TODO TODO TODO per wing
5185+
InitInp%AirDens = p%rotors(iR)%AirDens
51665186
InitInp%KinVisc = p%rotors(iR)%KinVisc
51675187
InitInp%MHK = p%rotors(iR)%MHK
51685188
InitInp%WtrDpth = p%rotors(iR)%WtrDpth
@@ -5225,9 +5245,6 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta
52255245
if(Failed()) return
52265246

52275247
enddo ! iB, blades
5228-
5229-
! Unsteady Aero Data
5230-
InitInp%UA_Flag = p%UA_Flag
52315248
call UA_CopyInitInput(InputFileData%UA_Init, InitInp%UA_Init, MESH_NEWCOPY, ErrStat2, ErrMsg2)
52325249

52335250
iW_incr = iW_incr+p%rotors(iR)%numBlades

modules/aerodyn/src/AeroDyn_Driver_Subs.f90

Lines changed: 19 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -541,7 +541,7 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, needInitIW, errStat, err
541541
InitInp%AD%defPvap = dvr%Pvap
542542
InitInp%AD%WtrDpth = dvr%WtrDpth
543543
InitInp%AD%MSL2SWL = dvr%MSL2SWL
544-
InitInp%AD%CompSeaSt = dvr%SS_InitInp%CompSeaSt
544+
InitInp%AD%CompSeaSt = dvr%SS_InitInp%CompSeaSt /=0
545545
! Init data per rotor
546546
allocate(InitInp%AD%rotors(dvr%numTurbines), stat=errStat)
547547
if (errStat/=0) then
@@ -1038,7 +1038,7 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg )
10381038
call ParseCom(FileInfo_In, CurLine, Line, errStat2, errMsg2, unEc); if (Failed()) return
10391039
call ParseVar(FileInfo_In, CurLine, "compInflow", dvr%IW_InitInp%compInflow , errStat2, errMsg2, unEc); if (Failed()) return
10401040
call ParseVar(FileInfo_In, CurLine, "InflowFile", dvr%IW_InitInp%InputFile, errStat2, errMsg2, unEc, IsPath=.true.); if (Failed()) return
1041-
if (dvr%IW_InitInp%compInflow==0) then
1041+
if (dvr%IW_InitInp%compInflow/=1) then
10421042
call ParseVar(FileInfo_In, CurLine, "HWindSpeed", dvr%IW_InitInp%HWindSpeed , errStat2, errMsg2, unEc); if (Failed()) return
10431043
call ParseVar(FileInfo_In, CurLine, "RefHt" , dvr%IW_InitInp%RefHt , errStat2, errMsg2, unEc); if (Failed()) return
10441044
call ParseVar(FileInfo_In, CurLine, "PLExp" , dvr%IW_InitInp%PLExp , errStat2, errMsg2, unEc); if (Failed()) return
@@ -1544,6 +1544,9 @@ subroutine Dvr_InitializeOutputs(nWT, out, numSteps, errStat, errMsg)
15441544
end do ! i
15451545
write (out%unOutFile(iWT),'()')
15461546
enddo
1547+
if (nWT==0) then
1548+
! Special case, no turbine, TODO
1549+
endif
15471550
endif
15481551

15491552
! --- Binary
@@ -1745,7 +1748,11 @@ subroutine Dvr_WriteOutputs(nt, t, dvr, out, yADI, SeaSt, errStat, errMsg)
17451748
END IF
17461749

17471750
! Packing all outputs excpet time into one array
1748-
nAD = size(yADI%AD%rotors(1)%WriteOutput)
1751+
if (size(yADI%AD%rotors)>=1) then
1752+
nAD = size(yADI%AD%rotors(1)%WriteOutput)
1753+
else
1754+
nAD = 0
1755+
endif
17491756
nIW = size(yADI%IW_WriteOutput)
17501757
nSS = size(SeaSt%y%WriteOutput)
17511758
nDV = out%nDvrOutputs
@@ -1873,11 +1880,15 @@ subroutine setVTKParameters(DVR_Outs, dvr, ADI, errStat, errMsg, dirname)
18731880
! Get radius for ground (blade length + hub radius):
18741881
GroundRad = MaxBladeLength + MaxTwrLength+ DVR_Outs%VTKHubRad
18751882
! write the ground or seabed reference polygon:
1876-
RefPoint(1:2) = dvr%WT(1)%originInit(1:2)
1877-
do iWT=2,dvr%numTurbines
1878-
RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2)
1879-
end do
1880-
RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines
1883+
if (dvr%numTurbines>0) then
1884+
RefPoint(1:2) = dvr%WT(1)%originInit(1:2)
1885+
do iWT=2,dvr%numTurbines
1886+
RefPoint(1:2) = RefPoint(1:2) + dvr%WT(iWT)%originInit(1:2)
1887+
end do
1888+
RefPoint(1:2) = RefPoint(1:2) / dvr%numTurbines
1889+
else
1890+
RefPoint(1:3) = 0.0_ReKi
1891+
endif
18811892

18821893
RefPoint(3) = 0.0_ReKi
18831894
RefLengths = GroundRad + sqrt((WorldBoxMax(1)-WorldBoxMin(1))**2 + (WorldBoxMax(2)-WorldBoxMin(2))**2)

modules/aerodyn/src/AeroDyn_Inflow.f90

Lines changed: 29 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -86,12 +86,18 @@ subroutine ADI_Init(InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
8686
InitOut%Ver = InitOut_AD%ver
8787
! Add writeoutput units and headers to driver, same for all cases and rotors!
8888
!TODO: this header is too short if we add more rotors. Should also add a rotor identifier
89-
call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return
89+
if (size(InitOut_AD%rotors)>=1) then
90+
call concatOutputHeaders(InitOut%WriteOutputHdr, InitOut%WriteOutputUnt, InitOut_AD%rotors(1)%WriteOutputHdr, InitOut_AD%rotors(1)%WriteOutputUnt, errStat2, errMsg2); if(Failed()) return
91+
endif
9092

9193
! --- Initialize grouped outputs
92-
!TODO: assumes one rotor
93-
p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts
94-
call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return
94+
if (size(InitOut_AD%rotors)>=1) then
95+
!TODO: assumes one rotor
96+
p%NumOuts = p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts + m%IW%p%NumOuts
97+
call AllocAry(y%WriteOutput, p%NumOuts, 'WriteOutput', errStat2, errMsg2); if (Failed()) return
98+
else
99+
p%NumOuts = m%IW%p%NumOuts
100+
endif
95101

96102
! --- Initialize outputs
97103
call AllocAry(y%IW_WriteOutput, size(m%IW%y%WriteOutput),'IW_WriteOutput', errStat2, errMsg2); if(Failed()) return
@@ -297,21 +303,26 @@ subroutine ADI_CalcOutput(t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg)
297303
y%PLExp = m%IW%PLExp
298304

299305
! --- Set outputs
300-
!TODO: this assumes one rotor!!!
301-
associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, &
302-
IW_NumOuts => m%IW%p%NumOuts)
303-
y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts)
304-
y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts)
305-
end associate
306-
307-
!----------------------------------------------------------------------------
308-
! Store hub height velocity calculated in CalcOutput
309-
!----------------------------------------------------------------------------
306+
if (size(p%AD%rotors)>=1) then
307+
!TODO: this assumes one rotor!!!
308+
associate(AD_NumOuts => p%AD%rotors(1)%NumOuts + p%AD%rotors(1)%BldNd_TotNumOuts, &
309+
IW_NumOuts => m%IW%p%NumOuts)
310+
y%WriteOutput(1:IW_NumOuts) = y%IW_WriteOutput(1:IW_NumOuts)
311+
y%WriteOutput(IW_NumOuts+1:p%NumOuts) = y%AD%rotors(1)%WriteOutput(1:AD_NumOuts)
312+
end associate
313+
314+
!----------------------------------------------------------------------------
315+
! Store hub height velocity calculated in CalcOutput
316+
!----------------------------------------------------------------------------
317+
318+
if (p%storeHHVel) then
319+
do iWT = 1, size(u%AD%rotors)
320+
y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1)
321+
end do
322+
endif
310323

311-
if (p%storeHHVel) then
312-
do iWT = 1, size(u%AD%rotors)
313-
y%HHVel(:,iWT) = m%AD%Inflow(1)%RotInflow(iWT)%InflowOnHub(:,1)
314-
end do
324+
else
325+
y%WriteOutput(1:p%NumOuts) = y%IW_WriteOutput(1:m%IW%p%NumOuts)
315326
endif
316327

317328
contains

0 commit comments

Comments
 (0)