Skip to content

Commit

Permalink
Merge pull request #31 from sandialabs/added_mass
Browse files Browse the repository at this point in the history 
Added Mass and Buoyancy Inputs
  • Loading branch information
@kevmoor
kevmoor authored Aug 29, 2024
2 parents 7a2f320 + 040c8dd commit 0af8257
Show file tree
Hide file tree
Showing 12 changed files with 239 additions and 71 deletions.
74 changes: 63 additions & 11 deletions src/DMS.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,16 +22,21 @@ end
"""
function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solvestep=false)

# Unpack Vars
B = turbine.B
k = 1.0
af = turbine.af
chord = turbine.chord[1]
thickness = turbine.thick[1]
ntheta = turbine.ntheta
suction = false
rho = env.rho
mu = env.mu
gravity = env.gravity
AM_flag = env.AM_flag
buoy_flag = env.buoy_flag
rotAccel_flag = env.rotAccel_flag
# winddir = env.winddir

dtheta = 2*pi/(ntheta) #Assuming discretization is fixed equidistant (but omega can change between each point)
Expand All @@ -47,6 +52,8 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
V_yt = env.V_y[idx]
V_zt = env.V_z[idx]
V_twist = env.V_twist[idx]
accel_flap = env.accel_flap[idx]
accel_edge = env.accel_edge[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 @@ -61,6 +68,8 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
V_y = env.V_y[1]
V_z = env.V_z[1]
V_twist = env.V_twist[1]
accel_flap = env.accel_flap[1]
accel_edge = env.accel_edge[1]
# V_delta = env.V_delta[1] # Does not apply since the model calculation is centered around the point of rotation
# V_sweep = env.V_sweep[1] # Does not apply since the model calculation is centered around the point of rotation
end
Expand Down Expand Up @@ -115,14 +124,53 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
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
if AM_flag
Vol_flap = pi * (chord/2)^2 * 1.0
Vol_edge = pi * ((thickness/10)/2)^2 * 1.0

if rotAccel_flag
accel_rot = omega^2 * r
else
accel_rot = 0.0
end

M_addedmass_flap = rho * env.AM_Coeff_Ca * Vol_flap
M_addedmass_edge = rho * env.AM_Coeff_Ca * Vol_edge

F_addedmass_flap = M_addedmass_flap * (accel_flap+accel_rot)
F_addedmass_edge = M_addedmass_edge * (accel_edge+accel_rot)

M_addedmass_Np = M_addedmass_flap*cos(twist) + M_addedmass_edge*sin(twist) # Go from the beam frame of reference to the normal and tangential direction #TODO: verify the directions
M_addedmass_Tp = M_addedmass_edge*cos(twist) - M_addedmass_flap*sin(twist)

# Instantaneous Forces (Unit Height) #Based on this, radial is inward and tangential is in direction of rotation
F_addedmass_Np = F_addedmass_flap*cos(twist) + F_addedmass_edge*sin(twist) # Go from the beam frame of reference to the normal and tangential direction #TODO: verify the directions
F_addedmass_Tp = F_addedmass_edge*cos(twist) - F_addedmass_flap*sin(twist)
else
M_addedmass_Np = 0.0
M_addedmass_Tp = 0.0
F_addedmass_Np = 0.0
F_addedmass_Tp = 0.0
end

if buoy_flag
section_volume = chord*thickness/2*1.0 # per unit length TODO: input volume
dcm = [cos(theta) -sin(theta) 0
sin(theta) cos(theta) 0
0 0 1]
F_buoy = dcm * gravity .* (rho*section_volume) #TODO: verify direction
else
F_buoy = [0.0, 0.0, 0.0]
end

# Instantaneous Forces (Unit Height) #Based on this, radial is inward positive and tangential is in direction of rotation positive
Ab = chord * 1.0 # planform area Assuming unit section height
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

Np = cn.*q_loc + -F_addedmass_Np
Rp = Np + F_buoy[2]# Ning Eq. 27 # Negate to match cactus frame of reference, note that delta cancels out
Zp = Np*tan(delta) + F_buoy[3] # Ning Eq. 27 # Negate to match cactus frame of reference
Tp = -rotation*(ct.*q_loc + -F_addedmass_Tp)/cos(delta) + F_buoy[1] # TODO: verify direction Ning Eq. 27 # Negate to match cactus frame of reference

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?
Expand All @@ -134,7 +182,7 @@ function streamtube(a,theta,turbine,env;output_all=false,Vxwake=nothing,solveste
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, M_addedmass_Np, M_addedmass_Tp, F_addedmass_Np, F_addedmass_Tp
else
return CD-CT # Residual, section 2.4
end
Expand Down Expand Up @@ -185,7 +233,11 @@ function DMS(turbine, env; w=0, idx_RPI=1:turbine.ntheta, solve=true)
cl = zeros(Real,ntheta)
cd_af = zeros(Real,ntheta)
Re = zeros(Real,ntheta)

M_addedmass_Np = zeros(Real,ntheta)
M_addedmass_Tp = zeros(Real,ntheta)
F_addedmass_Np = zeros(Real,ntheta)
F_addedmass_Tp = zeros(Real,ntheta)

# For All Upper
iter_RPI = 1
for i = 1:Int(ntheta/2)
Expand All @@ -197,7 +249,7 @@ function DMS(turbine, env; w=0, idx_RPI=1:turbine.ntheta, solve=true)
astar[i], _ = FLOWMath.brent(resid, 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], M_addedmass_Np[i], M_addedmass_Tp[i], F_addedmass_Np[i], F_addedmass_Tp[i] = streamtube(astar[i],thetavec[i],turbine,env;output_all=true)
end
end

Expand All @@ -221,13 +273,13 @@ function DMS(turbine, env; w=0, idx_RPI=1:turbine.ntheta, solve=true)
astar[i], _ = FLOWMath.brent(resid, 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], M_addedmass_Np[i], M_addedmass_Tp[i], F_addedmass_Np[i], F_addedmass_Tp[i] = streamtube(astar[i],thetavec[i],turbine,env;output_all=true,Vxwake)
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

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, M_addedmass_Np, M_addedmass_Tp, F_addedmass_Np, F_addedmass_Tp
end
32 changes: 24 additions & 8 deletions src/OWENSAero.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@ struct Turbine{TF1,TF2,TI1,TI2,TAF0,TAF1,TAF2,TAF3,TAF4,TAF5,TFN,TB,TAI}
r::TAF1
z::TF2
chord::TAF3
thick::TAF3
twist::TAF5
delta::TAF0
omega::TAF4
Expand All @@ -67,8 +68,8 @@ struct Turbine{TF1,TF2,TI1,TI2,TAF0,TAF1,TAF2,TAF3,TAF4,TAF5,TFN,TB,TAI}
helical_offset::TAI
end

Turbine(R,r,z,chord,twist,delta,omega,B,af,ntheta,r_delta_infl) = Turbine(R,r,z,chord,twist,delta,omega,B,af,ntheta,r_delta_infl,zeros(Real,size(R)),zeros(Real,size(R)),zeros(1))
Turbine(R,r,chord,twist,delta,omega,B,af,ntheta,r_delta_infl) = Turbine(R,r,1.0,chord,twist,delta,omega,B,af,ntheta,r_delta_infl,zeros(Real,size(R)),zeros(Real,size(R)),zeros(1))
Turbine(R,r,z,chord,twist,delta,omega,B,af,ntheta,r_delta_infl) = Turbine(R,r,z,chord,0.18,twist,delta,omega,B,af,ntheta,r_delta_infl,zeros(Real,size(R)),zeros(Real,size(R)),zeros(1))
Turbine(R,r,chord,twist,delta,omega,B,af,ntheta,r_delta_infl) = Turbine(R,r,1.0,chord,0.18,twist,delta,omega,B,af,ntheta,r_delta_infl,zeros(Real,size(R)),zeros(Real,size(R)),zeros(1))

"""
Environment(rho::TF,mu::TF,V_x::TAF #Vinf is Vx,V_y::TAF,V_z::TAF,V_twist::TAF,windangle::TF #radians,DSModel::TS,AModel::TS,aw_warm::TVF,steplast::TAI,idx_RPI::TAI,V_wake_old::TVF2,BV_DynamicFlagL::TAI,BV_DynamicFlagD::TAI,alpha_last::TAF2,suction::TB)
Expand Down Expand Up @@ -100,16 +101,20 @@ Contains specications for turbine slice environment/operating conditions as well
* `none`:
"""
struct Environment{TF,TB,TAF,TAF2,TS,TVF,TVF2,TAI,TAF3}
struct Environment{TF,TB,TAFx,TAFy,TAF2,TS,TVF,TVF2,TAI,TAF3,TAF4,TAF5}
rho::TF
mu::TF
V_x::TAF #Vinf is Vx
V_y::TAF
V_x::TAFx #Vinf is Vx
V_y::TAFy
V_z::TAF3
V_twist::TAF3
windangle::TF #radians
DSModel::TS
AModel::TS
AM_flag::TB
buoy_flag::TB
rotAccel_flag::TB
AM_Coeff_Ca::TF
aw_warm::TVF
steplast::TAI
idx_RPI::TAI
Expand All @@ -118,10 +123,13 @@ struct Environment{TF,TB,TAF,TAF2,TS,TVF,TVF2,TAI,TAF3}
BV_DynamicFlagD::TAI
alpha_last::TAF2
suction::TB
accel_flap::TAF4
accel_edge::TAF4
gravity::TAF5
end

Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,aw_warm) = Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,aw_warm,zeros(Int,1),zeros(Int,length(V_x)),deepcopy(V_x),zeros(Int,length(V_x)),zeros(Int,length(V_x)),zeros(Real,length(V_x)),false)
Environment(rho,mu,V_x,DSModel,AModel,aw_warm) = Environment(rho,mu,V_x,zeros(Real,size(V_x)),zeros(Real,size(V_x)),zeros(Real,size(V_x)),0.0,DSModel,AModel,aw_warm,zeros(Int,1),zeros(Int,length(V_x)),deepcopy(V_x),zeros(Int,length(V_x)),zeros(Int,length(V_x)),zeros(Real,length(V_x)),false)
Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,aw_warm) = Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,false,false,false,1.0,aw_warm,zeros(Int,1),zeros(Int,length(V_x)),deepcopy(V_x),zeros(Int,length(V_x)),zeros(Int,length(V_x)),zeros(Real,length(V_x)),false,zeros(Real,length(V_x)),zeros(Real,length(V_x)),[0.0,0.0,-9.81])
Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,AM_flag,buoy_flag,rotAccel_flag,AM_Coeff_Ca,aw_warm) = Environment(rho,mu,V_x,V_y,V_z,V_twist,windangle,DSModel,AModel,AM_flag,buoy_flag,rotAccel_flag,AM_Coeff_Ca,aw_warm,zeros(Int,1),zeros(Int,length(V_x)),deepcopy(V_x),zeros(Int,length(V_x)),zeros(Int,length(V_x)),zeros(Real,length(V_x)),false,zeros(Real,length(V_x)),zeros(Real,length(V_x)),[0.0,0.0,-9.81])
Environment(rho,mu,V_x,DSModel,AModel,aw_warm) = Environment(rho,mu,V_x,zeros(Real,size(V_x)),zeros(Real,size(V_x)),zeros(Real,size(V_x)),0.0,DSModel,AModel,false,false,false,1.0,aw_warm,zeros(Int,1),zeros(Int,length(V_x)),deepcopy(V_x),zeros(Int,length(V_x)),zeros(Int,length(V_x)),zeros(Real,length(V_x)),false,zeros(Real,length(V_x)),zeros(Real,length(V_x)),[0.0,0.0,-9.81])

"""
UnsteadyParams(RPI::TB,tau::TAF,ifw::TB,IECgust::TB,nominalVinf::TF,G_amp::TF,gustX0::TF,gustT::TF)
Expand Down Expand Up @@ -201,6 +209,14 @@ function steady(turbine, env; w=zeros(Real,2*turbine.ntheta), idx_RPI=1:2*turbin
end
end

@inline function safeakima(x,y,xpt)
if minimum(xpt)<(minimum(x)-abs(minimum(x))*0.1) || maximum(xpt)>(maximum(x)+abs(maximum(x))*0.1)
msg="Extrapolating on akima spline results in undefined solutions minimum(xpt)<minimum(x) $(minimum(xpt))<$(minimum(x)) or maximum(xpt)<maximum(x) $(maximum(xpt))>$(maximum(x))"
throw(OverflowError(msg))
end
return FLOWMath.akima(x,y,xpt)
end

include("DMS.jl")
include("./vawt-ac/src/airfoilread.jl") #TODO: switch for the CCBlade airfoil reading library
include("./vawt-ac/src/acmultiple.jl")
Expand Down
4 changes: 2 additions & 2 deletions src/Unsteady_Step.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ function Unsteady_Step(turbine,env,us_param,mystep)
aw_filtered = aw_warm[:] .* exp.(-[dt;dt] ./ tau_near) .+ awnew[:] .* (1.0 .- exp.(-[dt;dt] ./ tau_near))

# Get the actual performance using the filtered induction factors
CP, Th, Q_temp, Rp_temp, Tp_temp, Zp_temp, Vloc, CD, CT, a, awstar, alpha, cl, cd, thetavec, Re = steady(turbine, env; w=aw_filtered,idx_RPI,solve=false,ifw)
CP, Th, Q_temp, Rp_temp, Tp_temp, Zp_temp, Vloc, CD, CT, a, awstar, alpha, cl, cd, thetavec, Re, M_addedmass_Np, M_addedmass_Tp, F_addedmass_Np, F_addedmass_Tp = steady(turbine, env; w=aw_filtered,idx_RPI,solve=false,ifw)

if env.steplast[1] != mystep
env.aw_warm[:] = aw_filtered[:] #Mutate, do not link
Expand All @@ -121,5 +121,5 @@ function Unsteady_Step(turbine,env,us_param,mystep)
# env.idx_RPI[:] = idx_RPI_half #Mutate, do not link, TODO: fix this, is used for indexing outputs
theta_temp = thetavec[idx_RPI_half[1]]

return CP,Th,Q_temp, Rp_temp, Tp_temp, Zp_temp, Vloc, CD, CT, a, awstar, alpha, cl, cd, thetavec, Re
return CP,Th,Q_temp, Rp_temp, Tp_temp, Zp_temp, Vloc, CD, CT, a, awstar, alpha, cl, cd, thetavec, Re, M_addedmass_Np, M_addedmass_Tp, F_addedmass_Np, F_addedmass_Tp
end
Loading

0 comments on commit 0af8257

Please sign in to comment.