Added option in wake data output module to calculate streamwise wake deficit on a plane, averaged over the last revolution simulated. The description of this

plane is currently hardcoded in the wakedata module, but the calculation of the output is fairly general. Could eventually generalize the wake data output 
geometry, along with the turbine geometry, with a geometry preprocessor.

JCM

Author: jcokemurray

Doc/Input_Description

@@ -72,7 +72,8 @@ hAG (real): Height above ground at base of blades (ft)
 CTExcrM (real): Additional machine level excrescence torque based on tip speed and Rmax, CTExcrM = TorqueExcr / (1/2*rho*Utip^2*Rmax^3) (default: 0)
 
 Optional outputs:
-WakeOutFlag (int): 1 to write wakeline outputs for the element indicies noted below, 0 otherwise (default).
+WakeOutFlag (int): 1 to write wakeline outputs for the element indicies noted below, 2 to include average wake velocity deficit (on a plane with 
+	definition that's currently hardcoded in wakedata module), 0 otherwise (default).
 	WLI (int): Array (comma separated) of element end indicies for which to write wake line data. 	
 
 

Doc/Output_Description

@@ -57,6 +57,7 @@ te (-):			Torque coefficient from the element, Torque / (1/2*rho*Uinf^2*Af*R), w
 			in the input file).
 LogicFlags:		Twelve logic flags related to the dynamic stall models for model testing... (see BVLogicOutputs and LBLogicOutputs in the code)
 
+
 Wake Data Output (_WakeData.csv):
 
 Timestep:		Timestep counter.
@@ -70,3 +71,9 @@ Z/R:			Z location of the wake node normalized by turbine radius, R, which is the
 U/Uinf:			X velocity of the wake node normalized by freestream speed, Uinf.
 V/Uinf:			Y velocity of the wake node normalized by freestream speed, Uinf.
 W/Uinf:			Z velocity of the wake node normalized by freestream speed, Uinf.
+
+
+Streamwise Wake Deficit Data Output (_WakeDefData.csv):
+	
+2D array of streamwise wake deficit normalized by freestream speed. X (streamwise direction) varies down a column, Z varies across a row.
+

mod/wakedata.f95

@@ -3,10 +3,25 @@ MODULE wakedata
 	! Wake visualization data for WriteWakeData
 
         integer :: WakeOutFlag
-	real, allocatable :: WakeDefAve(:,:)
         integer, allocatable :: WakeLineInd(:) 
         integer :: NWakeInd     
         character(1000) :: WakeOutHead = 'Timestep,Element,X/R,Y/R,Z/R,U/Uinf,V/Uinf,W/Uinf'
+        
+        ! Wake deficit calculation performed if WakeOutFlag=2
+        ! JCM test: wake deficit output plane is currently hardcoded...
+        integer :: nxgrid = 141
+        integer :: nzgrid = 81
+        real :: ygrid = 1.2325
+        real :: xgridL = -2
+        real :: xgridU = 5
+        real :: zgridL = -2
+        real :: zgridU = 2
+        integer :: xcount, zcount, ntcount
+        real :: dxgrid, dzgrid
+        real :: PointGrid(3), IndVelGrid(3)
+        real, allocatable :: XGrid(:,:) 
+        real, allocatable :: ZGrid(:,:) 
+        real, allocatable :: SVDef(:,:) 
 
 
 	CONTAINS
@@ -15,7 +30,12 @@ SUBROUTINE wakedata_cns()
 
 		! Constructor for the arrays in this module
 
-                allocate(WakeLineInd(NWakeInd))              
+                allocate(WakeLineInd(NWakeInd))     
+                         
+                ! JCM test: wake deficit output
+                allocate(XGrid(nxgrid,nzgrid))
+                allocate(ZGrid(nxgrid,nzgrid))
+                allocate(SVDef(nxgrid,nzgrid))         
 
 	End SUBROUTINE
 

src/CACTUS.f95

@@ -55,7 +55,7 @@ PROGRAM CACTUS
         real :: cpave, cpave_last
         real :: delt, delty, deltb, deltr
 
