"Free-flight" boundary conditions

[luizpancini]

Hi,

Got another issue... it concerns problems with "free-flight" boundary conditions, say for a beam free to move in 3D space.

First, I verified what seems to be an issue even in a 2D constrained flight. The problem is the "flying spaghetti" from [^1], here is the GXBeam input code:

using GXBeam, LinearAlgebra

# SIMO & VU-QUOC's 2D flying spaghetti

# beam properties
L = 10; 
EA = 1e4; 
GA = EA; 
EI = 500; 
GJ = EI; 
rhoA = 1; 
rhoIs = 20; 
rhoI = 10;

# Number of elements
nelem = 20

# Frame for initial rotation of the beam
Cab = [0.6 0.0 0.8
       0.0 1.0 0.0
      -0.8 0.0 0.6]

# Discretization variables      
lengths, points, midpoints, Cab = discretize_beam(L, [0, 0, 0], nelem, frame=Cab)

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

# stiffness matrix for each beam element
stiffness = fill(Diagonal([EA, GA, GA, GJ, EI, EI]), nelem)

# mass matrix for each beam element
mass = fill(Diagonal([rhoA, rhoA, rhoA, rhoIs, rhoI, rhoI]), nelem)

# create assembly
assembly = Assembly(points, start, stop;
    stiffness = stiffness,
    mass = mass,
    frames = Cab,
    lengths = lengths,
    midpoints = midpoints)

 # external forces
 tau = 2.5
M = (t) -> begin
    if 0.0 <= t <= tau
        80
    else
        zero(t)
    end
end

# prescribed conditions
prescribed_conditions = (t) -> begin
    Dict(
        # left side: no out of plane movement
        1 => PrescribedConditions(uy=0, theta_z=0),
        # right side: no out of plane movement, applied forces 
        nelem+1 => PrescribedConditions(Fx=M(t)/10, My=M(t), uy=0, theta_z=0) 
    )
end

# simulation time
t = 0:1e-2:13

# analysis
system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions = prescribed_conditions)

# Get displacement outputs at both ends of the beam    
point = [1 nelem+1]
field = [:u, :u]
direction = [1, 3]
ylabel = ["\$u\$ (\$m\$)","\$w\$ (\$m\$)"]

u1 = [getproperty(state.points[point[1]], field[1])[direction[1]] for state in history]
w1 = [getproperty(state.points[point[1]], field[2])[direction[2]] for state in history]
uend = [getproperty(state.points[point[2]], field[1])[direction[1]] for state in history]
wend = [getproperty(state.points[point[2]], field[2])[direction[2]] for state in history]

# plot results
using Plots
pyplot()

# plot horizontal displacement
plot(xlabel = "Time (s)",
     ylabel = ylabel[1],
     grid = false,
     overwrite_figure=false
     )   
plot!(t, u1, label="Node 1")     
plot!(t, uend, label="Node end")
plot!(show=true)

# plot vertical displacement
plot(xlabel = "Time (s)",
     ylabel = ylabel[2],
     grid = false,
     overwrite_figure=false
     )   
plot!(t, w1, label="Node 1")     
plot!(t, wend, label="Node end")
plot!(show=true)

A motion plot of the problem is as follows:

As you can see, after about 7 seconds, the beam starts to exhibit an unexpected motion, moving towards the left and downwards, even though its momentum should make it keep going to the right indefinitely, with no vertical movement. The authors in the original paper[^1] do not make it clear if they modified their motion plots in any way, but a more recent work[^3] confirms that indeed the beam keeps moving to the right:

Setting the prescribed condition u_z = 0 for the beam's midpoint fixes the problem, but such boundary condition is not stated by either [^1] or [^3], neither should it be necessary. Comparing solutions, the deformed shapes of the beam seem to agree, but the positions in space do not.

I also investigated the 3D version of that problem[^2], in which an out-of-plane moment is also applied, a true free-flight case. Here is the GXBeam file:

using GXBeam, LinearAlgebra

# SIMO & VU-QUOC's 3D flying spaghetti

# beam properties
L = 10; 
EA = 1e4; 
GA = EA; 
EI = 500; 
GJ = EI; 
rhoA = 1; 
rhoIs = 20; 
rhoI = 10;

# Number of elements
nelem = 20

