Added Mass and Buoyancy Inputs

src/DMS.jl

@@ -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)
@@ -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
 
@@ -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
@@ -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?
@@ -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
@@ -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)
@@ -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
 
@@ -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

src/OWENSAero.jl

@@ -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
@@ -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)
@@ -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
@@ -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)
@@ -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")

src/Unsteady_Step.jl

@@ -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
@@ -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