Merge branch 'master' of github.com:sandialabs/OWENS.jl into dev

docs/src/literate/C_customizablePreprocessing.jl

@@ -22,7 +22,7 @@ import PyPlot
 import OWENSOpenFASTWrappers
 
 
-path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #bld_height_numadsplitdir(@__FILE__)[1]
+path = runpath = "/home/runner/work/OWENS.jl/OWENS.jl/docs/src/literate" #splitdir(@__FILE__)[1]
 
 Inp = OWENS.MasterInput("$runpath/sampleOWENS.yml")
 

examples/dev/RM2.jl

@@ -104,7 +104,7 @@ controlpts = [3.6479257474344826, 6.226656883619295, 9.082267631309085, 11.44933
 z_shape1 = collect(LinRange(0,41.9,length(controlpts)+2))
 x_shape1 = [0.0;controlpts;0.0]
 z_shape = collect(LinRange(0,41.9,60))
-x_shape = safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
+x_shape = OWENS.safeakima(z_shape1,x_shape1,z_shape)#[0.0,1.7760245854312287, 5.597183088188207, 8.807794161662574, 11.329376903432605, 13.359580331518579, 14.833606099357858, 15.945156349709, 16.679839160110422, 17.06449826588358, 17.10416552269884, 16.760632435904647, 16.05982913536134, 15.02659565585254, 13.660910465851046, 11.913532434360155, 9.832615229216344, 7.421713825584581, 4.447602800040282, 0.0]
 toweroffset = 4.3953443986241725
 SNL34_unit_xz = [x_shape;;z_shape]
 SNL34x = SNL34_unit_xz[:,1]./maximum(SNL34_unit_xz[:,1])
@@ -131,6 +131,7 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng,
     R = Blade_Radius, #windio
     shapeZ, #windio
     shapeX, #windio
+    shapeY=zero(shapeX), #windio
     ifw, #modeling options
     delta_t, #modeling options
     numTS, #modeling options
@@ -156,6 +157,9 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng,
     AModel, #AD, DMS, AC
     DSModel="BV",
     RPI=true,
+    AM_flag = true,
+    rotAccel_flag = true,
+    buoy_flag = true,
     cables_connected_to_blade_base = true,
     angularOffset = pi/2,
     meshtype = turbineType)

src/SetupTurbine.jl

@@ -48,6 +48,9 @@ function setupOWENS(OWENSAero,path;
     AModel="DMS",
     DSModel="BV",
     RPI=true,
+    AM_flag = false,
+    rotAccel_flag = false,
+    buoy_flag = false,
     cables_connected_to_blade_base = true,
     meshtype = "Darrieus",
     custommesh = nothing) #Darrieus, H-VAWT, ARCUS
@@ -435,6 +438,8 @@ function setupOWENS(OWENSAero,path;
         chord = safeakima(bld_height_numad_unit, numadIn_bld.chord,LinRange(bld_height_numad_unit[1],1,Nslices))
         airfoils = fill("nothing",Nslices)
 
+        twist = safeakima(bld_height_numad_unit, numadIn_bld.twist_d.*pi/180,LinRange(bld_height_numad_unit[1],1,Nslices))
+
         # Discretely assign the airfoils
         for (iheight_numad,height_numad) in enumerate(bld_height_numad_unit)
             for (iheight,height_slices) in enumerate(collect(LinRange(0,1,Nslices)))
@@ -445,15 +450,19 @@ function setupOWENS(OWENSAero,path;
         end
 
         OWENSAero.setupTurb(shapeX,shapeZ,B,chord,tsr,Vinf;AModel,DSModel,
-        afname = airfoils, #TODO: map to the numad input
+        afname = airfoils,
         bld_y = shapeY,
         rho,
+        twist, #TODO: verify twist is in same direction
         mu,
         eta,
         ifw, #TODO: propogate WindType
         turbsim_filename = windINPfilename,
         ifw_libfile,
         tau = [1e-5,1e-5],
+        AM_flag,
+        rotAccel_flag,
+        buoy_flag,
         ntheta,
         Nslices,
         RPI)

src/Unsteady.jl

@@ -384,7 +384,7 @@ function Unsteady(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
                             if inputs.AD15On
                                 aeroVals,aeroDOFs = run_aero_with_deformAD15(aero,deformAero,topMesh,topEl,u_j,udot_j,uddot_j,inputs,t[i],azi_j,Omega_j,OmegaDot_j)
                             else
-                                aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,u_j,inputs,numIterations,t[i],azi_j,Omega_j)
+                                aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,u_j,uddot_j,inputs,numIterations,t[i],azi_j,Omega_j,topModel.gravityOn)
                             end
                         end
                     end
@@ -577,7 +577,7 @@ function Unsteady(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
                                 if inputs.AD15On
                                     aeroVals,aeroDOFs = run_aero_with_deformAD15(aero,deformAero,topMesh,topEl,u_j,udot_j,uddot_j,inputs,t[i],azi_j,Omega_j,OmegaDot_j)
                                 else
-                                    aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,u_j,inputs,numIterations,t[i],azi_j,Omega_j)
+                                    aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,u_j,uddot_j,inputs,numIterations,t[i],azi_j,Omega_j,topModel.gravityOn)
                                 end
                             end
                         end
