Bug fix: Previously, the user could specify element normal and tangent vectors that make the spanwise vector (and element circulation) incorrectly oriented in the wake grid. Added code to reconcile the element spanwise vector direction with the wake grid spanwise direction.

Author: jcokemurray

mod/element.f95

@@ -55,7 +55,8 @@ MODULE element
 	real, allocatable :: tzBC(:)		! Tangential Z for each blade segment 	
         real, allocatable :: sxBC(:)		! Spanwise X for each blade segment 	
 	real, allocatable :: syBC(:)		! Spanwise Y for each blade segment 	
-	real, allocatable :: szBC(:)		! Spanwise Z for each blade segment 	
+	real, allocatable :: szBC(:)		! Spanwise Z for each blade segment 
+        real, allocatable :: CircSign(:)        ! Direction of segment circulation on wake grid	at positive lift
 	real, allocatable :: eArea(:)		! Element area to radius ratio for each element
 	real, allocatable :: eChord(:)		! Element chord to radius ratio for each element
         integer, allocatable :: iSect(:)        ! Array of indicies of the section table to apply to each blade element
@@ -116,7 +117,8 @@ SUBROUTINE element_cns(MaxSegEnds,MaxSegEndPerBlade)
 		allocate(tzBC(MaxSegEnds))		
 		allocate(sxBC(MaxSegEnds))		
 		allocate(syBC(MaxSegEnds))		
-		allocate(szBC(MaxSegEnds))	        
+		allocate(szBC(MaxSegEnds))
+                allocate(CircSign(MaxSegEnds))              	        
 		allocate(eArea(MaxSegEnds))	
 		allocate(eChord(MaxSegEnds))			
                 allocate(iSect(MaxSegEnds))              

src/BGeomSetup.f95

@@ -77,6 +77,11 @@ SUBROUTINE BGeomSetup()
                         yBC(nej)=0.5*(Blades(i)%QCy(j)+Blades(i)%QCy(j+1))                                                                                              
                         zBC(nej)=0.5*(Blades(i)%QCz(j)+Blades(i)%QCz(j+1)) 
 
+                        ! QC line 
+                        dx=Blades(i)%QCx(j+1)-Blades(i)%QCx(j)
+                        dy=Blades(i)%QCy(j+1)-Blades(i)%QCy(j)
+                        dz=Blades(i)%QCz(j+1)-Blades(i)%QCz(j)
+                        
                         ! Set element tangent vectors, txBC(MaxSegEnd)  
                         txBC(nej)=0.5*(Blades(i)%tx(j)+Blades(i)%tx(j+1))
                         tyBC(nej)=0.5*(Blades(i)%ty(j)+Blades(i)%ty(j+1))                                                                             
@@ -102,13 +107,23 @@ SUBROUTINE BGeomSetup()
                         nyBC(nej)=nyBC(nej)/VMag    
                         nzBC(nej)=nzBC(nej)/VMag 
 
-                        ! Calculate spanwise vector (s = t x n)
+                        ! Calculate element spanwise vector (s = t x n)
                         CALL cross(txBC(nej),tyBC(nej),tzBC(nej),nxBC(nej),nyBC(nej),nzBC(nej),sxBC(nej),syBC(nej),szBC(nej))     
                         ! Force normalize
                         VMag=sqrt(sxBC(nej)**2+syBC(nej)**2+szBC(nej)**2) 
                         sxBC(nej)=sxBC(nej)/VMag 
                         syBC(nej)=syBC(nej)/VMag    
                         szBC(nej)=szBC(nej)/VMag 
+                        ! Calc projection direction of element spanwise vector on QC line.
+                        ! Used to correctly orient circulation in wake grid.
+                        ! Note pos. lift creates circulation in neg. element spanwise direction
+                        dp=sxBC(nej)*dx+syBC(nej)*dy+szBC(nej)*dz
+                        if (dp>0) then
+                            CircSign(nej)=-1.0
+                        else
+                            CircSign(nej)=1.0
+                        end if
+                        
 
 		end do 
 

src/bsload.f95

@@ -20,7 +20,7 @@ SUBROUTINE bsload(nElem,IsBE,alpha,Re,umach,ur,CN,CT,Fx,Fy,Fz,te)
         real uAve, vAve, wAve, uFSAve, vFSAve, wFSAve, uBlade, vBlade, wBlade, urdn, urdc
         real xe5, ye5, ze5, xe75, ye75, ze75, uBlade5, vBlade5, wBlade5, uBlade75, vBlade75, wBlade75
         real urdn5, ur5, alpha5, Re5, urdn75, ur75, alpha75, Re75     
-        real CL, CD, CM25, FN, FT, MS, TRx, TRy, TRz, CLCirc   
+        real CL, CD, CM25, FN, FT, MS, TRx, TRy, TRz, CLCirc, CircDir
 
 
         ! Calculates aero loads on a blade element. Static and dynamic airfoil characteristics calculated here...
@@ -50,6 +50,9 @@ SUBROUTINE bsload(nElem,IsBE,alpha,Re,umach,ur,CN,CT,Fx,Fy,Fz,te)
         sye=syBC(nElem)                                                  
         sze=szBC(nElem)
 
+        ! Direction of circulation in wake grid at positive lift
+        CircDir=CircSign(nElem)
+        
         ! Airfoil section                                                 
         SectInd=isect(nElem)                                                     
 
@@ -114,7 +117,7 @@ SUBROUTINE bsload(nElem,IsBE,alpha,Re,umach,ur,CN,CT,Fx,Fy,Fz,te)
 
         ! Bound vortex strength from CL via Kutta-Joukowski analogy. 
         ! Save corresponding AOA as well                                                                                          
-        GB(nElem1)=CLCirc*ElemChordR*ur/2.0  
+        GB(nElem1)=CircDir*(CLCirc*ElemChordR*ur/2.0)
         AOA(nElem1)=alpha
         ! normalized time step used to update states in the LB model
         ds(nElem)=2.0*ur*DT/ElemChordR