[ANN] GXBeam.jl v0.3.0 - Geometrically Exact Beam Theory

taylormcd · November 11, 2021, 4:18pm

I’m pleased to announce the release of version 0.3 of the GXBeam package, which contains a thoroughly verified and validated pure Julia implementation of geometrically exact beam theory.

This version introduces:

Point masses/rigid bodies
Zero length beam elements
Massless beam elements
Native support for fictitious loads caused by body frame linear/angular acceleration
Native support for gravitational loads

I also added a new validation case to the examples which compares the computational results of GXBeam to experimental data for the Sandia 34-Meter Vertical Axis Wind Turbine, pictured here:
vawt

Here’s the results:

comparison

This example was originally presented to me by one of my former colleagues who now works at Sandia. They asked that the following citation accompany it:

Moore, K. and Ennis, B., “Aeroelastic Validation of the Offshore Wind Energy Simulator for Vertical-Axis Wind Turbines”, forthcoming 2022

Ribeiro · November 11, 2021, 9:45pm

Hi, Taylor. This is exciting stuff!
I’ve been doing some workarounds to make up for not having lumped masses and immediately tested the new features, had some trouble.
Here’s a tiny example:

using GXBeam, StaticArrays, LinearAlgebra
stiff = [1000000 0 0 -.5 -1 -50000; 0 3000000 0 0 0 0; 0 0 3000000 0 0 0; 0 0 0 7 .1 -.02; 0 0 0 0 5 .1; 0 0 0 0 0 3000]
stiff = Symmetric(stiff)
mass = [.02 0 0 0 -5e-7 -1e-7; 0 .02 0 5e-7 0 1e-4; 0 0 .02 1e-7 -1e-4 0; 0 0 0 1e-5 1e-8 2e-10; 0 0 0 0 6e-7 9e-9; 0 0 0 0 0 1e-5]
mass = Symmetric(mass)

nodes = [SA_F64[0,0,0], SA_F64[0,.25,0],SA_F64[0,.5,0]]

start = [1,2,2,3]
stop =  [2,2,3,3]

stiffness = [stiff,zeros(6,6),stiff,zeros(6,6)]
pointmass=Dict(2=>PointMass(mass),4=>PointMass(mass))

transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:4];

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);

ERROR: LoadError: SingularException(0)
Stacktrace:
 [1] lu(S::SparseArrays.SparseMatrixCSC{Float64, Int64}; check::Bool)
   @ SuiteSparse.UMFPACK C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\SuiteSparse\src\umfpack.jl:162
 [2] #safe_lu#309
   @ C:\Users\apo5\.julia\packages\GXBeam\iY3sc\src\interfaces\forwarddiff.jl:62 [inlined]
 [3] safe_lu
   @ C:\Users\apo5\.julia\packages\GXBeam\iY3sc\src\interfaces\forwarddiff.jl:62 [inlined]
 [4] eigenvalue_analysis!(system::System{Float64, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}}, assembly::Assembly{Float64, Vector{SVector{3, Float64}}, Vector{Int64}, Vector{GXBeam.Element{Float64}}}; prescribed_conditions::Dict{Int64, PrescribedConditions{Float64}}, distributed_loads::Dict{Int64, DistributedLoads{Float64}}, point_masses::Dict{Int64, PointMass{Float64}}, gravity::SVector{3, Int64}, linear::Bool, linearization_state::Nothing, update_linearization_state::Bool, method::Symbol, linesearch::LineSearches.BackTracking{Float64, Int64}, ftol::Float64, iterations::Int64, reset_state::Bool, find_steady_state::Bool, origin::SVector{3, Float64}, linear_velocity::SVector{3, Float64}, angular_velocity::SVector{3, Float64}, linear_acceleration::SVector{3, Float64}, angular_acceleration::SVector{3, Float64}, tvec::Float64, nev::Int64)
   @ GXBeam C:\Users\apo5\.julia\packages\GXBeam\iY3sc\src\analyses.jl:593
 [5] #eigenvalue_analysis#237
   @ C:\Users\apo5\.julia\packages\GXBeam\iY3sc\src\analyses.jl:460 [inlined]
 [6] top-level scope

Maybe I misunderstand how to use the lumped masses.
BTW, I’ll likely publish some FSI stuff using GXBeam in the near future. Any preferences on how to reference it?
Thanks!

taylormcd · November 12, 2021, 4:54pm