@@ -933,7 +933,7 @@ function structuralDynamicsTransientGX(topModel,mesh,Fexternal,ForceDof,system,a
     return (strainGX,curvGX), dispOut, FReaction_j,systemout
 end
 
-function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t_i,azi_j,Omega_j)
+function run_aero_with_deform(aero,deformAero,mesh,el,u_j,uddot_j,inputs,numIterations,t_i,azi_j,Omega_j,gravityOn)
 
     if inputs.tocp_Vinf == -1
         newVinf = -1
@@ -946,6 +946,7 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
         numDofPerNode = 6
 
         u_j_local = zeros(Int(max(maximum(mesh.structuralNodeNumbers))*6))
+        uddot_j_local = zeros(Int(max(maximum(mesh.structuralNodeNumbers))*6))
         for jbld = 1:length(mesh.structuralElNumbers[:,1])
             for kel = 1:length(mesh.structuralElNumbers[1,:])-1
                 # orientation angle,xloc,sectionProps,element order]
@@ -956,6 +957,7 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
                 dofList = [(node1-1)*numDofPerNode.+(1:6);(node2-1)*numDofPerNode.+(1:6)]
 
                 globaldisp = u_j[dofList]
+                globaldispddot = uddot_j[dofList]
 
                 twist = el.props[elNum].twist
                 sweepAngle = el.psi[elNum]
@@ -965,10 +967,12 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
                 twistAvg = rollAngle + 0.5*(twist[1] + twist[2])
                 lambda = OWENSFEA.calculateLambda(sweepAngle*pi/180.0,coneAngle*pi/180.0,twistAvg.*pi/180.0)
                 localdisp = lambda*globaldisp
+                localdispddot = lambda*globaldispddot
 
                 #asssembly
                 for m = 1:length(dofList)
                     u_j_local[dofList[m]] =  u_j_local[dofList[m]]+localdisp[m]
+                    uddot_j_local[dofList[m]] =  uddot_j_local[dofList[m]]+localdispddot[m]
                 end
 
             end
@@ -978,13 +982,17 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
         disp_y = [u_j[i] for i = 2:6:length(u_j)]
         disp_z = [u_j[i] for i = 3:6:length(u_j)]
         disp_twist = [u_j_local[i] for i = 4:6:length(u_j_local)]
+        dispddot_flap = [uddot_j_local[i] for i = 3:6:length(u_j_local)]
+        dispddot_edge = [uddot_j_local[i] for i = 2:6:length(u_j_local)] #TODO: verify these are correct
         # disp_twist2 = [u_j_local[i] for i = 5:6:length(u_j_local)]
         # disp_twist3 = [u_j_local[i] for i = 6:6:length(u_j_local)]
 
         bld_x = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
         bld_y = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
         bld_z = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
         bld_twist = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
+        accel_flap_in = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
+        accel_edge_in = copy(mesh.structuralNodeNumbers[:,1:end-1]).*0.0
 
         for jbld = 1:length(mesh.structuralNodeNumbers[:,1])
             bld_indices = Int.(mesh.structuralNodeNumbers[jbld,1:end-1])
@@ -994,6 +1002,8 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
             # flatten blade x,y
             bld_x[jbld,:] = sqrt.(bld_x[jbld,:].^2 .+bld_y[jbld,:].^2) #TODO: a better way via the blade offset azimuth?
             bld_twist[jbld,:] = -disp_twist[bld_indices] #the bending displacements are in radians
