From da0a16bdb784b3bbad2882d157b542fd1085d866 Mon Sep 17 00:00:00 2001 From: kevmoor Date: Thu, 22 Aug 2024 10:53:32 -0600 Subject: [PATCH 1/4] Propogate acceleration inputs and flags for added mass --- examples/dev/RM2.jl | 5 ++++- src/SetupTurbine.jl | 9 ++++++++- src/Unsteady.jl | 27 +++++++++++++++++++++++---- src/Unsteady_Land.jl | 2 +- 4 files changed, 36 insertions(+), 7 deletions(-) diff --git a/examples/dev/RM2.jl b/examples/dev/RM2.jl index 4ec7f9cd..aba30e6d 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,8 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng, AModel, #AD, DMS, AC DSModel="BV", RPI=true, + AM_flag = true, + rotAccel_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..72dfc7d9 100644 --- a/src/SetupTurbine.jl +++ b/src/SetupTurbine.jl @@ -48,6 +48,8 @@ function setupOWENS(OWENSAero,path; AModel="DMS", DSModel="BV", RPI=true, + AM_flag = false, + rotAccel_flag = false, cables_connected_to_blade_base = true, meshtype = "Darrieus", custommesh = nothing) #Darrieus, H-VAWT, ARCUS @@ -435,6 +437,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 +449,18 @@ 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, ntheta, Nslices, RPI) diff --git a/src/Unsteady.jl b/src/Unsteady.jl index 48feb862..6e051c40 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) 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) 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) 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,19 @@ 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 # 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) #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..34530219 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) end end end From 368c5392c19c0a8adfd379749c11ce78b90a5343 Mon Sep 17 00:00:00 2001 From: kevmoor Date: Thu, 22 Aug 2024 14:17:15 -0600 Subject: [PATCH 2/4] Buoyancy runs --- examples/dev/RM2.jl | 1 + src/SetupTurbine.jl | 2 ++ src/Unsteady.jl | 20 ++++++++++++++++---- src/Unsteady_Land.jl | 2 +- 4 files changed, 20 insertions(+), 5 deletions(-) diff --git a/examples/dev/RM2.jl b/examples/dev/RM2.jl index aba30e6d..0f3c05ae 100644 --- a/examples/dev/RM2.jl +++ b/examples/dev/RM2.jl @@ -159,6 +159,7 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections,AD15bldNdIdxRng, 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 72dfc7d9..a8b4458f 100644 --- a/src/SetupTurbine.jl +++ b/src/SetupTurbine.jl @@ -50,6 +50,7 @@ function setupOWENS(OWENSAero,path; 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 @@ -461,6 +462,7 @@ function setupOWENS(OWENSAero,path; 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 6e051c40..7093cbce 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,uddot_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,top_grav_setting) 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,uddot_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,top_grav_setting) 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,uddot_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 @@ -1025,6 +1025,17 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,uddot_j,inputs,numIter 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, @@ -1032,7 +1043,8 @@ function run_aero_with_deform(aero,deformAero,mesh,el,u_j,uddot_j,inputs,numIter bld_z, bld_twist, accel_flap_in, - accel_edge_in) #TODO: implement deformation induced velocities + accel_edge_in, + gravity) #TODO: implement deformation induced velocities aeroVals,aeroDOFs = aero(t_i,azi_j) # println(maximum(abs.(aeroVals))) diff --git a/src/Unsteady_Land.jl b/src/Unsteady_Land.jl index 34530219..46a04d14 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,topdata.uddot_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 From c17f20bb071110756e7fcc6c0f70e102d7d37fce Mon Sep 17 00:00:00 2001 From: kevmoor Date: Thu, 22 Aug 2024 14:30:24 -0600 Subject: [PATCH 3/4] AD interface still runs --- src/Unsteady_Land.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Unsteady_Land.jl b/src/Unsteady_Land.jl index 46a04d14..c227ad7e 100644 --- a/src/Unsteady_Land.jl +++ b/src/Unsteady_Land.jl @@ -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 From b2c735dcc1e1959841b834000002881c3670f6ed Mon Sep 17 00:00:00 2001 From: kevmoor Date: Fri, 30 Aug 2024 09:24:29 -0600 Subject: [PATCH 4/4] Update items caught by CICD --- docs/src/literate/C_customizablePreprocessing.jl | 2 +- src/Unsteady.jl | 4 ++-- src/utilities.jl | 2 +- test/34mROMtest.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) 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/src/Unsteady.jl b/src/Unsteady.jl index 7093cbce..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,uddot_j,inputs,numIterations,t[i],azi_j,Omega_j,top_grav_setting) + 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,uddot_j,inputs,numIterations,t[i],azi_j,Omega_j,top_grav_setting) + 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 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)