I looked into this and the current solution is to use two separate nodes for the zero length elements rather than using the same node. I’ll probably work on getting your version of this example to work in a future minor release. Here’s how your example would change.

using GXBeam, StaticArrays, LinearAlgebra
stiff = [1000000 0 0 -.5 -1 -50000; 0 3000000 0 0 0 0; 0 0 3000000 0 0 0; 0 0 0 7 .1 -.02; 0 0 0 0 5 .1; 0 0 0 0 0 3000]
stiff = Symmetric(stiff)
mass = [.02 0 0 0 -5e-7 -1e-7; 0 .02 0 5e-7 0 1e-4; 0 0 .02 1e-7 -1e-4 0; 0 0 0 1e-5 1e-8 2e-10; 0 0 0 0 6e-7 9e-9; 0 0 0 0 0 1e-5]
mass = Symmetric(mass)

nodes = [SA_F64[0,0,0], SA_F64[0,.25,0], SA_F64[0,.25,0], SA_F64[0,.5,0], SA_F64[0,.5,0]]

start = 1:4
stop =  2:5

stiffness = [stiff,zeros(6,6),stiff,zeros(6,6)]
pointmass=Dict(2=>PointMass(mass),4=>PointMass(mass))

transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:4];

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);

I’d be curious to know how your lumped mass workarounds compare to what I’ve implemented. I don’t have a proper test for the lumped masses, but since they effectively just modify the beam element mass matrix there isn’t much to mess up with their implementation.

With regard to a citation, I don’t currently have a reference that goes with the code so I’ll get back to you on that.

PetrKryslUCSD · November 12, 2021, 7:23pm

Sorry: what is the point of zero-length beams? Implementation of various joints?

taylormcd · November 12, 2021, 7:25pm

They can be convenient for constructing joints. They also allow point masses to be effectively placed between or at the end of elements, like in the above test case.

Ribeiro · November 14, 2021, 8:10pm

Hi, Taylor.
Your method works, in the sense that it gives me an answer, but it’s really different from what I expect.
Here’s a demo:

using GXBeam, StaticArrays, LinearAlgebra
nodes = [SA_F64[0,i,0] for i in 0:.1:1]

nElements = length(nodes)-1
start = 1:nElements
stop =  2:(nElements+1)
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:nElements];

stiffness = [stiff for i in 1:nElements]
massmatrix = [mass./0.1 for i in 1:nElements] # here I divide by the element length, which I don't do for the lumped mass

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation, mass=massmatrix);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10);
frequencies = imag(λ[1:2:10])/2π


nodes = sort(vcat(nodes,nodes[2:end]))

nElements = length(nodes)-1
start = 1:nElements
stop =  2:(nElements+1)
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:nElements];

stiffness = [i%2==0 ? 0*stiff : stiff for i in 1:nElements]
pointmass = Dict(2 => PointMass(mass))
for i in 4:2:nElements
    pointmass[i] = PointMass(mass)
end
pointmass[nElements] = PointMass(mass./2) # last lumped mass is half of the others, as it represents the last half of an element

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);
frequencies = imag(λ[1:2:10])/2π

First and then second answers:

frequencies
5-element Vector{Float64}:
  2.778377726671569
 18.025106038785804
 27.589923115564932
 53.92849024854462
 66.28027798190845

5-element Vector{Float64}:
 1.144955956537015        
 1.4777550582703745       
 3.4326136503753086       
 4.396883716032072        
 5.981719555816847

I didn’t check, but the first seems to make more sense than the second and I’d expect them to be close.
Am I doing something wrong? I’d guess so =].
Thanks a lot!

taylormcd · November 15, 2021, 6:01pm

Looks like there’s a bug in my math for the point masses. I’ll get a new release with the fix released today. I’ll post an update here when I’m done. Thanks for finding this.

PetrKryslUCSD · November 15, 2021, 7:14pm

Out of curiosity: would your code handle a fast spinning Lagrangian top modeled as a beam?

taylormcd · November 15, 2021, 10:03pm

I’m not sure. I suppose it depends on how it is modeled. The only limitation I can think of that impacts this use case is the fact that only rigid joints are supported currently.

taylormcd · November 15, 2021, 10:31pm

I released a new version with the point mass bug fix. I also added a unit test. Keep in mind that the mass matrix for point masses are defined in the body frame rather than the element frame, which changes your example a little bit since you need to rotate the mass matrices for the lumped masses. Here’s the updated version:

using GXBeam, StaticArrays, LinearAlgebra
stiff = [1000000 0 0 -.5 -1 -50000; 0 3000000 0 0 0 0; 0 0 3000000 0 0 0; 0 0 0 7 .1 -.02; 0 0 0 0 5 .1; 0 0 0 0 0 3000]
stiff = Symmetric(stiff)
mass = [.02 0 0 0 -5e-7 -1e-7; 0 .02 0 5e-7 0 1e-4; 0 0 .02 1e-7 -1e-4 0; 0 0 0 1e-5 1e-8 2e-10; 0 0 0 0 6e-7 9e-9; 0 0 0 0 0 1e-5]
mass = Symmetric(mass)

nodes = [SA_F64[0,i,0] for i in 0:.1:1]

nElements = length(nodes)-1
start = 1:nElements
stop =  2:(nElements+1)
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:nElements];

stiffness = [stiff for i in 1:nElements]
massmatrix = [mass./0.1 for i in 1:nElements] # here I divide by the element length, which I don't do for the lumped mass

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation, mass=massmatrix);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10);
frequencies = imag(λ[1:2:10])/2π


nodes = sort(vcat(nodes,nodes[2:end]))

nElements = length(nodes)-1
start = 1:nElements
stop =  2:(nElements+1)
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:nElements];

stiffness = [i%2==0 ? 0*stiff : stiff for i in 1:nElements]

pointmass = Dict(2 => PointMass(GXBeam.transform_properties(mass, transformation[2]')))
for i in 4:2:nElements
    pointmass[i] = PointMass(GXBeam.transform_properties(mass, transformation[i]'))
end
pointmass[nElements] = PointMass(GXBeam.transform_properties(mass, transformation[nElements]')./2) # last lumped mass is half of the others, as it represents the last half of an element

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);
frequencies = imag(λ[1:2:10])/2π

First and then second answers:

5-element Vector{Float64}:
  2.778377726671564
 18.02510603873221
 27.589923228529937
 53.92848744691813
 66.28028876304779

5-element Vector{Float64}:
  2.762457025865554
 17.31061493577968
 27.432545128177775
 48.88599842911449
 66.07594276651227

Here’s how to get the exact same results using the point mass matrices instead of element mass matrices, which is informative for showing how they’re implemented.

using GXBeam, StaticArrays, LinearAlgebra
stiff = [1000000 0 0 -.5 -1 -50000; 0 3000000 0 0 0 0; 0 0 3000000 0 0 0; 0 0 0 7 .1 -.02; 0 0 0 0 5 .1; 0 0 0 0 0 3000]
stiff = Symmetric(stiff)
mass = [.02 0 0 0 -5e-7 -1e-7; 0 .02 0 5e-7 0 1e-4; 0 0 .02 1e-7 -1e-4 0; 0 0 0 1e-5 1e-8 2e-10; 0 0 0 0 6e-7 9e-9; 0 0 0 0 0 1e-5]
mass = Symmetric(mass)

nodes = [SA_F64[0,i,0] for i in 0:.1:1]

nElements = length(nodes)-1
start = 1:nElements
stop =  2:(nElements+1)
transformation = [[0 -1 0; 1 0 0; 0 0 1] for _ in 1:nElements];

stiffness = [stiff for i in 1:nElements]
massmatrix = [mass./0.1 for i in 1:nElements] # here I divide by the element length, which I don't do for the lumped mass

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation, mass=massmatrix);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10);
frequencies = imag(λ[1:2:10])/2π

pointmass = Dict(1 => PointMass(GXBeam.transform_properties(mass, transformation[1]')))
for i in 2:nElements
    pointmass[i] = PointMass(GXBeam.transform_properties(mass, transformation[i]'))
end

assembly = GXBeam.Assembly(nodes, start, stop, stiffness=stiffness, frames=transformation);
prescribed_conditions = Dict(1 => GXBeam.PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0));
system, λ, V, converged = GXBeam.eigenvalue_analysis(assembly, prescribed_conditions = prescribed_conditions, nev = 10, point_masses=pointmass);
frequencies = imag(λ[1:2:10])/2π

Ribeiro · November 16, 2021, 9:38am

Yay, it works now!
And this did improve my results! Previously I had up to 5% differences (depending on the mode) between my model and a reference one using lumped masses and now I’m well under 1%.
Thanks for the update and for fixing the bug!

taylormcd · November 16, 2021, 2:56pm

That’s great to hear! I hope GXBeam continues to work well for you.