byuflowlab / gxbeam.jl Goto Github PK

View Code? Open in Web Editor NEW

87.0 5.0 18.0 72.48 MB

Pure Julia Implementation of Geometrically Exact Beam Theory

License: MIT License

Julia 98.64% TeX 1.36%

gebt beam-elements finite-element-analysis julia nonlinear

gxbeam.jl's Introduction

GXBeam

Pure Julia Implementation of Geometrically Exact Beam Theory

Author: Taylor McDonnell

GXBeam is a pure Julia implementation of Geometrically Exact Beam Theory, originally based on the open source code GEBT and its associated papers[1][2], though it has since been augmented with a number of additional features.

As a sample of one of the many things this package can do, here's a time domain simulation of the dynamic response of a joined wing subjected to a simulated gust, scaled up in order to visualize the deflections:

And here's a dynamic simulation of a wind turbine subjected to a sinusoidal tip load.

Package Features

Performs multiple types of analyses including:
- Linear/Nonlinear static analyses
- Linear/Nonlinear steady-state analyses
- Linear/Nonlinear eigenvalue analyses (by linearizing about a steady state condition)
- Linear/Nonlinear time-marching dynamic analyses
Accurately models arbitrary systems of interconnected highly flexible composite beams.
- Captures all geometric nonlinearities due to large deflections and rotations (subject to a small strain assumption)
- Models angular displacements of any magnitude using only three parameters
- Uses the full 6x6 Timoshenko beam stiffness matrix
Calculate section compliance and inertia matrices and compute strain recovery
- Uses quadrilateral finite elements rather than classical lamiante theory for much better accuracy and cross coupling
- Allows for general geometry with inhomogenous properties and anisotropic behavior (computes full 6x6 matrix)
- Ply materials are general orthotropic
- Provides convenience method for paramterizing airfoil layups
- Strain recovery functionality to compute stresses and strains on the mesh using the results of the beam solution
- Tsai-wu failure criteria
Models time-varying distributed forces/moments including
- Point and distributed loads which remain fixed in the body-frame
- Point and distributed loads which rotate with the structure
- Loads due to known body frame velocities and accelerations
- Gravitational loads acting on beam elements and point masses
- Loads resulting from stiffness-proportional structural damping
Optional DifferentialEquations interface.
- Constant mass matrix differential algebraic equation formulation
- Fully implicit differential algebraic equation formulation
Provides derivatives with ForwardDiff (including overloading internal solvers with implicit analytic methods)
Result visualization using WriteVTK
Verified and validated against published analytical and computational results. See the examples in the documentation.

Installation

Enter the package manager by typing ] and then run the following:

pkg> add GXBeam

Performance

This code has been optimized to be highly performant. In our tests we found that GXBeam outperforms GEBT by a significant margin across all analysis types, as seen in the following table. More details about the specific cases which we test may be found by inspecting the input files and scripts for these tests in the benchmark folder.

Package	Steady Analysis	Eigenvalue Analysis	Time Marching Analysis
GEBT	13.722 ms	33.712 ms	26.870 s
GXBeam	4.716 ms	18.478 ms	9.019 s

Usage

See the documentation

Limitations

By using the simplest possible shape functions (constant or linear shape functions), this package avoids using numerical quadrature except when integrating applied distributed loads (which can be pre-integrated). As a result, element properties are approximated as constant throughout each beam element and a relatively large number of beam elements may be necessary to achieve grid-independent results. More details about the convergence of this package may be found in the examples.

This package does not currently model cross section warping, and therefore should not be used to model open cross sections (such as I, C, or L-beams). The one exception to this rule is if the beam's width is much greater than its height, in which case the beam may be considered to be strip-like (like a helicopter blade).

This package relies on the results of linear cross-sectional analyses. Most notably, it does not model the nonlinear component of the Trapeze effect, which is the tendency of a beam to untwist when subjected to axial tension. This nonlinear effect is typically most important when modeling rotating structures such as helicopter blades due to the presence of large centrifugal forces. It is also more important when modeling strip-like beams than for modeling closed cross-section beams due to their low torsional rigidity.

Related Codes

