Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Amass buoy #25

Closed
wants to merge 5 commits into from
Closed

Amass buoy #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
71bb8b2
added mass and buoyancy options for underwater turbine
ckelley2 Apr 22, 2024
da29fca
delete DMS_ck, features in DMS
Apr 22, 2024
ffe12a0
added **/.DS_Store to .gitignore
ckelley2 Apr 29, 2024
2579430
added **/.DS_Store to .gitignore
ckelley2 Apr 29, 2024
0b8af44
changes to DMS and advanceTurbine to use blade accelerations for adde…
ckelley2 May 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,3 +15,4 @@ theta*
*.docx
*.csv
*Manual*
**/.DS_Store
158 changes: 94 additions & 64 deletions src/DMS.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

import FLOWMath
import Statistics
import Statistics:mean
import Statistics: mean

"""
INTERNAL streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solvestep=false)
Expand All @@ -21,8 +21,8 @@ else
end

"""
function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solvestep=false)

function streamtube(a, theta, turbine, env; output_all=false, Vxwake=nothing, solvestep=false, add_mass=false, buoy=false)
# Unpack Vars
B = turbine.B
k = 1.0
Expand All @@ -32,12 +32,14 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
suction = false
rho = env.rho
mu = env.mu
# uddot = xx.x[1]

# winddir = env.winddir

dtheta = 2*pi/(ntheta) #Assuming discretization is fixed equidistant (but omega can change between each point)
dtheta = 2 * pi / (ntheta) #Assuming discretization is fixed equidistant (but omega can change between each point)

if length(turbine.r)>1
idx = round(Int, (theta+dtheta/2)/dtheta)
if length(turbine.r) > 1
idx = round(Int, (theta + dtheta / 2) / dtheta)
r = turbine.r[idx]
twist = turbine.twist[idx]
delta = turbine.delta[idx]
Expand All @@ -46,6 +48,10 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
V_xt = env.V_x[idx] # Vinf is V_xt, t is for turbine direction f.o.r.
V_yt = env.V_y[idx]
V_zt = env.V_z[idx]

# CK. Modify to include blade acceleration from structural model
# A_r = env.A_r[idx] # radial acceleration (unsteady, not including centripetal) from structural model. Outward acceleration should have inward, resiting, apparent mass force

V_twist = env.V_twist[idx]
# V_delta = env.V_delta[idx] # Does not apply since the model calculation is centered around the point of rotation
# V_sweep = env.V_sweep[idx] # Does not apply since the model calculation is centered around the point of rotation
Expand All @@ -65,34 +71,34 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
# V_sweep = env.V_sweep[1] # Does not apply since the model calculation is centered around the point of rotation
end

if Vxwake!=nothing
if Vxwake != nothing
V_xt = Vxwake
end

if suction
# Drag Coefficient: Eq. 5
if a<=0.7
CD = 4/3 * a * (3-a)/(1+a)
elseif a > 0.7 && a <=1
CD = 0.889 - (0.0203-(a-0.143)^2)/0.6427
if a <= 0.7
CD = 4 / 3 * a * (3 - a) / (1 + a)
elseif a > 0.7 && a <= 1
CD = 0.889 - (0.0203 - (a - 0.143)^2) / 0.6427
else
error("a > 1 is undefined")
end
else
# Drag Coefficient: Eq. 2
if a<=0.4
CD = 4*a*(1-a)
if a <= 0.4
CD = 4 * a * (1 - a)
else
CD = 0.889 - (0.0203-(a-0.143)^2)/0.6427
CD = 0.889 - (0.0203 - (a - 0.143)^2) / 0.6427
end
end

# Cycle-average Thrust Coefficient
# TSR = abs(omega)*r/sqrt(V_x^2+V_y^2) # Tip Speed Ratio, section 2.3, however since this is being used for relative local velocities, we use the local radius r, and not the nominal radius R
# Vt = V_x*(1-a)*cos(theta)+V_y*sin(theta)+V_y*(a)*cos(theta)+TSR