# Frame for initial rotation of the beam
Cab = [0.6 0.0 0.8
       0.0 1.0 0.0
      -0.8 0.0 0.6]

# Discretization variables      
lengths, points, midpoints, Cab = discretize_beam(L, [0, 0, 0], nelem, frame=Cab)

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

# stiffness matrix for each beam element
stiffness = fill(Diagonal([EA, GA, GA, GJ, EI, EI]), nelem)

# mass matrix for each beam element
mass = fill(Diagonal([rhoA, rhoA, rhoA, rhoIs, rhoI, rhoI]), nelem)

# create assembly
assembly = Assembly(points, start, stop;
    stiffness = stiffness,
    mass = mass,
    frames = Cab,
    lengths = lengths,
    midpoints = midpoints)

 # forces
 M0 = 200
 tau = 2.5
M = (t) -> begin
    if 0.0 <= t <= tau
        M0*t/tau
    elseif tau < t <= 2*tau
        2*M0*(1.0-t/(2*tau))
    else
        zero(t)
    end
end

# prescribed conditions
prescribed_conditions = (t) -> begin
    Dict(
        # forces on right side
        nelem+1 => PrescribedConditions(Fx=M(t)/10, My=M(t), Mz=M(t)/2) 
    )
end

# simulation time (problems get serious after about 3 seconds)
t = 0:1e-3:3.14

# analysis
system, history, converged = time_domain_analysis(assembly, t;
    prescribed_conditions = prescribed_conditions)

# Get displacement outputs at both ends of the beam      
point = [1 nelem+1]
field = [:u, :u, :theta]
direction = [1, 3, 2]
ylabel = ["\$u\$ (\$m\$)","\$w\$ (\$m\$)","Rodrigues parameter \$\\theta_y\$"]

u1 = [getproperty(state.points[point[1]], field[1])[direction[1]] for state in history]
w1 = [getproperty(state.points[point[1]], field[2])[direction[2]] for state in history]
uend = [getproperty(state.points[point[2]], field[1])[direction[1]] for state in history]
wend = [getproperty(state.points[point[2]], field[2])[direction[2]] for state in history]
thetay_end = [getproperty(state.points[point[2]], field[3])[direction[3]] for state in history]

# plot results
using Plots
pyplot()

# plot u
plot(xlabel = "Time (s)",
     ylabel = ylabel[1],
     grid = true,
     overwrite_figure=false
     )   
plot!(t, u1, label="Node 1")     
plot!(t, uend, label="Node end")
plot!(show=true)

# plot w
plot(xlabel = "Time (s)",
     ylabel = ylabel[2],
     grid = true,
     overwrite_figure=false
     )  
plot!(t, w1, label="Node 1")     
plot!(t, wend, label="Node end")
plot!(show=true)

# plot theta_y
plot(xlabel = "Time (s)",
     ylabel = ylabel[3],
     grid = true,
     overwrite_figure=false
     )  
plot!(t, thetay_end, label="")     
plot!(show=true)

The plot of the theta_y rotation parameter reveals an unstable behavior of the beam, and the solution only converges up until around 3.16 seconds:

No discretization or time step refinement could remedy the numerical instability. A plane view motion plot makes it clear that there is an issue in the simulation:

Since the program appears to work perfectly in cases where the structure is restrained in some way, I was led to believe the issue only happens in these free-flight cases. Any ideas of what may be going on? Maybe there is a need to integrate the "a" frame states over time as well?

[^1]: Simo, J. C., & Vu-Quoc, L. (1986). On the dynamics of flexible beams under large overall motions—The plane case: Part II. Journal of Applied Mechanics, 53, 855-863. [^2]: Simo, J. C., & Vu-Quoc, L. (1988). On the dynamics in space of rods undergoing large motions—a geometrically exact approach. Computer Methods in Applied Mechanics and Engineering, 66(2), 125-161. [^3]: Palacios, R., & Cea, A. (2019). Nonlinear modal condensation of large finite element models: Application of Hodges’s intrinsic theory. AIAA Journal, 57(10), 4255-4268.

[taylormcd]

Isn't it usually better to decouple the rigid body motion and the structural motion? That was just my understanding, but perhaps that's a non-issue here. For the first case I might need to revisit the rotation rescaling I use for large angles to make sure it's correct. For the second case I think adding a small amount of structural damping might alleviate the problem (currently there is no damping). I'm working on adding that now.