GEBT: Open source geometrically exact beam theory code developed in Fortran as a companion to the proprietary cross sectional analysis tool VABS. The theory for this code is provided in references 1 and 2. GXBeam was originally developed based on this package and its associated papers, but has since been augmented with additional features.

BeamDyn: Open source geometrically exact beam theory code developed in Fortran by NREL as part of the OpenFAST project. This code was also developed based on GEBT, but uses Legendre spectral finite elements. This allows for exponential rather than algebraic convergence when the solution is smooth. This makes this code a good candidate for use when analyzing beams with smoothly varying properties. Unfortunately, the code is limited to analyzing a single beam, rather than an assembly of beams.

The cross sectional analysis uses the same underlying theory as in BECAS, but was written to be fast and optimization-friendly. VABS and PreComp are other popular tools for composite cross sectional analysis. The former is not freely available, whereas the latter is lower fidelity as it is based on classical laminate theory.

Contributing

Contributions are welcome and encouraged. If at any point you experience issues or have suggestions related to this package, create a new Github issue so we can discuss it. If you're willing to help solve an issue yourself, we encourage you to create a fork of this repository and submit a pull request with the requested change. Pull requests should generally also add a unit test in test/runtests.jl to ensure that issues do not reoccur along with future changes.

References

[1] Yu, W., & Blair, M. (2012). GEBT: A general-purpose nonlinear analysis tool for composite beams. Composite Structures, 94(9), 2677-2689.

[2] Wang, Q., & Yu, W. (2017). Geometrically nonlinear analysis of composite beams using Wiener-Milenković parameters. Journal of Renewable and Sustainable Energy, 9(3), 033306.

[3] Hodges, D. (2006). Nonlinear Composite Beam Theory. American Institute of Aeronautics and Astronautics.

gxbeam.jl's People

Contributors

Stargazers

Watchers

Forkers

zhuwenguo kongweizhen yang7857854 budhyant yosinlpet ben-isaac christianhvejsel vagrawal91 simo-11 hwfluid ss1870 moe93 downl3824 faisal-bhuiyan axla-io xue-cheng createkaero alpaan

gxbeam.jl's Issues

Zero pivots in factorization error when computing compliance matrix

Here's a basic example that leads to a factorization error with zero pivots.

using GXBeam

xaf = [1.0, 0.99669031, 0.98707731, 0.97177624, 0.95120013, 0.92554805, 0.89513234, 0.86042844, 0.82194932, 0.78024605, 0.73589923, 0.68948289, 0.64155561, 0.59265387, 0.5432852, 0.49390561, 0.44491485, 0.39667271, 0.34949257, 0.30376269, 0.25995386, 0.21850111, 0.17980252, 0.14421934, 0.11207729, 0.08364666, 0.05914156, 0.03873701, 0.02260069, 0.01080931, 0.00329342, 0.0, 0.00480487, 0.0111091, 0.02423919, 0.04202841, 0.0645839, 0.09186189, 0.1236919, 0.15984299, 0.20001108, 0.24382591, 0.29085133, 0.34057636, 0.39242047, 0.4457723, 0.49999633, 0.55444002, 0.60844101, 0.66133467, 0.71246137, 0.7611783, 0.80686551, 0.84892722, 0.8867981, 0.91995456, 0.9479178, 0.97026156, 0.98661057, 0.99662097, 1.0]
yaf = [0.0, 0.00049003, 0.00207033, 0.00475426, 0.00825965, 0.01241627, 0.01723294, 0.0226708, 0.02862516, 0.034964, 0.0415088, 0.0480415, 0.05432933, 0.06013032, 0.06520439, 0.06931828, 0.07229087, 0.07399667, 0.0744182, 0.07365704, 0.07178323, 0.06883168, 0.064851, 0.05990381, 0.05406639, 0.04742811, 0.04011465, 0.03230144, 0.02418269, 0.01593784, 0.00791465, 0.0, -0.00634666, -0.0102234, -0.01556325, -0.02006888, -0.02353976, -0.02598343, -0.02742948, -0.02793392, -0.02758763, -0.026503, -0.02481832, -0.02268597, -0.02024403, -0.01760693, -0.01488444, -0.01217351, -0.00956485, -0.00713548, -0.0049519, -0.0030634, -0.00151205, -0.00032506, 0.00048717, 0.00093779, 0.00105862, 0.00091038, 0.00057426, 0.00018251, 0.0]