Vn = (V_xt*(1-a)*sin(theta)-V_yt*cos(theta)+V_yt*(a)*sin(theta))*cos.(delta) + V_zt*sin(delta)
Vt = V_xt*((1-a)*cos(theta))+V_yt*sin(theta)+V_yt*(1-a)*cos(theta)+abs(omega)*r
Vn = (V_xt * (1 - a) * sin(theta) - V_yt * cos(theta) + V_yt * (a) * sin(theta)) * cos.(delta) + V_zt * sin(delta)
Vt = V_xt * ((1 - a) * cos(theta)) + V_yt * sin(theta) + V_yt * (1 - a) * cos(theta) + abs(omega) * r
Vloc = sqrt(Vn^2 + Vt^2)
phi = atan(Vn, Vt)
alpha = phi - twist
Expand All @@ -101,42 +107,65 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
# phi = atan(((1-a)*sin(theta)).*cos.(delta) , ((1-a)*cos(theta)+TSR)) # Eq. 7 but correct for slope from the 3D to 2D frame of reverence, also must use atan2 to determine which quadrant and get the signs correct
# Vloc = Vinf .* ((1-a)^2 + 2*(1-a)*TSR*cos(theta).*cos.(delta) + TSR^2).^0.5 # Eq. 8

Re = rho*Vloc*chord/mu
dt = dtheta/abs.(omega)
Re = rho * Vloc * chord / mu
dt = dtheta / abs.(omega)
v_sound = 343.0 #m/s #TODO: calculate this using Atmosphere.jl
mach = Vloc/v_sound
mach = Vloc / v_sound
if env.DSModel == "BV"
cl, cd_af = af(alpha,Re,mach,env,V_twist,chord,dt,Vloc;solvestep,idx)
cl, cd_af = af(alpha, Re, mach, env, V_twist, chord, dt, Vloc; solvestep, idx)
elseif env.DSModel == "LB"
error("LB Dynamic Stall Model Not Implemented Yet")
else
cl, cd_af = af(alpha,Re,mach)
cl, cd_af = af(alpha, Re, mach)
end

# CK, hard coded parameter for added mass, temporary
# r_ddot = uddot-omega^2*r
r_ddot = omega^2*r
Ca = 1

ct = (cd_af * cos(phi) - cl * sin(phi)) * 1 # Eq. 9
if add_mass
cn = cd_af * sin(phi) + cl * cos(phi) - r_ddot * Ca * pi / 2 * chord / Vloc^2 # modified cn for added mass acceleration
else
cn = cd_af * sin(phi) + cl * cos(phi) # Eq. 10
end
ct = cd_af*cos(phi) - cl*sin(phi) # Eq. 9
cn = cd_af*sin(phi) + cl*cos(phi) # Eq. 10

Ab = chord * 1.0 # Assuming unit section height

# Instantaneous Forces (Unit Height) #Based on this, radial is inward and tangential is in direction of rotation
q_loc = 0.5 * rho * Ab * Vloc^2 # From Eq. 11
Rp = cn.*q_loc # Ning Eq. 27 # Negate to match cactus frame of reference
Tp = -rotation*ct.*q_loc/cos(delta) # Ning Eq. 27 # Negate to match cactus frame of reference
Zp = cn.*q_loc*tan(delta) # Ning Eq. 27 # Negate to match cactus frame of reference

# CK, buoyancy parameters
airfoil_tc = 0.21
airfoil_area = airfoil_tc * chord^2 / 1.46 # for all NACA airfoils, Area = t/c/1.46
airfoil_volume = airfoil_area # assuming per unit length force
F_z_buoy = airfoil_volume * 1000 * 9.81
# CK Need to project buoyancy force into all other force directions. assume gravity is in z direction for now, and turbine is an underwater vawt
if buoy
Rp = cn .* q_loc + 0 # Ning Eq. 27 # Negate to match cactus frame of reference
Tp = -rotation * ct .* q_loc / cos(delta) + 0 # Ning Eq. 27 # Negate to match cactus frame of reference
Zp = cn .* q_loc * tan(delta) + F_z_buoy # Ning Eq. 27 # Negate to match cactus frame of reference
else
Rp = cn .* q_loc # Ning Eq. 27 # Negate to match cactus frame of reference
Tp = -rotation * ct .* q_loc / cos(delta) # Ning Eq. 27 # Negate to match cactus frame of reference
Zp = cn .* q_loc * tan(delta) # Ning Eq. 27 # Negate to match cactus frame of reference
end