On Fri, Mar 25, 2022, 11:19 AM Luiz Guilherme Pancini dos Santos < @.***> wrote:

Hi,

Got another issue... it concerns problems with "free-flight" boundary conditions, say for a beam free to move in 3D space.

First, I verified what seems to be an issue even in a 2D constrained flight. The problem is the "flying spaghetti" from 1 <#m_7144067857463109944_user-content-fn-1-5f385703876c0c5f5dd529f2b994984d>, here is the GXBeam input code:

using GXBeam, LinearAlgebra

SIMO & VU-QUOC's 2D flying spaghetti

beam properties

L = 10;

EA = 1e4;

GA = EA;

EI = 500;

GJ = EI;

rhoA = 1;

rhoIs = 20;

rhoI = 10;

Number of elements

nelem = 20

Frame for initial rotation of the beam

Cab = [0.6 0.0 0.8

   0.0 1.0 0.0

  -0.8 0.0 0.6]

Discretization variables

lengths, points, midpoints, Cab = discretize_beam(L, [0, 0, 0], nelem, frame=Cab)

index of endpoints of each beam element

start = 1:nelem

stop = 2:nelem+1

stiffness matrix for each beam element

stiffness = fill(Diagonal([EA, GA, GA, GJ, EI, EI]), nelem)

mass matrix for each beam element

mass = fill(Diagonal([rhoA, rhoA, rhoA, rhoIs, rhoI, rhoI]), nelem)

create assembly

assembly = Assembly(points, start, stop;

stiffness = stiffness,

mass = mass,

frames = Cab,

lengths = lengths,

midpoints = midpoints)

external forces

tau = 2.5

M = (t) -> begin

if 0.0 <= t <= tau

    80

else

    zero(t)

end

end

prescribed conditions

prescribed_conditions = (t) -> begin

Dict(

    # left side: no out of plane movement

    1 => PrescribedConditions(uy=0, theta_z=0),

    # right side: no out of plane movement, applied forces

    nelem+1 => PrescribedConditions(Fx=M(t)/10, My=M(t), uy=0, theta_z=0)

)

end

simulation time

t = 0:1e-2:13

analysis

system, history, converged = time_domain_analysis(assembly, t;

prescribed_conditions = prescribed_conditions)

Get displacement outputs at both ends of the beam

point = [1 nelem+1]

field = [:u, :u]

direction = [1, 3]

ylabel = ["$u$ ($m$)","$w$ ($m$)"]

u1 = [getproperty(state.points[point[1]], field[1])[direction[1]] for state in history]

w1 = [getproperty(state.points[point[1]], field[2])[direction[2]] for state in history]

uend = [getproperty(state.points[point[2]], field[1])[direction[1]] for state in history]

wend = [getproperty(state.points[point[2]], field[2])[direction[2]] for state in history]

plot results

using Plots pyplot()

plot horizontal displacement

plot(xlabel = "Time (s)",

 ylabel = ylabel[1],

 grid = false,

 overwrite_figure=false

 )

plot!(t, u1, label="Node 1") plot!(t, uend, label="Node end") plot!(show=true)

plot vertical displacement

plot(xlabel = "Time (s)",

 ylabel = ylabel[2],

 grid = false,

 overwrite_figure=false

 )

plot!(t, w1, label="Node 1") plot!(t, wend, label="Node end") plot!(show=true)

A motion plot of the problem is as follows:

[image: flexible_beam_2D] https://user-images.githubusercontent.com/40968757/160163079-86b6a8b3-6997-43e8-a643-1a40532fde43.gif

As you can see, after about 7 seconds, the beam starts to exhibit an unexpected motion, moving towards the left and downwards, even though its momentum should make it keep going to the right indefinitely, with no vertical movement. The authors in the original paper1 <#m_7144067857463109944_user-content-fn-1-5f385703876c0c5f5dd529f2b994984d> do not make it clear if they modified their motion plots in any way, but a more recent work2 <#m_7144067857463109944_user-content-fn-3-5f385703876c0c5f5dd529f2b994984d> confirms that indeed the beam keeps moving to the right:

[image: Palacios_2D_spaghetti] https://user-images.githubusercontent.com/40968757/160164378-e6cc0c13-86f2-4237-b71c-3603d64edc49.png