chord = 0.0863
twist = 17.7*pi/180
pitch_axis = 0.25 / chord

uni = Material(1.75e11, 8.0e9, 8.0e9, 5.0e9, 5.0e9, 5.0e9, 0.3, 0.3, 0.3, 1600.0)
weave = Material(8.5e10, 8.5e10, 8.5e10, 5.0e9, 5.0e9, 5.0e9, 0.1, 0.1, 0.1, 1600.0)
materials = [uni, weave]

thickness = [0.001, 0.01853, 0.001]
orientation = [45.0,0.0,45.0]*pi/180
material_ids = [2;1;2]
layup1 = Layer.(materials[material_ids], thickness, orientation)

segments = [layup1]

xbreak = [0.0, 1.0]
web_location = [] #no webs
webs = []

nodes, elements = afmesh(xaf, yaf, chord, twist, pitch_axis, xbreak, web_location, segments, webs)
cache = initialize_cache(nodes, elements, eltype(chord))

compliance, _, _ = compliance_matrix(nodes, elements; cache)
mass, _ = mass_matrix(nodes, elements)

Here is the stacktrace:

ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] ldlt!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool)
   @ SuiteSparse.CHOLMOD /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1472
 [2] #ldlt!#11
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
 [3] ldlt!
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
 [4] factorize(A::LinearAlgebra.Symmetric{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}})
   @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1634
 [5] linearsolve2(AM::SparseArrays.SparseMatrixCSC{Float64, Int64}, B2::SparseArrays.SparseMatrixCSC{Float64, Int64})
   @ GXBeam ~/.julia/dev/GXBeam/src/section.jl:559
 [6] compliance_matrix(nodes::Vector{Node{Float64}}, elements::Vector{MeshElement{Vector{Int64}, Float64}}; cache::GXBeam.SectionCache{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Int64}}, gxbeam_order::Bool)
   @ GXBeam ~/.julia/dev/GXBeam/src/section.jl:736
 [7] gxbeam_test2(i::Int64)
   @ Main ~/Dropbox/Research-Ning/projects/Prop_out/src/debugging/gxbeam_test.jl:106
 [8] top-level scope
   @ REPL[30]:1

I'm using GXBeam 0.4.1 and Julia 1.6.6.

I can try to figure out what's going on but from a brief look into it it seems others have gotten this error before (here and here) on older versions of Julia, where it was determined it was a bug in Julia itself that was fixed in subsequent releases. After seeing this I tried it again with Julia 1.8.0 with the same result.

@andrewning did you see this at all during the code development?

Load ramping

Hi, just wondered if you've considered a load ramping method for the nonlinear solve? I'm wondering if it would make the solve more robust, and could it potentially result in convergence with fewer iterations?

NaN compliance matrix for more cambered airfoil

If I run the following example for composite material properties it results in a NaN-filled compliance matrix. I've taken this straight from the docs, only changing the airfoil coordinates. (Proposed solution below)

using GXBeam