Th = q_loc * (ct*cos(theta) + cn*sin(theta)/cos(delta)) # Eq. 11 but with delta correction
Th = q_loc * (ct * cos(theta) + cn * sin(theta) / cos(delta)) # Eq. 11 but with delta correction

Q = q_loc * r * -ct # Eq. 12 but with Local radius for local torque, Negate the force for reaction torque, in the power frame of reference?

Ast = 1.0 * r * dtheta *abs(sin(theta)) # Section 2.0, local radius for local area, however do not allow negative area, so use absolute value
Ast = 1.0 * r * dtheta * abs(sin(theta)) # Section 2.0, local radius for local area, however do not allow negative area, so use absolute value

q_inf = 0.5 * rho * Ast * (V_xt^2+V_yt^2) # this is in the radial frame of reference V_x.^2 , also (sqrt(V_x^2+V_y^2)).^2 is reduced
q_inf = 0.5 * rho * Ast * (V_xt^2 + V_yt^2) # this is in the radial frame of reference V_x.^2 , also (sqrt(V_x^2+V_y^2)).^2 is reduced

CT = (k * B/(2*pi) * dtheta*Th) ./ q_inf # Eq. 13
CT = (k * B / (2 * pi) * dtheta * Th) ./ q_inf # Eq. 13

if output_all
return Th, Q, Rp, Tp, Zp, Vloc, CD, CT, alpha, cl, cd_af, Re
return Th, Q, Rp, Tp, Zp, Vloc, CD, CT, alpha, cl, cd_af, Re, cn
else
return CD-CT # Residual, section 2.4
return CD - CT # Residual, section 2.4
end
end

Expand All @@ -147,87 +176,88 @@ see ?steady for detailed i/o description

