From 9fa8614c3f1da8998e9d446493b0ad770bf45f70 Mon Sep 17 00:00:00 2001 From: kevmoor Date: Wed, 11 Sep 2024 16:58:27 -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 | 2 + src/OWENSFEA.jl | 1 + src/modal.jl | 163 +++++++++++++++++++++++++++++++++++++++++++++++ src/utilities.jl | 79 +++++++++++++++++++++++ 4 files changed, 245 insertions(+) diff --git a/.gitignore b/.gitignore index e61a538..bf542c7 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,5 @@ *.vtu *.pvd *.pdf +*.vscode +*.mat \ No newline at end of file diff --git a/src/OWENSFEA.jl b/src/OWENSFEA.jl index 7fb750e..4df521d 100644 --- a/src/OWENSFEA.jl +++ b/src/OWENSFEA.jl @@ -7,6 +7,7 @@ import DelimitedFiles import Printf import MAT import StaticArrays +import GXBeam export mapACloads, calculateStructureMassProps, createJointTransform, calculateReducedDOFVector, constructReducedDispVectorMap, calculateBCMap, calculateElementMass!, ConcMassAssociatedWithElement, applyBC, calculateReactionForceAtNode, calculateStrainForElements, findElementsAssociatedWithNodeNumber, applyGeneralConcentratedTerms diff --git a/src/modal.jl b/src/modal.jl index 501312f..93147c5 100644 --- a/src/modal.jl +++ b/src/modal.jl @@ -1,3 +1,166 @@ +""" +frequencies = autoCampbellDiagram(FEAinputs,mymesh,myel,system,assembly; + minRPM = 0.0, + maxRPM = 40.0, + NRPM = 9, # int + vtksavename = nothing, + saveModes = [1,3,5], #must be int + saveRPM = [1,3,5], #must be int + mode_scaling = 500.0, + ) + +Automated Campbell Diagram Generator, this function runs the model with centrifugal stiffening for the specified RPM levels. +If FEAinputs.analysisType == "GX" and vtksavename are specified, it will output paraview mode shape files at the specified save name. + +#Inputs +* `FEAinputs::OWENSFEA.FEAModel`: The FEA modeling options +* `mymesh::OWENSFEA.Mesh`: a previously generated turbine mesh +* `myel::OWENSFEA.El`: the element properties associated with that mesh +* `system::GXBeam.System`: the converted GXBeam system from the mesh and el +* `assembly::GXBeam.AssemblyState`: the converted GXBeam assembly from the mesh and el +* `sections::Array{Float64, 3}`: the 3D point cloud to be converted to VTK format +* `minRPM::Float64`: minimum RPM to be run, e.x. 0.0 +* `maxRPM::Float64`: maximum RPM to be run e.x. 40.0 +* `NRPM::Int64`: number of linear discretizations of RPM e.x. 9 must be int +* `vtksavename::string`: filename (with path if desired) of the VTK outputs if GX. Set to "nothing" to not save. +* `saveModes::Array{Int64}`: The modes to save in the VTK outputs e.x. [1,3,5] must be int +* `saveRPM::Array{Int64}`: The RPMs to save in the VTK outputs e.x. [1,3,5] must be int +* `mode_scaling::Float64`: The mode scaling in the VTK outputs e.x. 500.0 + +#Outputs +* `frequency::Array{Float64}`: The output modal frequencies +""" +function autoCampbellDiagram(FEAinputs,mymesh,myel,system,assembly,sections; + minRPM = 0.0, + maxRPM = 40.0, + NRPM = 9, # int + vtksavename = nothing, + saveModes = [1,3,5], #must be int + saveRPM = [1,3,5], #must be int + mode_scaling = 500.0, + ) + + rotSpdArrayRPM = LinRange(minRPM,maxRPM,NRPM) + + # Campbell Diagram generation + if FEAinputs.analysisType!="GX" + rotorSpeedArrayHZ = rotSpdArrayRPM ./ 60.0 + freq = zeros(length(rotorSpeedArrayHZ),FEAinputs.numModes) + for irpm=1:length(rotorSpeedArrayHZ) + println("$irpm of $(length(rotorSpeedArrayHZ))") + rotorSpeed = rotorSpeedArrayHZ[irpm] + local Omega = rotorSpeed + displ = zeros(mymesh.numNodes*6) + 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=modal(FEAinputs, + mymesh,myel;displ,Omega,OmegaStart,returnDynMatrices=false) + + freq[irpm,:] = freqtemp[1:FEAinputs.numModes] + end + + return freq + else + ######################################## + ############ GXBeam #################### + ######################################## + + prescribed_conditions = setPrescribedConditions(mymesh;pBC=FEAinputs.BC.pBC)#,Fexternal,ForceDof) + + # 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 --- # + + # gravity vector + gravity = [0, 0, -9.81] #TODO: from FEAinputs + + # number of modes + nmode = FEAinputs.numModes + + # number of eigenvalues + nev = 2*nmode + + # storage for results + freq2 = zeros(length(rotSpdArrayRPM), nmode) + λ_save = Array{Array{ComplexF64,1},1}(undef,length(rotSpdArrayRPM)) + eigenstates_save = Array{Array{GXBeam.AssemblyState,1},1}(undef,length(rotSpdArrayRPM)) + Up = [] + # perform an analysis for each rotation rate + for (irpm,rpm) in enumerate(rotSpdArrayRPM) + println("$irpm of $(length(rotorSpeedArrayHZ))") + # set turbine rotation + angular_velocity = [0, 0, rpm*(2*pi)/60] + + # eigenvalues and (right) eigenvectors + system, λ, V, converged = GXBeam.eigenvalue_analysis!(system, assembly; + prescribed_conditions = prescribed_conditions, + angular_velocity = angular_velocity, + gravity = gravity, + nev = nev + ) + + # check convergence + @assert converged + + if irpm > 1 + # construct correlation matrix + C = Up*system.M*V + + # correlate eigenmodes + perm, corruption = GXBeam.correlate_eigenmodes(C) + + # re-arrange eigenvalues + λ = λ[perm] + + # update left eigenvector matrix + Up = GXBeam.left_eigenvectors(system, λ, V) + Up = Up[perm,:] + else + # update left eigenvector matrix + Up = GXBeam.left_eigenvectors(system, λ, V) + end + + # save frequencies + freq2[irpm,:] = [imag(λ[k])/(2*pi) for k = 1:2:nev] + + state = GXBeam.AssemblyState(system, assembly; + prescribed_conditions = prescribed_conditions) + eigenstates = [GXBeam.AssemblyState(V[:,k],system, assembly; + prescribed_conditions = prescribed_conditions) for k = 1:nev] + λ_save[irpm] = λ + eigenstates_save[irpm] = eigenstates + end + + if !isnothing(vtksavename) #TODO: map the OWENS state into the gx state for filesaving + state = GXBeam.AssemblyState(system, assembly; prescribed_conditions=prescribed_conditions) + + try #this should error if someone on windows uses backslash '\' + lastforwardslash = findlast(x->x=='/',vtksavename) + filepath = vtksavename[1:lastforwardslash-1] + if !isdir(filepath) + mkdir(filepath) + end + catch + @info "Please manually create the directory to house $vtksavename" + end + for isaveRPM in saveRPM + for isavemode in saveModes + GXBeam.write_vtk("$(vtksavename)_RPM$(rotSpdArrayRPM[isaveRPM])_Mode$(isavemode)", assembly, state, + λ_save[isaveRPM][isavemode], eigenstates_save[isaveRPM][isavemode]; sections,mode_scaling) + end + end + end + + return freq2 + end +end + """ modal(feamodel,mesh,el;Omega=0.0,displ=zeros(mesh.numNodes*6),OmegaStart=0.0,returnDynMatrices=false) diff --git a/src/utilities.jl b/src/utilities.jl index aec0ca1..bc61acb 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1856,3 +1856,82 @@ function applyConcentratedTerms(numNodes, numDOFPerNode; filename="none",data=[] return NodalTerms(concLoad,concStiff,concMass,concDamp) end + +""" + prescribed_conditions = setPrescribedConditions(mesh;pBC=zeros(2,2),Fexternal=[],ForceDof=[]) + +Internal, maps OWENS boundary conditions and applied forces to the GXBeam PrescribedConditions input + +#Input +* `mesh::FEAModel.mesh`: Input turbine mesh +* `pBC::FEAModel.BC.pBC`: Boundary conditions +* `Fexternal::Array{Float64}`: Applied forces to the mesh, 1D array ordered node 1 Dof 1-6, node 2 Dof 1-2, etc. +* `ForceDof::Array{Float64}`: Degrees of freedom aligned with Fexternal, currently is unused, assumes Fexternal uses the full DOF array. + +#Output +* `prescribed_conditions::GXBeam.PrescribedConditions`: the boundary conditions/applied forces used by GXBeam see ?GXBeam.PrescribedConditions +""" +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 \ No newline at end of file