# airfoil coordinates in the docs --> work great
# xaf = [1.00000000, 0.99619582, 0.98515158, 0.96764209, 0.94421447, 0.91510964, 0.88074158, 0.84177999, 0.79894110, 0.75297076, 0.70461763, 0.65461515, 0.60366461, 0.55242353, 0.50149950, 0.45144530, 0.40276150, 0.35589801, 0.31131449, 0.26917194, 0.22927064, 0.19167283, 0.15672257, 0.12469599, 0.09585870, 0.07046974, 0.04874337, 0.03081405, 0.01681379, 0.00687971, 0.00143518, 0.00053606, 0.00006572, 0.00001249, 0.00023032, 0.00079945, 0.00170287, 0.00354717, 0.00592084, 0.01810144, 0.03471169, 0.05589286, 0.08132751, 0.11073805, 0.14391397, 0.18067874, 0.22089879, 0.26433734, 0.31062190, 0.35933893, 0.40999990, 0.46204424, 0.51483073, 0.56767889, 0.61998250, 0.67114514, 0.72054815, 0.76758733, 0.81168064, 0.85227225, 0.88883823, 0.92088961, 0.94797259, 0.96977487, 0.98607009, 0.99640466, 1.00000000]
# yaf = [0.00000000, 0.00017047, 0.00100213, 0.00285474, 0.00556001, 0.00906779, 0.01357364, 0.01916802, 0.02580144, 0.03334313, 0.04158593, 0.05026338, 0.05906756, 0.06766426, 0.07571157, 0.08287416, 0.08882939, 0.09329359, 0.09592864, 0.09626763, 0.09424396, 0.09023579, 0.08451656, 0.07727756, 0.06875796, 0.05918984, 0.04880096, 0.03786904, 0.02676332, 0.01592385, 0.00647946, 0.00370956, 0.00112514, -0.00046881, -0.00191488, -0.00329201, -0.00470585, -0.00688469, -0.00912202, -0.01720842, -0.02488211, -0.03226730, -0.03908459, -0.04503763, -0.04986836, -0.05338180, -0.05551392, -0.05636585, -0.05605816, -0.05472399, -0.05254383, -0.04969990, -0.04637175, -0.04264894, -0.03859653, -0.03433153, -0.02996944, -0.02560890, -0.02134397, -0.01726049, -0.01343567, -0.00993849, -0.00679919, -0.00402321, -0.00180118, -0.00044469, 0.00000000]

# my airfoil coordinates
xaf = [1.0, 0.99669031, 0.98707731, 0.97177624, 0.95120013, 0.92554805, 0.89513234, 0.86042844, 0.82194932, 0.78024605, 0.73589923, 0.68948289, 0.64155561, 0.59265387, 0.5432852, 0.49390561, 0.44491485, 0.39667271, 0.34949257, 0.30376269, 0.25995386, 0.21850111, 0.17980252, 0.14421934, 0.11207729, 0.08364666, 0.05914156, 0.03873701, 0.02260069, 0.01080931, 0.00329342, 0.00180844, 0.00078264, 0.00016559, 2.4e-6, 9.545e-5, 0.0005152, 0.00124174, 0.00277989, 0.00480487, 0.0111091, 0.02423919, 0.04202841, 0.0645839, 0.09186189, 0.1236919, 0.15984299, 0.20001108, 0.24382591, 0.29085133, 0.34057636, 0.39242047, 0.4457723, 0.49999633, 0.55444002, 0.60844101, 0.66133467, 0.71246137, 0.7611783, 0.80686551, 0.84892722, 0.8867981, 0.91995456, 0.9479178, 0.97026156, 0.98661057, 0.99662097, 1.0]
yaf = [0.0, 0.00049003, 0.00207033, 0.00475426, 0.00825965, 0.01241627, 0.01723294, 0.0226708, 0.02862516, 0.034964, 0.0415088, 0.0480415, 0.05432933, 0.06013032, 0.06520439, 0.06931828, 0.07229087, 0.07399667, 0.0744182, 0.07365704, 0.07178323, 0.06883168, 0.064851, 0.05990381, 0.05406639, 0.04742811, 0.04011465, 0.03230144, 0.02418269, 0.01593784, 0.00791465, 0.00553261, 0.00337629, 0.0013778, 0.00015663, -0.00093355, -0.00196476, -0.0030302, -0.00467095, -0.00634666, -0.0102234, -0.01556325, -0.02006888, -0.02353976, -0.02598343, -0.02742948, -0.02793392, -0.02758763, -0.026503, -0.02481832, -0.02268597, -0.02024403, -0.01760693, -0.01488444, -0.01217351, -0.00956485, -0.00713548, -0.0049519, -0.0030634, -0.00151205, -0.00032506, 0.00048717, 0.00093779, 0.00105862, 0.00091038, 0.00057426, 0.00018251, 0.0]

chord = 1.9
twist = 0.0*pi/180
paxis = 0.4750 / chord

uni = Material(37.00e9, 9.00e9, 9.00e9, 4.00e9, 4.00e9, 4.00e9, 0.28, 0.28, 0.28, 1.86e3)
double = Material(10.30e9, 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.83e3)
gelcoat = Material(1e1, 1e1, 1e1, 1.0, 1.0, 1.0, 0.30, 0.30, 0.30, 1.83e3)
nexus = Material(10.30e9, 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.664e3)
balsa = Material(0.01e9, 0.01e9, 0.01e9, 2e5, 2e5, 2e5, 0.30, 0.30, 0.30, 0.128e3)
mat = [uni, double, gelcoat, nexus, balsa]