-        character(80) :: InputFN, SFOutputFN, RevOutputFN, TSOutputFN, ELOutputFN, RegOutputFN, WakeOutputFN, FNBase                                   
+        character(80) :: InputFN, SFOutputFN, RevOutputFN, TSOutputFN, ELOutputFN, RegOutputFN, WakeOutputFN, WakeDefOutputFN, FNBase                                   
 
         ! Pi definition
         pi = 4.0*atan(1.0)
@@ -124,9 +124,15 @@ PROGRAM CACTUS
         end if
 
         ! Optional wake line data output
-        if (WakeOutFlag == 1) then
+        if (WakeOutFlag > 0) then
                 WakeOutputFN=trim(FNBase)//'_WakeData.csv'
                 OPEN(12, FILE=WakeOutputFN)
+                
+                if (WakeOutFlag > 1) then
+                        ! wake deficit surface output
+                        WakeDefOutputFN=trim(FNBase)//'_WakeDefData.csv'
+                        OPEN(13, FILE=WakeDefOutputFN)
+                end if
         end if                                                                    
 
         ! If wake update interval set to a negative number, set next wake update iteration to -1 (no wake velocity updates will be performed)
@@ -333,7 +339,7 @@ PROGRAM CACTUS
                         end if 
 
                         ! Write current wake data for viewing in Matlab
-                        if (WakeOutFlag == 1) then
+                        if (WakeOutFlag > 0) then
                                 Call WriteWakeData()
                         end if  
 

src/WriteWakeData.f95

@@ -7,11 +7,26 @@ SUBROUTINE WriteWakeData()
 	use wallsoln 
 	use configr
 
+        implicit none
+        
 	integer :: tCount, tCountMax, wcount
 
         ! Write header
         if (NT==1) then
                 write(12,*) trim(WakeOutHead)
+                
+                ! JCM test: wake deficit setup
+                ntcount=0
+                
+                dxgrid=(xgridU-xgridL)/(nxgrid-1)
+                dzgrid=(zgridU-zgridL)/(nzgrid-1)
+                do xcount=1,nxgrid                                                    
+                        do zcount=1,nzgrid
+                                XGrid(xcount,zcount)=xgridL+(xcount-1)*dxgrid
+                                ZGrid(xcount,zcount)=zgridL+(zcount-1)*dzgrid
+                                SVDef(xcount,zcount)=0.0
+                        end do
+                end do
         end if
 
 	! Write wake positions and velocity for each wake line on last rev
@@ -29,7 +44,42 @@ SUBROUTINE WriteWakeData()
                                 ! Dont suppress carriage return on last column
                                 write(12,'(E13.7)') W(tCount,WakeLineInd(wcount))
 	                end do
-                end do  
+                end do 
+                 
+                
+                if (WakeOutFlag > 1) then
+                
+                        ! Output blade, wake, and wall induced streamwise velocity deficit on a plane.
+                        
+                        ! Averaged over last revolution
+                        do xcount=1,nxgrid                                                    
+                                do zcount=1,nzgrid
+        
+                                        ! Calculate wall and wake induced velocities at grid
+                                        PointGrid=[XGrid(xcount,zcount),ygrid,ZGrid(xcount,zcount)]
+                                        Call CalcIndVel(NT,ntTerm,NBE,NB,NE,PointGrid,IndVelGrid)
+                                        SVDef(xcount,zcount)=SVDef(xcount,zcount)+IndVelGrid(1)/nti
+                                        
+                                end do
+                        end do
+                        ntcount=ntcount+1
+                        
+                        ! Write on last iter
+                        if (ntcount == nti) then
+                        do xcount=1,nxgrid        
+                                        do zcount=1,nzgrid  
+                                                if (zcount < nzgrid) then
+                                                        write(13,'(E13.7,",",$)') SVDef(xcount,zcount)
+                                                else 
+                                                        ! Dont suppress carriage return on last column
+                                                        write(13,'(E13.7)') SVDef(xcount,zcount)
+                                                end if
+                                        end do
+                                end do 
+                        end if
+                
+                end if
+                 
         end if     
 
 Return