+            accel_flap_in[jbld,:] = dispddot_flap[bld_indices] # TODO: verify sign
+            accel_edge_in[jbld,:] = dispddot_edge[bld_indices] # TODO: verify sign
             # The local structural FOR follows right hand rule, so at the bottom of the blade, the x-vector is pointing outward, and a positive
             # rotation about x would make the blade twist into the turbine.  In AC and DMS, if we are looking at say Andrew's 2016 paper, fig 10,
             # the blade has looped up and is pointing at us, so positive twist would increase the aoa.  In the AC and DMS equations, aoa is decreased by twist
@@ -1011,10 +1021,31 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t
         bld_x = -1
         bld_z = -1
         bld_twist = -1
+        accel_flap_in = -1
+        accel_edge_in = -1
+    end
+
+    # TODO: gravity vector changing with platform motion
+    if eltype(gravityOn) == Bool && gravityOn == true
+        gravity = [0.0,0.0,-9.81]
+    elseif eltype(gravityOn) == Bool && gravityOn == false
+        gravity = [0.0,0.0,0.0]
+    end
+
+    if eltype(gravityOn) == Float64
+        gravity = gravityOn
     end
 
     # println("Calling Aero $(Omega_j*60) RPM $newVinf Vinf")
-    deformAero(azi_j;newOmega=Omega_j*2*pi,newVinf,bld_x,bld_z,bld_twist) #TODO: implement deformation induced velocities
+    deformAero(azi_j;newOmega=Omega_j*2*pi,
+    newVinf,
+    bld_x,
+    bld_z,
+    bld_twist,
+    accel_flap_in,
+    accel_edge_in,
+    gravity) #TODO: implement deformation induced velocities
+
     aeroVals,aeroDOFs = aero(t_i,azi_j)
     # println(maximum(abs.(aeroVals)))
     return aeroVals,aeroDOFs

src/Unsteady_Land.jl

@@ -400,7 +400,7 @@ function Unsteady_Land(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
                             aeroVals = aeroVals[1]
                             aeroDOFs = aeroDOFs[1]
                         else
-                            aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,topdata.u_j,inputs,numIterations,t[i],topdata.azi_j,topdata.Omega_j)
+                            aeroVals,aeroDOFs = run_aero_with_deform(aero,deformAero,topMesh,topEl,topdata.u_j,topdata.uddot_j,inputs,numIterations,t[i],topdata.azi_j,topdata.Omega_j,topModel.gravityOn)
                         end
                     end
                 end
@@ -421,6 +421,7 @@ function Unsteady_Land(inputs;topModel=nothing,topMesh=nothing,topEl=nothing,
                     if inputs.AD15On
                         # AD15 is in global frame, so no frame conversion???
                         topdata.topFexternal = aeroVals
+                        full_aeroDOFs = aeroDOFs
                     else
                         if length(size(aeroVals))==1 || size(aeroVals)[2]==1 #i.e. the standard aero force input as a long array
                             # Fill in forces and dofs if they were specified not in full arrays TODO: make this more efficient

src/utilities.jl

@@ -330,7 +330,7 @@ function owens(owensfile,analysisType;
 
         aeroForces(t,azi) = externalForcing(t+delta_t,aerotimeArray,aeroForceValHist,aeroForceDof)
 
-        deformAero(azi;newOmega=-1,newVinf=-1,bld_x=-1,bld_z=-1,bld_twist=-1) = 0.0 #placeholder function
+        deformAero(azi;newOmega=-1,newVinf=-1,bld_x=-1,bld_z=-1,bld_twist=-1,accel_flap_in=-1,accel_edge_in=-1,gravity=-1) = 0.0 #placeholder function
         Unsteady(model;topModel=feamodel,topMesh=mesh,topEl=el,bin,aero=aeroForces,deformAero)
 
         return model

test/34mROMtest.jl

@@ -33,7 +33,7 @@ for i=1:length(myel2["props"])
 end
 myel = OWENSFEA.El(sectionPropsArray,myel2["elLen"],myel2["psi"],myel2["theta"],myel2["roll"],myel2["rotationalEffects"])
 
-deformTurb2(a;newOmega=-1,newVinf=-1, bld_x=-1, bld_z=-1, bld_twist=-1) = 0
+deformTurb2(a;newOmega=-1,newVinf=-1, bld_x=-1, bld_z=-1, bld_twist=-1,accel_flap_in=-1,accel_edge_in=-1,gravity=-1) = 0
 
 delta_t = 0.1
 function aeroForcesDMS(t,azi)