Solving the `tipmoment` example with DifferentialEquations.jl

[axla-io]

Hi and thanks for an interesting package! I'm interested in trying out the interface with DifferentialEquations.jl on the tipmoment example. My code is shown below. It doesn't converge. I was wondering if I'm making any obvious mistake here. I'm very new to GXBeam.jl

using GXBeam
using LinearAlgebra
using DifferentialEquations

L = 12 # inches
h = w = 1 # inches
E = 30e6 # lb/in^4 Young's Modulus

A = h * w
Iyy = w * h^3 / 12
Izz = w^3 * h / 12

# bending moment (applied at end)
# note that solutions for λ > 1.8 do not converge
λ = 0.4
m = pi * E * Iyy / L
M = λ * m

# create points
nelem = 16
x = range(0, L, length=nelem + 1)
y = zero(x)
z = zero(x)
points = [[x[i], y[i], z[i]] for i in axes(x, 1)]

# index of endpoints for each beam element
start = 1:nelem
stop = 2:nelem+1

# compliance matrix for each beam element
# ones on the diagonal so that the matrix can be inverted, needed for DifferentialEquations
compliance = fill(Diagonal([1 / (E * A), 1.0, 1.0, 1.0, 1 / (E * Iyy), 1 / (E * Izz)]), nelem)

# mass matrix for each beam element added for DifferentialEquations
ρ = 287 # Approximate density of steel

mass = fill(Diagonal([ρ * A, ρ * A, ρ * A, 2 * ρ * Iyy, ρ * Iyy, ρ * Iyy]), nelem)

# create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop, compliance=compliance, mass=mass)

tspan = (0.0, 1.0)

# create dictionary of prescribed conditions
prescribed_conditions = Dict(
    # fixed left side
    1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
    # moment on right side
    nelem + 1 => PrescribedConditions(Mz=M)
)

# define named tuple with parameters
p = (; prescribed_conditions)

# run initial condition analysis to get consistent set of initial conditions
system, converged = initial_condition_analysis(assembly, tspan[1];
    prescribed_conditions=prescribed_conditions,
    constant_mass_matrix=false,
    structural_damping=true)

# construct DAEProblem
prob = DAEProblem(system, assembly, tspan, p;
    structural_damping=true)

# solve DAEProblem
@time sol = solve(prob, DABDF2(), saveat=[tspan[2]])
sol.retcode

The solver does not converge, instead I get:

┌ Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=1.0e-6. Aborting. There is either an error in your model specification or the true solution is unstable.

[taylormcd]

You'll probably get better results if you start from a steady state solution and then gradually increase the loads.

[taylormcd]

One of the hard things about time domain simulations is getting the right initial conditions.

[axla-io]

Thanks! Is it possible to define a time dependent loading or do you mean to repeatedly solve for equilibrium with increasing loads each time?

[taylormcd]

It's totally possible to do a time dependent loading. Just pass in the prescribed conditions as a function of time. See the wind turbine blade example or the joined wing example.