diff --git a/docs/src/literate/C_customizablePreprocessing.jl b/docs/src/literate/C_customizablePreprocessing.jl index 05aa24fa..acfee8aa 100644 --- a/docs/src/literate/C_customizablePreprocessing.jl +++ b/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") diff --git a/examples/dev/RM2.jl b/examples/dev/RM2.jl index 4ec7f9cd..0f3c05ae 100644 --- a/examples/dev/RM2.jl +++ b/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) diff --git a/src/SetupTurbine.jl b/src/SetupTurbine.jl index ee005f44..a8b4458f 100644 --- a/src/SetupTurbine.jl +++ b/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) diff --git a/src/Unsteady.jl b/src/Unsteady.jl index 48feb862..8d6b9ace 100644 --- a/src/Unsteady.jl +++ b/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,6 +982,8 @@ 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)] @@ -985,6 +991,8 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,inputs,numIterations,t 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 diff --git a/src/Unsteady_Land.jl b/src/Unsteady_Land.jl index 9d286e59..c227ad7e 100644 --- a/src/Unsteady_Land.jl +++ b/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 diff --git a/src/utilities.jl b/src/utilities.jl index b514b05d..8402bf26 100644 --- a/src/utilities.jl +++ b/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 diff --git a/test/34mROMtest.jl b/test/34mROMtest.jl index 0240af6c..570b0eae 100644 --- a/test/34mROMtest.jl +++ b/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)