t = 0.001
theta = 20*pi/180
layer = Layer(balsa, t, theta)

idx = [3, 4, 2]  # material index
t = [0.000381, 0.00051, 18*0.00053]
theta = [0, 0, 20]*pi/180
layup1 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2]
t = [0.000381, 0.00051, 33*0.00053]
theta = [0, 0, 20]*pi/180
layup2 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2, 1, 5, 1, 2]
t = [0.000381, 0.00051, 17*0.00053, 38*0.00053, 1*0.003125, 37*0.00053, 16*0.00053]
theta = [0, 0, 20, 30, 0, 30, 20]*pi/180
layup3 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2, 5, 2]
t = [0.000381, 0.00051, 17*0.00053, 0.003125, 16*0.00053]
theta = [0, 0, 20, 0, 0]*pi/180
layup4 = Layer.(mat[idx], t, theta)

segments = [layup1, layup2, layup3, layup4]

xbreak = [0.0, 0.0041, 0.1147, 0.5366, 1.0]

idx = [1, 5, 1]
t = [38*0.00053, 0.003125, 38*0.00053]
theta = [0, 0, 0]*pi/180
web = Layer.(mat[idx], t, theta)

webs = [web, web]

webloc = [0.15, 0.5]
nodes, elements = afmesh(xaf, yaf, chord, twist, paxis, xbreak, webloc, segments, webs)

S, sc, tc = compliance_matrix(nodes, elements)
M, mc = mass_matrix(nodes, elements)

I traced it to line 492 of afmesh.jl. There, yiu and yil are searched to find their intersection. Because the nodes don't match up exactly (I'm presuming this is the reason), the values are not compared to each other but rather to the value of zero. In the example, this works great. Below I've plotted yiu (blue) and yil (red), and sure enough they intersect when both surfaces cross zero.

But in my example, this isn't the case (see below). In fact, my second point in yil happens to be barely positive. An easy workaround for that is to limit the possible range to the back half or so of the surface. But even then, the intersection doesn't match up with where they cross zero...

I have an idea on how to get around this - I don't know if this would impact the rest of the code much, but I think it does what the current version is trying to do. The idea is to search backwards from the TE, as it appears that in both cases the nodes match up with each other (i.e. node 40 on top is approximately lined up with node 40 on the bottom, etc.), so instead of comparing both sets of nodes against zero you could compare to each other, such as

# index from the end for both top and bottom 
idx_intersection_from_TE = findlast(xt->xt<0, (reverse(yiu[end-20:end]) .- reverse(yil[end-20:end]))) 
# Note: need to specify a specific range so the number of elements is the same - instead of 20 it could 
# be half of the total length of one of them, I don't think it really matters as long as it's long enough 
# to include the intersection

# since they have different lengths, compute the index from the start separately for both surfaces
iu = length(yiu) - idx_intersection_from_TE + 1
il = length(yil) - idx_intersection_from_TE + 1

This might be overkill or perhaps not robust to every case, but it at least works for these two scenarios.

Thoughts @andrewning?

time_domain_analysis does not seem to work with point_masses

Hi, Taylor.
I'm trying to do my time domain analyses using point masses now instead of the inertia matrix, as the point masses improved my modal analyses (as per https://discourse.julialang.org/t/ann-gxbeam-jl-v0-3-0-geometrically-exact-beam-theory/71337/11).
However, time time_domain_analysis doesn't seem to converge with point masses. Here's an example of a case that runs just fine without point masses, and does not converge with them:

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),
							 length(nodes) => GXBeam.PrescribedConditions(Fz=30, My=-0.2));
t=0:.05:.05
system,history,converged = time_domain_analysis(assembly, t;
                                                prescribed_conditions = prescribed_conditions);

Here converged==true.

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),
							 length(nodes) => GXBeam.PrescribedConditions(Fz=30, My=-0.2));
t=0:.05:.05
system,history,converged = time_domain_analysis(assembly, t;
                                                prescribed_conditions = prescribed_conditions,
                                                point_masses = pointmass);

Here converged==false.

Any help would be much appreciated!