Setting the prescribed condition u_z = 0 for the beam's midpoint fixes the problem, but such boundary condition is not stated by either 1 <#m_7144067857463109944_user-content-fn-1-5f385703876c0c5f5dd529f2b994984d> or 2 <#m_7144067857463109944_user-content-fn-3-5f385703876c0c5f5dd529f2b994984d>, neither should it be necessary. Comparing solutions, the deformed shapes of the beam seem to agree, but the positions in space do not.

I also investigated the 3D version of that problem3 <#m_7144067857463109944_user-content-fn-2-5f385703876c0c5f5dd529f2b994984d>, in which an out-of-plane moment is also applied, a true free-flight case. Here is the GXBeam file:

using GXBeam, LinearAlgebra

SIMO & VU-QUOC's 3D flying spaghetti

beam properties

L = 10;

EA = 1e4;

GA = EA;

EI = 500;

GJ = EI;

rhoA = 1;

rhoIs = 20;

rhoI = 10;

Number of elements

nelem = 20

Frame for initial rotation of the beam

Cab = [0.6 0.0 0.8

   0.0 1.0 0.0

  -0.8 0.0 0.6]

Discretization variables

lengths, points, midpoints, Cab = discretize_beam(L, [0, 0, 0], nelem, frame=Cab)

index of endpoints of each beam element

start = 1:nelem

stop = 2:nelem+1

stiffness matrix for each beam element

stiffness = fill(Diagonal([EA, GA, GA, GJ, EI, EI]), nelem)

mass matrix for each beam element

mass = fill(Diagonal([rhoA, rhoA, rhoA, rhoIs, rhoI, rhoI]), nelem)

create assembly

assembly = Assembly(points, start, stop;

stiffness = stiffness,

mass = mass,

frames = Cab,

lengths = lengths,

midpoints = midpoints)

forces

M0 = 200

tau = 2.5

M = (t) -> begin

if 0.0 <= t <= tau

    M0*t/tau

elseif tau < t <= 2*tau

    2*M0*(1.0-t/(2*tau))

else

    zero(t)

end

end

prescribed conditions

prescribed_conditions = (t) -> begin

Dict(

    # forces on right side

    nelem+1 => PrescribedConditions(Fx=M(t)/10, My=M(t), Mz=M(t)/2)

)

end

simulation time (problems get serious after about 3 seconds)

t = 0:1e-3:3.14

analysis

system, history, converged = time_domain_analysis(assembly, t;

prescribed_conditions = prescribed_conditions)

Get displacement outputs at both ends of the beam

point = [1 nelem+1]

field = [:u, :u, :theta]

direction = [1, 3, 2]

ylabel = ["$u$ ($m$)","$w$ ($m$)","Rodrigues parameter $\theta_y$"]

u1 = [getproperty(state.points[point[1]], field[1])[direction[1]] for state in history]

w1 = [getproperty(state.points[point[1]], field[2])[direction[2]] for state in history]

uend = [getproperty(state.points[point[2]], field[1])[direction[1]] for state in history]

wend = [getproperty(state.points[point[2]], field[2])[direction[2]] for state in history]

thetay_end = [getproperty(state.points[point[2]], field[3])[direction[3]] for state in history]

plot results

using Plots pyplot()

plot u

plot(xlabel = "Time (s)",

 ylabel = ylabel[1],

 grid = true,

 overwrite_figure=false

 )

plot!(t, u1, label="Node 1") plot!(t, uend, label="Node end") plot!(show=true)

plot w

plot(xlabel = "Time (s)",

 ylabel = ylabel[2],

 grid = true,

 overwrite_figure=false

 )

plot!(t, w1, label="Node 1") plot!(t, wend, label="Node end") plot!(show=true)

plot theta_y

plot(xlabel = "Time (s)",

 ylabel = ylabel[3],

 grid = true,

 overwrite_figure=false

 )

plot!(t, thetay_end, label="") plot!(show=true)

The plot of the theta_y rotation parameter reveals an unstable behavior of the beam, and the solution only converges up until around 3.16 seconds:

[image: 3D_spaghetti_thetay] https://user-images.githubusercontent.com/40968757/160165634-5858b8f0-fa7f-46eb-9548-b576a7501965.png

No discretization or time step refinement could remedy the numerical instability. A plane view motion plot makes it clear that there is an issue in the simulation:

