diff --git a/glue-codes/fast-farm/src/FASTWrapper.f90 b/glue-codes/fast-farm/src/FASTWrapper.f90 index c2da0ce92b..7643dc8561 100644 --- a/glue-codes/fast-farm/src/FASTWrapper.f90 +++ b/glue-codes/fast-farm/src/FASTWrapper.f90 @@ -535,7 +535,6 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg) REAL(R8Ki) :: theta(3) REAL(R8Ki) :: orientation(3,3) REAL(ReKi) :: tmp_sz_z, tmp_sz_y - REAL(ReKi) :: R_rotor ! Rotor radius INTEGER(IntKi) :: j, k ! loop counters @@ -625,12 +624,6 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg) call setErrStat(ErrStat2,ErrMsg2,ErrStat2,ErrMsg,RoutineName) if (ErrStat >= AbortErrLev) return end do - - ! Rotor radius - R_rotor = 0.0_ReKi - do k=1,size(m%ADRotorDisk) ! loop on blades - R_rotor = max(R_rotor, TwoNorm(m%ADRotorDisk(k)%Position(:,p%nr) - p0) ) - enddo ! --- Ct and Cq on polar grid (goes beyond rotor radius) if (EqualRealNos(y%DiskAvg_Vx_Rel,0.0_ReKi)) then @@ -658,14 +651,14 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg) do k=1,size(m%ADRotorDisk) ! loop on blades force contribution num = num + dot_product(y%xHat_Disk, m%ADRotorDisk(k)%Moment(:,j) ) end do - y%AzimAvg_Cq(j) = num / (denom * R_rotor) + y%AzimAvg_Cq(j) = 2.0_ReKi*num / (denom * y%D_rotor) end do end if ! --- Variables needed to orient wake planes in "skew" coordinate system ! chi_skew and psi_skew - y%chi_skew = Calc_Chi0(m%Turbine%AD%m%rotors(1)%V_diskAvg, m%turbine%AD%m%rotors(1)%V_dot_x) * R2D ! AeroDyn_IO + y%chi_skew = Calc_Chi0(m%Turbine%AD%m%rotors(1)%V_diskAvg, m%turbine%AD%m%rotors(1)%V_dot_x) ! AeroDyn_IO ! TODO place me in an AeroDyn Function like Calc_Chi0 ! Construct y_hat, orthogonal to x_hat when its z component is neglected (in a projected horizontal plane) @@ -678,8 +671,8 @@ SUBROUTINE FWrap_CalcOutput(p, u, y, m, ErrStat, ErrMsg) zHat_plane(1:3) = zHat_plane/TwoNorm(zHat_plane) zHat_Disk = m%Turbine%AD%Input(1)%rotors(1)%HubMotion%Orientation(3,:,1) ! TODO TODO, shoudn't rotate - tmp_sz_y = dot_product(zHat_Disk,yHat_plane) - tmp_sz_z = -1.0* dot_product(zHat_Disk,zHat_plane) + tmp_sz_y = dot_product(zHat_Disk,yHat_plane) + tmp_sz_z = -1.0_ReKi* dot_product(zHat_Disk,zHat_plane) if ( EqualRealNos(tmp_sz_y,0.0_ReKi) .and. EqualRealNos(tmp_sz_z,0.0_ReKi) ) then y%psi_skew = 0.0_ReKi else diff --git a/modules/awae/src/AWAE.f90 b/modules/awae/src/AWAE.f90 index 3e12c83e59..b4d2bea3d9 100644 --- a/modules/awae/src/AWAE.f90 +++ b/modules/awae/src/AWAE.f90 @@ -210,8 +210,6 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg) real(ReKi) :: xplane_sq, yplane_sq, xysq_Z(3), xzplane_X(3) real(ReKi), ALLOCATABLE :: tmp_xhat_plane(:,:), tmp_yhat_plane(:,:), tmp_zhat_plane(:,:) real(ReKi), ALLOCATABLE :: tmp_Vx_wake(:), tmp_Vz_wake(:), tmp_Vy_wake(:) - real(ReKi) :: Vx_debug - integer(IntKi) :: ILo integer(IntKi) :: maxPln, tmpPln, maxN_wake integer(IntKi) :: i,np1,errStat2 character(*), parameter :: RoutineName = 'LowResGridCalcOutput' @@ -242,7 +240,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg) ! Loop over the entire grid of low resolution ambient wind data to compute: ! 1) the disturbed flow at each point and 2) the averaged disturbed velocity of each wake plane - !$OMP PARALLEL DO PRIVATE(nx_low,ny_low,nz_low, nXYZ_low, n_wake, xhatBar_plane, x_end_plane,nt,np,ILo,x_start_plane,delta,deltad,Vx_debug,& + !$OMP PARALLEL DO PRIVATE(nx_low,ny_low,nz_low, nXYZ_low, n_wake, xhatBar_plane, x_end_plane,nt,np,x_start_plane,delta,deltad,& !$OMP& p_tmp_plane,tmp_vec,r_vec_plane,y_tmp_plane,z_tmp_plane,xhatBar_plane_norm, Vx_wake_tmp, Vr_wake_tmp,& !$OMP& nw,Vr_term,Vx_term,tmp_x,tmp_y,tmp_z,& !$OMP& xHat_plane, yHat_plane, zHat_plane,& @@ -277,8 +275,6 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg) x_end_plane = tmp_x + tmp_y + tmp_z do np = 0, maxPln - ! Reset interpolation counter - ILo = 0 ! TODO, curently not used np1 = np + 1 ! Construct the endcaps of the current wake plane volume x_start_plane = x_end_plane @@ -333,7 +329,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg) ! Increment number of wakes contributing to current grid point n_wake = n_wake + 1 - ! Store unit vectors for projection, TODO Compute Vr_term directly + ! Store unit vectors for projection tmp_xhat_plane(:,n_wake) = xHat_plane tmp_yhat_plane(:,n_wake) = yHat_plane tmp_zhat_plane(:,n_wake) = zHat_plane @@ -427,7 +423,7 @@ subroutine LowResGridCalcOutput(n, u, p, y, m, errStat, errMsg) Vsum_low = 0.0_ReKi iwsum = 0 - n_r_polar = FLOOR((p%C_Meander*u%D_wake(np,nt))/(2.0_ReKi*p%dpol)) ! TODO change me Dwake*sqrt(2)? + n_r_polar = FLOOR((p%C_Meander*u%D_wake(np,nt))/(2.0_ReKi*p%dpol)) do nr = 0,n_r_polar @@ -567,7 +563,7 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg) real(ReKi) :: p_tmp_plane(3) real(ReKi) :: tmp_vec(3) real(ReKi) :: delta, deltad - integer(IntKi) :: np1, ILo + integer(IntKi) :: np1 integer(IntKi) :: maxPln integer(IntKi) :: n_high_low character(*), parameter :: RoutineName = 'HighResGridCalcOutput' @@ -612,8 +608,6 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg) x_end_plane = dot_product(u%xhat_plane(:,0,nt2), (p%Grid_high(:,nXYZ_high,nt) - u%p_plane(:,0,nt2)) ) do np = 0, maxPln !p%NumPlanes-2 - ! Reset interpolation counter - ILo = 0 ! TODO, curently not used np1 = np + 1 ! Construct the endcaps of the current wake plane volume x_start_plane = x_end_plane @@ -663,7 +657,7 @@ subroutine HighResGridCalcOutput(n, u, p, y, m, errStat, errMsg) ! Increment number of wakes contributing to current grid point n_wake = n_wake + 1 - ! Store unit vectors for projection, TODO Compute Vr_term directly + ! Store unit vectors for projection m%xhat_plane(:,n_wake) = xHat_plane m%yhat_plane(:,n_wake) = yHat_plane m%zhat_plane(:,n_wake) = zHat_plane diff --git a/modules/awae/src/AWAE_Registry.txt b/modules/awae/src/AWAE_Registry.txt index 36c82fe870..815db3e2e9 100644 --- a/modules/awae/src/AWAE_Registry.txt +++ b/modules/awae/src/AWAE_Registry.txt @@ -212,8 +212,8 @@ typedef ^ OutputType ReKi Vx_wind_disk {:} - # Define inputs that are contained on the mesh here: typedef ^ InputType ReKi xhat_plane {:}{:}{:} - - "Orientations of wake planes, normal to wake planes, for each turbine" - typedef ^ InputType ReKi p_plane {:}{:}{:} - - "Center positions of wake planes for each turbine" m -typedef ^ InputType ReKi Vx_wake {:}{:}{:}{:} - - "Axial wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s -typedef ^ InputType ReKi Vy_wake {:}{:}{:}{:} - - "Transverse horizonal wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s -typedef ^ InputType ReKi Vz_wake {:}{:}{:}{:} - - "Transverse nominally vertical wake velocity deficit at wake planes, distributed acrros the plane, for each turbine (ny,nz,np,nWT)" m/s +typedef ^ InputType ReKi Vx_wake {:}{:}{:}{:} - - "Axial wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s +typedef ^ InputType ReKi Vy_wake {:}{:}{:}{:} - - "Transverse horizonal wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s +typedef ^ InputType ReKi Vz_wake {:}{:}{:}{:} - - "Transverse nominally vertical wake velocity deficit at wake planes, distributed across the plane, for each turbine (ny,nz,np,nWT)" m/s typedef ^ InputType ReKi D_wake {:}{:} - - "Wake diameters at wake planes for each turbine" m diff --git a/modules/wakedynamics/src/WakeDynamics.f90 b/modules/wakedynamics/src/WakeDynamics.f90 index 2464a55983..76823b6107 100644 --- a/modules/wakedynamics/src/WakeDynamics.f90 +++ b/modules/wakedynamics/src/WakeDynamics.f90 @@ -1103,25 +1103,6 @@ subroutine WD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, errStat, errMsg ) call Axisymmetric2Cartesian(y%Vx_wake(:,i), y%Vr_wake(:,i), p%r, p%y, p%z, y%Vx_wake2(:,:,i), y%Vy_wake2(:,:,i), y%Vz_wake2(:,:,i)) enddo - ! DEBUG check - do i = 0, min(n+1,p%NumPlanes-1) - if( any(abs(y%Vx_wake2(0:,0,i)-y%Vx_wake(:,i))>1e-20)) then - print*,'Problem Vx',i - print*,'y%Vx_wake2(0:,0,i)',y%Vx_wake2(0:,0,i) - print*,'y%Vx_wake(:,i)',y%Vx_wake(:,i) - STOP - endif - !y%Vy_wake2(0:,0,:) - if (any(abs(y%Vy_wake2(0:,0,i)-y%Vr_wake(:,i))>1e-9)) then - print*,'Problem Vr',i - print*,'y%Vx_wake2(0:,0,i)',y%Vy_wake2(0:,0,i) - print*,'y%Vx_wake(:,i)',y%Vr_wake(:,i) - STOP - endif - enddo - - - end subroutine WD_CalcOutput @@ -1211,8 +1192,9 @@ subroutine InitStatesWithInputs(numPlanes, numRadii, u, p, xd, m, errStat, errMs xd%Vx_rel_disk_filt = u%Vx_rel_disk - ! Initialze Ct_azavg_filt and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here. - xd%Ct_azavg_filt (:) = u%Ct_azavg(:) + ! Initialze Ct_azavg_filt, Cq_azavg_filt, and Vx_wake; Vr_wake is already initialized to zero, so, we don't need to do that here. + xd%Ct_azavg_filt (:) = u%Ct_azavg(:) + xd%Cq_azavg_filt (:) = u%Cq_azavg(:) call NearWakeCorrection( xd%Ct_azavg_filt, xd%Cq_azavg_filt, xd%Vx_rel_disk_filt, p, m, xd%Vx_wake(:,0), xd%D_rotor_filt(0), errStat, errMsg ) xd%Vx_wake(:,1) = xd%Vx_wake(:,0)