Double multiple streamtube model
"""
function DMS(turbine, env; w=0, idx_RPI=1:turbine.ntheta, solve=true)
function DMS(turbine, env; w=0, idx_RPI=1:turbine.ntheta, solve=true, add_mass=false, buoy=false)
#TODO: Inputs documentation
a_in = w
ntheta = turbine.ntheta
#Convert global F.O.R. winds to turbine frame
windangle = env.windangle

V_xtemp = env.V_x.*cos(windangle)+env.V_y.*sin(windangle) # Vinf is V_x, t is for turbine direction f.o.r.
V_ytemp = -env.V_x.*sin(windangle)+env.V_y.*cos(windangle)
V_xtemp = env.V_x .* cos(windangle) + env.V_y .* sin(windangle) # Vinf is V_x, t is for turbine direction f.o.r.
V_ytemp = -env.V_x .* sin(windangle) + env.V_y .* cos(windangle)

env.V_x[:] = V_xtemp #TODO: ensure this doesn't mess up unsteady methods with persistent backend memory
env.V_y[:] = V_ytemp

# Unpack the rest

dtheta = 2*pi/(ntheta)
dtheta = 2 * pi / (ntheta)
# thetavec = collect(dtheta/2:dtheta:2*pi-dtheta/2)
thetavec = collect(dtheta/2:dtheta:2*pi)

astar = zeros(Real,ntheta*2)#env.aw_warm[1:ntheta] #zeros(ntheta)
astar = zeros(Real, ntheta * 2)#env.aw_warm[1:ntheta] #zeros(ntheta)

if a_in == 0
idx_RPI = 1:ntheta
else
astar[1:ntheta] = a_in[1:ntheta]
end
Th = zeros(Real,ntheta)
Q = zeros(Real,ntheta)
Rp = zeros(Real,ntheta)
Tp = zeros(Real,ntheta)
Zp = zeros(Real,ntheta)
Vloc = zeros(Real,ntheta)
CD = zeros(Real,ntheta)
CT = zeros(Real,ntheta)
alpha = zeros(Real,ntheta)
cl = zeros(Real,ntheta)
cd_af = zeros(Real,ntheta)
Re = zeros(Real,ntheta)
Th = zeros(Real, ntheta)
Q = zeros(Real, ntheta)
Rp = zeros(Real, ntheta)
Tp = zeros(Real, ntheta)
Zp = zeros(Real, ntheta)
Vloc = zeros(Real, ntheta)
CD = zeros(Real, ntheta)
CT = zeros(Real, ntheta)
alpha = zeros(Real, ntheta)
cl = zeros(Real, ntheta)
cd_af = zeros(Real, ntheta)
Re = zeros(Real, ntheta)
cn = zeros(Real, ntheta)

# For All Upper
iter_RPI = 1
for i = 1:Int(ntheta/2)
for i = 1:Int(ntheta / 2)
# Solve Residual
if idx_RPI[iter_RPI] == i
iter_RPI+=1
iter_RPI += 1
if solve
resid(a) = streamtube(a,thetavec[i],turbine,env;solvestep=true) #solvestep makes this independent solve not mess up the dynamic stall variables during the solve
astar[i], _ = FLOWMath.brent(resid, 0, .999)
resid(a) = streamtube(a, thetavec[i], turbine, env; solvestep=true) #solvestep makes this independent solve not mess up the dynamic stall variables during the solve
astar[i], _ = FLOWMath.brent(resid, 0, 0.999)
end

Th[i], Q[i], Rp[i], Tp[i], Zp[i], Vloc[i], CD[i], CT[i], alpha[i], cl[i], cd_af[i], Re[i] = streamtube(astar[i],thetavec[i],turbine,env;output_all=true)
Th[i], Q[i], Rp[i], Tp[i], Zp[i], Vloc[i], CD[i], CT[i], alpha[i], cl[i], cd_af[i], Re[i], cn[i] = streamtube(astar[i], thetavec[i], turbine, env; output_all=true, add_mass, buoy)
end
end

Vxsave = mean(env.V_x)

# For All Lower
for i = Int(ntheta/2+1):ntheta
for i = Int(ntheta / 2 + 1):ntheta
# Update Downstream Inflow based on upstream solution
if idx_RPI[iter_RPI] == i
iter_RPI+=1
a_used = astar[Int(ntheta/2)-(i-Int(ntheta/2))+1] # We index in a circle, but the streamtubes go straight through, swapping at 180 deg
iter_RPI += 1
a_used = astar[Int(ntheta / 2)-(i-Int(ntheta / 2))+1] # We index in a circle, but the streamtubes go straight through, swapping at 180 deg
if env.suction
Vxwake = (1.0-a_used)/(1.0+a_used)*env.V_x[i] # Eq. 6
Vxwake = (1.0 - a_used) / (1.0 + a_used) * env.V_x[i] # Eq. 6
else
Vxwake = (1.0-2.0*a_used)*env.V_x[i] # Eq. 3
Vxwake = (1.0 - 2.0 * a_used) * env.V_x[i] # Eq. 3
end

# Solve Residual
if solve
resid(a) = streamtube(a,thetavec[i],turbine,env;Vxwake,solvestep=true) #Solve wrt negative theta on the back side?
astar[i], _ = FLOWMath.brent(resid, 0, .999)
resid(a) = streamtube(a, thetavec[i], turbine, env; Vxwake, solvestep=true) #Solve wrt negative theta on the back side?
astar[i], _ = FLOWMath.brent(resid, 0, 0.999)
end

Th[i], Q[i], Rp[i], Tp[i], Zp[i], Vloc[i], CD[i], CT[i], alpha[i], cl[i], cd_af[i], Re[i] = streamtube(astar[i],thetavec[i],turbine,env;output_all=true,Vxwake)
Th[i], Q[i], Rp[i], Tp[i], Zp[i], Vloc[i], CD[i], CT[i], alpha[i], cl[i], cd_af[i], Re[i], cn[i] = streamtube(astar[i], thetavec[i], turbine, env; output_all=true, Vxwake, add_mass, buoy)
end
end

# Aggregate Turbine Performance
k = 1.0
CP = sum((k * turbine.B/(2*pi) * dtheta*abs.(Q)*mean(abs.(turbine.omega))) / (0.5*env.rho*1.0*2*turbine.R*Vxsave^3)) # Eq. 14, normalized by nominal radius R
CP = sum((k * turbine.B / (2 * pi) * dtheta * abs.(Q) * mean(abs.(turbine.omega))) / (0.5 * env.rho * 1.0 * 2 * turbine.R * Vxsave^3)) # Eq. 14, normalized by nominal radius R

return CP, Th, Q, Rp, Tp, Zp, Vloc, CD, CT, mean(astar[1:ntheta]), astar, alpha, cl, cd_af, thetavec.-windangle, Re
return CP, Th, Q, Rp, Tp, Zp, Vloc, CD, CT, mean(astar[1:ntheta]), astar, alpha, cl, cd_af, thetavec .- windangle, Re, cn
end
18 changes: 17 additions & 1 deletion src/advanceTurbine.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,9 @@ end
function deformTurb(azi;newOmega=-1,newVinf=-1,bld_x=-1,
bld_z=-1,
bld_twist=-1,
bld_x_accel=0,
bld_y_accel=0,
bld_z_accel=0,
steady=false) # each of these is size ntheta x nslices


Expand All @@ -224,13 +227,22 @@ steady=false) # each of these is size ntheta x nslices
if bld_x!=-1 && bld_z!=-1 && bld_twist!=-1
bld_x_temp = zeros(length(bld_x[:,1]),length(z3D))
bld_twist_temp = zeros(length(bld_x[:,1]),length(z3D))
bld_x_accel_temp = zeros(length(bld_x[:,1]),length(z3D))
bld_y_accel_temp = zeros(length(bld_x[:,1]),length(z3D))
bld_z_accel_temp = zeros(length(bld_x[:,1]),length(z3D))
for ibld = 1:length(bld_x[:,1])
bld_x_temp[ibld,:] = FLOWMath.akima(bld_z[ibld,:],bld_x[ibld,:],z3D.-1.0) #TODO: get rid of this -1.0, which is from inflowwind not liking evaluation at the boundary
bld_twist_temp[ibld,:] = FLOWMath.akima(bld_z[ibld,:],bld_twist[ibld,:],z3D.-1.0)
bld_x_accel_temp[ibld,:] = FLOWMath.akima(bld_z[ibld,:],bld_x_accel[ibld,:],z3D.-1.0)
bld_y_accel_temp[ibld,:] = FLOWMath.akima(bld_z[ibld,:],bld_y_accel[ibld,:],z3D.-1.0)
bld_z_accel_temp[ibld,:] = FLOWMath.akima(bld_z[ibld,:],bld_z_accel[ibld,:],z3D.-1.0)
end
bld_x = bld_x_temp
bld_twist = bld_twist_temp
bld_z = zeros(Real,size(bld_x)) #TODO: a better way to do this.
bld_x_accel = bld_x_accel_temp
bld_y_accel = bld_y_accel_temp
bld_z_accel = bld_z_accel_temp
for ibld = 1:length(bld_x[:,1])
bld_z[ibld,:] = z3D.-1.0
end
Expand Down Expand Up @@ -291,6 +303,9 @@ steady=false) # each of these is size ntheta x nslices
# turbslices[islice].chord = chord[islice]
turbslices[islice].delta[bld_idx] = delta3D[ibld,islice]
turbslices[islice].twist[bld_idx] = startingtwist[islice]+bld_twist[ibld,islice]
turbslices[islice].bld_x_accel[bld_idx] = bld_x_accel[ibld,islice]
turbslices[islice].bld_y_accel[bld_idx] = bld_y_accel[ibld,islice]
turbslices[islice].bld_z_accel[bld_idx] = bld_z_accel[ibld,islice]
end

if newOmega != -1
Expand Down Expand Up @@ -347,6 +362,7 @@ Runs a previously initialized aero model (see ?setupTurb) in the unsteady mode (

"""
function advanceTurb(tnew;ts=2*pi/(turbslices[1].omega[1]*turbslices[1].ntheta),azi=-1.0,verbosity=0,alwaysrecalc=nothing,last_step=nothing)

global timelast
global us_param
global turbslices
Expand Down Expand Up @@ -530,7 +546,7 @@ Runs a previously initialized aero model (see ?setupTurb) in the steady state mo

"""
function steadyTurb(;omega = -1,Vinf = -1)

println("AERO FORCES BEING CALCULATED STEADY")
# global us_param #TODO: add turbulence lookup option for steady?
global turbslices
global envslices
Expand Down
Binary file removed test/.DS_Store
View file
Open in desktop
Binary file not shown.