Initial conditions

[luizpancini]

Hi Taylor,

I can't get this simulation of a simply-supported beam initially displaced to excite only its 2nd mode to run:

using GXBeam, LinearAlgebra

# Simply-supported beam subject to initial conditions

# beam properties
L = 1
EI = 1
I = 1
A = 1
ρ = 1

# create points
nelem = 20
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]
lengths, points, midpoints, Cab = discretize_beam(L,[0, 0, 0],nelem)

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

# compliance matrix for each beam element
compliance = fill(Diagonal([0, 0, 0, 0, 1/EI, 0]), nelem)

# mass matrix for each beam element
mass = fill(Diagonal([ρ*A, ρ*A, ρ*A, 2*ρ*I, ρ*I, ρ*I]), nelem)

# create assembly 
assembly = Assembly(points, start, stop; compliance=compliance, mass=mass)

# prescribed conditions
prescribed_conditions = (t) -> begin
    Dict(
        # simply supported left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_z=0),
        # simply supported right side
        nelem+1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_z=0) 
    )
end

# initial conditions
delta = 1e-3    # small number to stay in the linear regime
u0 = [[0,0,delta*sin(2*pi*midpoints[i][1]/L)] for i = 1:length(midpoints)]
theta0 = [[0,-delta*2*pi/L*cos(2*pi*midpoints[i][1]/L),0] for i = 1:length(midpoints)]

# simulation time
t = 0:1e-3:0.25

# solve
system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions = prescribed_conditions,
    u0=u0,
    initialize=true,
    iterations=50) #theta0=theta0

# Get displacements at 1/4 of the beam's length
point = nelem/4+1
field = [:u]
direction = [3]
w14_of_t = [getproperty(state.points[point[1]], field[1])[direction[1]] for state in history]

# Plot
using Plots
pyplot()
plot(xlabel = "Time (s)",
     ylabel = "\$w\$ (\$m\$)",
     grid = false,
     overwrite_figure=false
     )   
plot!(t, w14_of_t, label="")     
plot!(show=true)

# analytical solution
x_quarter = x[nelem/4]
omega2 = (2*pi/L)^2*sqrt(EI/(ρ*A);
w14_analytic = delta*cos(omega2*t)*sin(2*pi*x_quarter/L)

The initial conditions solver does not converge (reduced to 50 iterations or it just takes forever), for this one or any other type of initial conditions. As you can see, the initial displacements were provided as inputs, but setting both displacements and angles does not help either.

I managed to solve the problem in a separate code, by setting the initial conditions as prescribed conditions (boundary conditions) and solving the resulting problem as a static one. The solution states vector is then copied as a consistent set of initial states. The displacement solution should be like this:

I was just wondering if I'm doing something wrong with the inputs to the time domain analysis, and say that I'd be glad to help fix the issue (if there really is one).

[taylormcd]

I took a look and came up with the following solution. Note that the initialization procedure doesn't converge. In this case however, this is not a problem because the initial state and rate variables that are actually used do converge. The initialization procedure should probably be modified to account for these cases.

using GXBeam, LinearAlgebra

# simply-supported beam with excited second bending mode 

# beam properties
L = 1
EA = 1
EI = 1
ρA = 1

# generate points
nelem = 40
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

# define connectivity
start = 1:nelem
stop = 2:nelem+1

# compliance matrix for each beam element
compliance = fill(Diagonal([0.0, 0.0, 0.0, 0.0, 1/(EI), 0.0]), nelem)

# mass matrix for each beam element
mass = fill(Diagonal([ρA, ρA, ρA, 0.0, 0.0, 0.0]), nelem)

# create assembly 
assembly = Assembly(points, start, stop; compliance=compliance, mass=mass)

# excite the second bending mode
delta = 1e-3 # displacement magnitude
uz =  delta*sin.(2*pi*x/L) # linear displacement

# prescribe the vertical displacement of each point
prescribed_conditions = Dict(i => PrescribedConditions(uz = uz[i]) for i = 1:length(points))

# simply supported left side
prescribed_conditions[1] = PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_z=0)

# simply supported right side
prescribed_conditions[nelem+1] = PrescribedConditions(uz=0)

# solve for static operating conditions
system, converged = static_analysis(assembly; prescribed_conditions)

# postprocess results
state = AssemblyState(system, assembly; prescribed_conditions)

# extract initial conditions from the state vector
u0 = getproperty.(state.elements, :u)
theta0 = getproperty.(state.elements, :theta)

# set new prescribed conditions
prescribed_conditions = Dict( 
    # simply supported left side
    1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_z=0),
    # simply supported right side
    nelem+1 => PrescribedConditions(uz=0)
)

# excited mode natural frequency
ω = (2*pi/L)^2*sqrt((EI)/(ρA))

# solution time vector
t = range(0, 2*pi/ω, step=0.001)

# initialize state variables (note that this doesn't converge)
system, converged = initial_condition_analysis(assembly, t[1]; prescribed_conditions,
    u0=u0, theta0=theta0)

# perform time domain analysis
system, history, converged = time_domain_analysis!(system, assembly, t; prescribed_conditions, 
    initialize = false, reset_state=false)

# write visualization file
write_vtk("excited", assembly, history, t; scaling = 100)

# Get displacements at 1/4 of the beam's length
ipoint = div(nelem, 4) + 1

# calculated deflection
w14_gxbeam = [state.points[ipoint].u[3]/L for state in history]

# analytic deflection
w14_analytic = delta/L*cos.(ω*t)*sin(2*pi*x[ipoint]/L)

# Plot
using Plots
pyplot()
plot(
    title="Deflection at x/L = $(x[ipoint]/L)",
    xlabel = "Time (s)",
    ylabel = "\$w/L\$",
    grid = false,
    legend = :bottomleft,
    overwrite_figure=false
    )   
plot!(t, w14_analytic, label="Analytic")    
plot!(t, w14_gxbeam, label="GXBeam")     
plot!(show=true)

[luizpancini]

Ok, that seems to be the same procedure I applied here to work around the problem!