Merge pull request #21 from sandialabs/f/AutoCampbell

Add and break out functions that automatically map OWENS boundary con…

.gitignore

@@ -4,3 +4,5 @@
 *.vtu
 *.pvd
 *.pdf
+*.vscode
+*.mat

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
 

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)
 

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