byuflowlab / GXBeam.jl

Pure Julia Implementation of Geometrically Exact Beam Theory
MIT License
87 stars 18 forks source link

time_domain_analysis does not seem to work with point_masses #71

Closed RibeiroAndre closed 2 years ago

RibeiroAndre commented 2 years ago

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!

taylormcd commented 2 years ago

It looks like the initial condition analysis doesn't converge. I'll have to look into that at some point. For now, you can just set initialize=false to skip the initialization.

RibeiroAndre commented 2 years ago

Thanks! As I'm doing FSI I'll need to wait for your next update, then, as I need to use initial conditions. Thanks again!