From 9dfdbdf1a5ff5f2f213689bcc492c9f5aa8b7350 Mon Sep 17 00:00:00 2001 From: kevmoor Date: Wed, 11 Sep 2024 17:17:16 -0600 Subject: [PATCH] Add and break out functions that automatically map OWENS boundary conditions and applied forces to GX and automate the campbell diagram generation process --- .gitignore | 6 +- docs/src/index.md | 14 +- docs/src/installation.md | 3 +- examples/SNL34m/SNL34m_Inputs.yml | 4 +- examples/literate/A_simplyRunningOWENS.jl | 2 +- examples/literate/B_detailedInputs.jl | 8 +- src/OWENS.jl | 1 + src/SetupTurbine.jl | 8 +- src/Unsteady.jl | 70 +--------- src/meshing_utilities.jl | 2 +- test/Fig4_5_campbell2.jl | 159 ++++------------------ 11 files changed, 58 insertions(+), 219 deletions(-) diff --git a/.gitignore b/.gitignore index 6e3740a..de77660 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,8 @@ new* *.out *.h5 *.mat -*.csv \ No newline at end of file +*.csv +*.png +*.gif +*.mp4 +*turbsimfiles/ \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 2c624f4..20faebd 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -9,19 +9,23 @@ OWENS is an ontology, or way of coupling modular aerodynamic, structural, hydrod Here are several examples of OWENS use cases, current and past, including the Sandia 34m research turbine. -![SNL34m.](./assets/SNL34m.png){} +![SNL34m.](./assets/SNL34m.png){#fig:34m +width="50%"} Then here is an example of a helical design. Note that arbitrary numbers of struts can be specified in the automatic meshing functions. You can also write your own generalized mesh using the internal building blocks, but it is not thouroughly documented. -![helical.](./assets/helical.png){} +![helical.](./assets/helical.png){#fig:34m +width="50%"} The generalized meshing was modified to include HAWT concepts, like this bi-wing concept. OWENS is capable of axial flow turbines/HAWTs, but it is not a mature feature, and no where near as developed as OpenFAST (i.e. for regular HAWTs it is recommended to use that software). -![biwing.](./assets/biwing.png){} +![biwing.](./assets/biwing.png){#fig:34m +width="50%"} Then, floating turbines are a possibility, though this feature adds another dimension to the nonlinear time stepping convergance and in turn a fair amount of time. Future work is to make this general interface and functionality an easy to use feature (right now it needs a high level of experience to use). -![arcus.](./assets/arcus.png){} +![arcus.](./assets/arcus.png){#fig:34m +width="50%"} # OWENS under the hood @@ -44,7 +48,7 @@ Turbulent inflow is provided by OWENSOpenFASTWrappers.jl and the inflowwind and Rainflow counting was provided by Rainflow.jl, however, this package became orphained and was pulled into the OWENS code base. -![arcus.](./assets/OWENS_Processes.png){} +drawing ## Installation Please follow the instructions on the setup page diff --git a/docs/src/installation.md b/docs/src/installation.md index fbe9af8..794b814 100644 --- a/docs/src/installation.md +++ b/docs/src/installation.md @@ -115,8 +115,7 @@ Additionally, if you find that your ssh is erroring when you try to install pack PubkeyAcceptedAlgorithms +ssh-ed25519 # Install Optional OpenFAST Dependices -Note that this is optional as it is automatically done by julia in the OWENSOpenFASTWrappers.jl deps/build.jl. For Windows, please follow the OpenFAST Windows instructions on the openfast site for the branch referenced here. - +If your system is already set up such that it is capable of compiling OpenFAST, and you are on mac or linux, then you may skip this and rely on the automatically compiled version that are created when the OWENSOpenFAST libraries are installed by Julia. mkdir coderepos cd coderepos # Install openfast coupled libraries !NOTE!: if you change the location of the compiled libraries, you may need to update the rpath variable, or recompile. diff --git a/examples/SNL34m/SNL34m_Inputs.yml b/examples/SNL34m/SNL34m_Inputs.yml index fcac463..770c0ac 100644 --- a/examples/SNL34m/SNL34m_Inputs.yml +++ b/examples/SNL34m/SNL34m_Inputs.yml @@ -30,12 +30,12 @@ controlParameters: AeroParameters: Nslices: 35 # number of OWENSAero discritizations #TODO: AD parameters ntheta: 30 # number of OWENSAero azimuthal discretizations - AModel: AD # AD, DMS, AC + AModel: DMS # AD, DMS, AC adi_lib: nothing adi_rootname: "/SNL34m" structuralParameters: - structuralModel: ROM #GX, TNB, ROM + structuralModel: GX #GX, TNB, ROM nonlinear: false #TODO: propogate ntelem: 20 #tower elements in each nbelem: 60 #blade elements in each diff --git a/examples/literate/A_simplyRunningOWENS.jl b/examples/literate/A_simplyRunningOWENS.jl index 9c0a0dc..0a11ec0 100644 --- a/examples/literate/A_simplyRunningOWENS.jl +++ b/examples/literate/A_simplyRunningOWENS.jl @@ -22,7 +22,7 @@ import OWENS runpath = path = "/home/runner/work/OWENS.jl/OWENS.jl/examples/literate" # to run locally, change to splitdir(@__FILE__)[1] -##runpath = path = "/Users/kevmoor/Documents/coderepos/OWENS_Toolkit/OWENS.jl/examples/literate"#splitdir(@__FILE__)[1] +## runpath = path = splitdir(@__FILE__)[1] Inp = OWENS.MasterInput("$runpath/sampleOWENS.yml") diff --git a/examples/literate/B_detailedInputs.jl b/examples/literate/B_detailedInputs.jl index f289adf..355d8ae 100644 --- a/examples/literate/B_detailedInputs.jl +++ b/examples/literate/B_detailedInputs.jl @@ -53,8 +53,8 @@ if ifw_libfile == "nothing" end Blade_Height = Inp.Blade_Height Blade_Radius = Inp.Blade_Radius -numTS = Inp.numTS -delta_t = Inp.delta_t +numTS = 1000#Inp.numTS +delta_t = 0.05#Inp.delta_t NuMad_geom_xlscsv_file_twr = "$(path)$(Inp.NuMad_geom_xlscsv_file_twr)" NuMad_mat_xlscsv_file_twr = "$(path)$(Inp.NuMad_mat_xlscsv_file_twr)" NuMad_geom_xlscsv_file_bld = "$(path)$(Inp.NuMad_geom_xlscsv_file_bld)" @@ -110,8 +110,8 @@ mass_breakout_blds,mass_breakout_twr,system,assembly,sections,AD15bldNdIdxRng, A NuMad_geom_xlscsv_file_strut, NuMad_mat_xlscsv_file_strut, Htwr_base=towerHeight, - strut_twr_mountpoint = [0.25,0.75], - strut_bld_mountpoint = [0.25,0.75], + strut_twr_mountpoint = [0.2,0.8], + strut_bld_mountpoint = [0.2,0.8], ntelem, nbelem, ncelem, diff --git a/src/OWENS.jl b/src/OWENS.jl index 11a00e9..bc0c85a 100644 --- a/src/OWENS.jl +++ b/src/OWENS.jl @@ -32,6 +32,7 @@ using WriteVTK const module_path = splitdir(@__FILE__)[1] # Path to this module Modal = OWENSFEA.modal +AutoCampbellDiagram = OWENSFEA.autoCampbellDiagram FEAModel = OWENSFEA.FEAModel Mesh = OWENSFEA.Mesh El = OWENSFEA.El diff --git a/src/SetupTurbine.jl b/src/SetupTurbine.jl index 93a9bf4..4aa89f6 100644 --- a/src/SetupTurbine.jl +++ b/src/SetupTurbine.jl @@ -55,6 +55,8 @@ function setupOWENS(OWENSAero,path; meshtype = "Darrieus", custommesh = nothing) #Darrieus, H-VAWT, ARCUS + custom_mesh_outputs = [] + if AModel=="AD" AD15On = true else @@ -127,7 +129,7 @@ function setupOWENS(OWENSAero,path; verbosity=0, # 0 nothing, 1 basic, 2 lots: amount of printed information connectBldTips2Twr) elseif custommesh != nothing - mymesh, myort, myjoint, AD15bldNdIdxRng, AD15bldElIdxRng = custommesh(;Htwr_base, + mymesh, myort, myjoint, AD15bldNdIdxRng, AD15bldElIdxRng, custom_mesh_outputs = custommesh(;Htwr_base, Htwr_blds, Hbld = H, #blade height R, # m bade radius @@ -637,13 +639,13 @@ function setupOWENS(OWENSAero,path; stiff_twr, stiff_bld,bld_precompinput, bld_precompoutput,plyprops_bld,numadIn_bld,lam_U_bld,lam_L_bld, twr_precompinput,twr_precompoutput,plyprops_twr,numadIn_twr,lam_U_twr,lam_L_twr,aeroForcesAD,deformAeroAD, - mass_breakout_blds,mass_breakout_twr,system,assembly,sections,AD15bldNdIdxRng, AD15bldElIdxRng + mass_breakout_blds,mass_breakout_twr,system,assembly,sections,AD15bldNdIdxRng, AD15bldElIdxRng, custom_mesh_outputs else return mymesh,myel,myort,myjoint,sectionPropsArray,mass_twr, mass_bld, stiff_twr, stiff_bld,bld_precompinput, bld_precompoutput,plyprops_bld,numadIn_bld,lam_U_bld,lam_L_bld, twr_precompinput,twr_precompoutput,plyprops_twr,numadIn_twr,lam_U_twr,lam_L_twr,aeroForcesACDMS,deformAeroACDMS, - mass_breakout_blds,mass_breakout_twr,system,assembly,sections,AD15bldNdIdxRng, AD15bldElIdxRng + mass_breakout_blds,mass_breakout_twr,system,assembly,sections,AD15bldNdIdxRng, AD15bldElIdxRng, custom_mesh_outputs end end diff --git a/src/Unsteady.jl b/src/Unsteady.jl index 8d6b9ac..45f10bc 100644 --- a/src/Unsteady.jl +++ b/src/Unsteady.jl @@ -764,72 +764,6 @@ end function structuralDynamicsTransientGX(topModel,mesh,Fexternal,ForceDof,system,assembly,t,Omega_j,OmegaDot_j,delta_t,numIterations,i,strainGX,curvGX) - function setPrescribedConditions(mesh;pBC=zeros(2,2),Fexternal=[],ForceDof=[]) - Fx = Fexternal[1:6:end] - Fy = Fexternal[2:6:end] - Fz = Fexternal[3:6:end] - Mx = Fexternal[4:6:end] - My = Fexternal[5:6:end] - Mz = Fexternal[6:6:end] - prescribed_conditions = Dict() - for inode = 1:mesh.numNodes - ux = nothing - uy = nothing - uz = nothing - theta_x = nothing - theta_y = nothing - theta_z = nothing - Fx_follower = nothing - Fy_follower = nothing - Fz_follower = nothing - Mx_follower = nothing - My_follower = nothing - Mz_follower = nothing - if inode in pBC[:,1] - for iBC = 1:length(pBC[:,1]) - if pBC[iBC,1] == inode - if pBC[iBC,2] == 1 - ux = pBC[iBC,3] - elseif pBC[iBC,2] == 2 - uy = pBC[iBC,3] - elseif pBC[iBC,2] == 3 - uz = pBC[iBC,3] - elseif pBC[iBC,2] == 4 - theta_x = pBC[iBC,3] - elseif pBC[iBC,2] == 5 - theta_y = pBC[iBC,3] - elseif pBC[iBC,2] == 6 - theta_z = pBC[iBC,3] - end - end - end - end - if isnothing(ux) && !isempty(Fexternal) && Fx[inode]!=0.0 - Fx_follower = Fx[inode] - end - if isnothing(uy) && !isempty(Fexternal) && Fy[inode]!=0.0 - Fy_follower = Fy[inode] - end - if isnothing(uz) && !isempty(Fexternal) && Fz[inode]!=0.0 - Fz_follower = Fz[inode] - end - if isnothing(theta_x) && !isempty(Fexternal) && Mx[inode]!=0.0 - Mx_follower = Mx[inode] - end - if isnothing(theta_y) && !isempty(Fexternal) && My[inode]!=0.0 - My_follower = My[inode] - end - if isnothing(theta_z) && !isempty(Fexternal) && Mz[inode]!=0.0 - Mz_follower = Mz[inode] - end - - prescribed_conditions[inode] = GXBeam.PrescribedConditions(;ux,uy,uz,theta_x,theta_y,theta_z,Fx_follower,Fy_follower,Fz_follower,Mx_follower,My_follower,Mz_follower) - - end - return prescribed_conditions - end - - linear_velocity = [0.0,0.0,0.0] angular_velocity = [0.0,0.0,Omega_j*2*pi] linear_acceleration = [0.0,0.0,0.0] @@ -852,12 +786,12 @@ function structuralDynamicsTransientGX(topModel,mesh,Fexternal,ForceDof,system,a initialize = false if i == 1 && numIterations == 1 # Do initial solve without external loads - prescribed_conditions = setPrescribedConditions(mesh;pBC=topModel.BC.pBC) + prescribed_conditions = GyricFEA.setPrescribedConditions(mesh;pBC=topModel.BC.pBC) system, state, converged = GXBeam.steady_state_analysis!(system,assembly; prescribed_conditions, gravity,angular_velocity,linear=false,reset_state=true) end - prescribed_conditions = setPrescribedConditions(mesh;pBC=topModel.BC.pBC,Fexternal,ForceDof) + prescribed_conditions = GyricFEA.setPrescribedConditions(mesh;pBC=topModel.BC.pBC,Fexternal,ForceDof) initial_state = GXBeam.AssemblyState(system, assembly; prescribed_conditions) systemout, history, converged = GXBeam.time_domain_analysis!(deepcopy(system),assembly, tvec; diff --git a/src/meshing_utilities.jl b/src/meshing_utilities.jl index 2a70f51..c7468a2 100644 --- a/src/meshing_utilities.jl +++ b/src/meshing_utilities.jl @@ -267,7 +267,7 @@ function mesh_beam_centered(;L1 = 6.0, #first section of beam length end """ -return mymesh, myort, myjoint, AD15bldNdIdxRng, AD15bldElIdxRng = create_mesh_struts(;Htwr_base = 15.0, +mymesh, myort, myjoint, AD15bldNdIdxRng, AD15bldElIdxRng = create_mesh_struts(;Htwr_base = 15.0, Htwr_blds = 147.148-15.0, Hbld = 147.148-15.0, #blade height R = 54.014, # m bade radius diff --git a/test/Fig4_5_campbell2.jl b/test/Fig4_5_campbell2.jl index c064078..f9ace4b 100644 --- a/test/Fig4_5_campbell2.jl +++ b/test/Fig4_5_campbell2.jl @@ -26,8 +26,6 @@ path = runpath = splitdir(@__FILE__)[1] nothing -# Unpack inputs, or you could directly input them here and bypass the file - verbosity = 1 turbineType = "Darrieus" @@ -151,8 +149,6 @@ mass_breakout_blds,mass_breakout_twr,system, assembly, sections = OWENS.setupOWE # PyPlot.xlabel("x") # PyPlot.ylabel("y") -starttime = time() - # node, dof, bc top_idx = Int(myjoint[7,2]) pBC = [1 1 0 @@ -172,10 +168,8 @@ top_idx 2 0] ############################################## # Modal Test ############################################# -displ = zeros(mymesh.numNodes*6) -numModes = 16 -mymodel = OWENSFEA.FEAModel(;analysisType = "M", +FEAinputs = OWENSFEA.FEAModel(;analysisType = "M", # outFilename = "$path/data/outplat.out", joint = myjoint, platformTurbineConnectionNodeNumber = 1, @@ -186,138 +180,37 @@ mymodel = OWENSFEA.FEAModel(;analysisType = "M", spinUpOn = true, iterationType = "NR", numNodes = mymesh.numNodes, - numModes) # number of modes to calculate) - - - -# Campbell Diagram generation -rotSpdArrayRPM = 0.0:5.0:40 # rpm -rotorSpeedArrayHZ = rotSpdArrayRPM ./ 60.0 -centStiff = true # centripetal stiffening -NperRevLines = 8 -freq = zeros(length(rotorSpeedArrayHZ),numModes) -for i=1:length(rotorSpeedArrayHZ) - println("$i of $(length(rotorSpeedArrayHZ))") - rotorSpeed = rotorSpeedArrayHZ[i] - local Omega = rotorSpeed - global displ - OmegaStart = Omega - - freqtemp,damp,imagCompSign,U_x_0,U_y_0,U_z_0,theta_x_0,theta_y_0,theta_z_0,U_x_90,U_y_90,U_z_90,theta_x_90,theta_y_90,theta_z_90,displ=OWENS.Modal(mymodel,mymesh,myel;displ,Omega,OmegaStart,returnDynMatrices=false) - - freq[i,:] = freqtemp[1:numModes] -end - -elapsedtime = time() - starttime -println(freq[1,:]) + numModes = 16) # number of modes to calculate) -######################################## -############ GXBeam #################### -######################################## +starttime = time() +freq = OWENS.AutoCampbellDiagram(FEAinputs,mymesh,myel,system,assembly,sections; + minRPM = 0.0, + maxRPM = 40.0, + NRPM = 9, # int + ) +freqOWENS = [freq[:,i] for i=1:4:FEAinputs.numModes-2] +elapsedtime = time() - starttime starttime2 = time() - -# create dictionary of prescribed conditions -prescribed_conditions = Dict( - # fixed base - 1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0), - # fixed top, but free to rotate around z-axis - top_idx => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0), -) - -# --- Perform Analysis --- # -# revolutions per minute -rpm = 0:5:40 - -# gravity vector -gravity = [0, 0, -9.81] - -# number of modes -nmode = 10 - -# number of eigenvalues -nev = 2*nmode - -# storage for results -freq2 = zeros(length(rpm), nmode) -λ = [] -eigenstates =[] -# perform an analysis for each rotation rate -for (i,rpm) in enumerate(rpm) - - global system, Up, λ, eigenstates - - # set turbine rotation - angular_velocity = [0, 0, rpm*(2*pi)/60] - - # eigenvalues and (right) eigenvectors - system, λ, V, converged = eigenvalue_analysis!(system, assembly; - prescribed_conditions = prescribed_conditions, - angular_velocity = angular_velocity, - gravity = gravity, - nev = nev - ) - - # check convergence - @assert converged - - if i > 1 - # construct correlation matrix - C = Up*system.M*V - - # correlate eigenmodes - perm, corruption = correlate_eigenmodes(C) - - # re-arrange eigenvalues - λ = λ[perm] - - # update left eigenvector matrix - Up = left_eigenvectors(system, λ, V) - Up = Up[perm,:] - else - # update left eigenvector matrix - Up = left_eigenvectors(system, λ, V) - end - - # save frequencies - freq2[i,:] = [imag(λ[k])/(2*pi) for k = 1:2:nev] - - state = AssemblyState(system, assembly; - prescribed_conditions = prescribed_conditions) - eigenstates = [AssemblyState(V[:,k],system, assembly; - prescribed_conditions = prescribed_conditions) for k = 1:nev] - -end - +FEAinputs.analysisType = "GX" +freq2 = OWENS.AutoCampbellDiagram(FEAinputs,mymesh,myel,system,assembly,sections; + minRPM = 0.0, + maxRPM = 40.0, + NRPM = 9, # int + vtksavename="$path/campbellVTK/SNL34m", + saveModes = [1,3,5], #must be int + mode_scaling = 500.0, + ) +freqGX = [freq2[:,i] for i=1:2:FEAinputs.numModes-6-2] elapsedtime2 = time() - starttime2 println("OWENS: $elapsedtime") println("GX: $elapsedtime2") -state = AssemblyState(system, assembly; prescribed_conditions=prescribed_conditions) - -filename = "$path/vtk/Campbell5_34m" -try #this should error if someone on windows uses backslash '\' - lastforwardslash = findlast(x->x=='/',filename) - filepath = filename[1:lastforwardslash-1] - if !isdir(filepath) - mkdir(filepath) - end -catch - @info "Please manually create the directory to house $filename" -end - -write_vtk(filename, assembly, state, - λ[5], eigenstates[5]; sections,mode_scaling = 500.0) - - -freqOWENS = [freq[:,i] for i=1:4:numModes-2] -freqGX = [freq2[:,i] for i=1:2:numModes-6-2] - for imode = 1:length(freqOWENS) for ifreq = 1:length(freqOWENS[imode]) - println("imode $imode, ifreq $ifreq") + # println("imode $imode, ifreq $ifreq") atol = freqGX[imode][ifreq]*0.065 @test isapprox(freqOWENS[imode][ifreq],freqGX[imode][ifreq];atol) end @@ -337,9 +230,11 @@ end # # PyPlot.rc("axes", prop_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"]) # plot_cycle=["#348ABD", "#A60628", "#009E73", "#7A68A6", "#D55E00", "#CC79A7"] +# NperRevLines = 8 +# rotSpdArrayRPM = LinRange(0.0, 40.0, 9) # int # PyPlot.figure() -# for i=1:1:numModes-2 -# PyPlot.plot(rotSpdArrayRPM,freq[:,i],color=plot_cycle[1],"b-") #plot mode i at various rotor speeds +# for i=1:1:FEAinputs.numModes-2 +# PyPlot.plot(rotSpdArrayRPM,freq[:,i],color=plot_cycle[1],"-") #plot mode i at various rotor speeds # end # # Plot 34m experimental data @@ -381,8 +276,8 @@ end # # PyPlot.savefig("$(path)/../figs/34mCampbell.pdf",transparent = true) # # Add to figure -# for i=1:2:numModes-6-2 -# PyPlot.plot(rpm,freq2[:,i],color=plot_cycle[2],"-") #plot mode i at various rotor speeds +# for i=1:2:FEAinputs.numModes-6-2 +# PyPlot.plot(rotSpdArrayRPM,freq2[:,i],color=plot_cycle[2],"-") #plot mode i at various rotor speeds # end # PyPlot.plot(0,0,color=plot_cycle[2],"-",label="GXBeam") # PyPlot.legend(fontsize=8.5,loc = (0.09,0.8),ncol=2,handleheight=1.8, labelspacing=0.03)