[image: 3D_spaghetti] https://user-images.githubusercontent.com/40968757/160165826-89f6e3d9-b6cd-445d-b76f-23cbc1f4f559.gif

Since the program appears to work perfectly in cases where the structure is restrained in some way, I was led to believe the issue only happens in these free-flight cases. Any ideas of what may be going on? Footnotes

Simo, J. C., & Vu-Quoc, L. (1986). On the dynamics of flexible beams under large overall motions—The plane case: Part II. Journal of Applied Mechanics, 53, 855-863. ↩ <#m_7144067857463109944_user-content-fnref-1-5f385703876c0c5f5dd529f2b994984d> ↩2 <#m_7144067857463109944_user-content-fnref-1-2-5f385703876c0c5f5dd529f2b994984d> ↩3 <#m_7144067857463109944_user-content-fnref-1-3-5f385703876c0c5f5dd529f2b994984d> 2.

Palacios, R., & Cea, A. (2019). Nonlinear modal condensation of large finite element models: Application of Hodges’s intrinsic theory. AIAA Journal, 57(10), 4255-4268. ↩ <#m_7144067857463109944_user-content-fnref-3-5f385703876c0c5f5dd529f2b994984d> ↩2 <#m_7144067857463109944_user-content-fnref-3-2-5f385703876c0c5f5dd529f2b994984d> 3.

Simo, J. C., & Vu-Quoc, L. (1988). On the dynamics in space of rods undergoing large motions—a geometrically exact approach. Computer Methods in Applied Mechanics and Engineering, 66(2), 125-161. ↩ <#m_7144067857463109944_user-content-fnref-2-5f385703876c0c5f5dd529f2b994984d>

— Reply to this email directly, view it on GitHub https://github.com/byuflowlab/GXBeam.jl/issues/75, or unsubscribe https://github.com/notifications/unsubscribe-auth/AEDIQX7LVO3D4XZ3PXO43ZLVBXYQ7ANCNFSM5RUX36GQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[luizpancini]

My understanding is that they can't be decoupled if the structural deformation is large. For instance, when the beam deforms under loads and alters its shape, its mass moment of inertias about its center of gravity will change, thus changing its "rigid-body" dynamics. If the deformations are small, then these changes can be neglected and the motions can be decoupled.

I think the rescaling is fine, tested it for rotations over thousands of radians for a driven robot arm.

Structural damping will be a nice addition.

[luizpancini]

My guess is that the "a" body-frame states, x0, v0, ω0, a0, α0 need to be correctly updated at each time step by time integration of Euler's differential equations of motion (with e.g. a Runge-Kutta method). That would also require updating the inertia tensor of the entire beam assembly every time step, because of beam deformation. When the motion of the "a" frame is specified over time, like in constrained problems, this entire step can be skipped, but in free-flight conditions the values of these states actually influence the solution.

[taylormcd]

I currently treat x0, v0, ω0, a0, and α0 as inputs rather than state variables. What I meant to express in my earlier comment was that Euler's differential equations of motion probably need to be included in order for the rigid body motion to be modeled correctly (which is the same conclusion you came to). I should not have used the word "decoupled", since the rigid body motion definitely has to be coupled with the structural model in order for the motion to be modeled correctly.

[taylormcd]

This is a nice test case though for when I do add in the rigid body motion. Hopefully that will happen soon.

[luizpancini]

As an update on these problems, the bug was solved for the 2D case! I realized that the response of the beam was correct until it tumbled for a full rotation, when the rescaling of the rotation parameters needed to be done, which led me to believe it had to be root of the issue. Then, I have seen that you had corrected the rotation tensor and tangent operator derivatives on your new reformulated branch of the code, in order to account for the derivatives of the rotation parameters scaling w.r.t. the parameters themselves and time. Applying these changes to the matrices C_0, Q_0, Qinv_0, C_t and C_t_0 has solved the problem, as can be seen by comparing the new response:

with that of [^1]:

However, that has not been enough to solve the bug in the 3D case, it still crashes, even with the reformulated "damping" branch, so there is still something wrong to correct. Anyways, at least it indicates that integrating the body-fixed frame origin states will not be necessary, somehow Hodges' theory already accounts for it in the equations of motion.

[^1]: Simo, J. C., & Vu-Quoc, L. (1986). On the dynamics of flexible beams under large overall motions—The plane case: Part II. Journal of Applied Mechanics, 53, 855-863.