# Flutter instability of Beck's column with nonconservative force
(created in wxMaxima 16.12.0)

## Calculation of the mass, stiffness, geometric stiffness
and load stiffness matrices

Restart, delete previous variables and functions:

kill(all)$

Load the draw package:

load(draw)$

Give an accessible temporary directory:

maxima_tempdir:"d:/maxima/";

Number of elements:

Ne:4$

Problem data:

data:[l=1,A=82.0e-6, Iz=200.0e-10,E11=205.0e9, ρ=7800, L=1/Ne];

Definition of the material stiffness matrix:

Ke:Iz*E11/L^3* matrix(
 [12  ,  6*L   , -12  ,  6*L   ], 
 [6*L ,  4*L^2 , -6*L ,  2*L^2 ], 
 [-12 , -6*L   ,  12  , -6*L   ], 
 [6*L ,  2*L^2 , -6*L ,  4*L^2 ] );

Definition of the lumped mass matrix:

Me_lu:ρ*A*L/2*matrix(
 [1,0,0,0], 
 [0,0,0,0], 
 [0,0,1,0], 
 [0,0,0,0]
);

Definition of the geometric stiffness matrix for Euler-buckling:

KGe:-F/30/L*matrix(
 [ 36  ,  3*L   , -36  ,  3*L   ], 
 [ 3*L ,  4*L^2 , -3*L , -L^2   ], 
 [-36  , -3*L   ,  36  , -3*L   ], 
 [ 3*L , -L^2   , -3*L ,  4*L^2 ] );

Definition of the load stiffness matrix for nonconservative (follower) force:

KLe:F*matrix(
 [0,-1,0,0], 
 [0, 0,0,0], 
 [0, 0,0,1], 
 [0, 0,0,0] );

## Matrix assembly, calculation of structural matrices

Definition of the structural matrix dimension:

Dim:2*(Ne+1)$

Index mapping, element - DOF relation

dof:makelist(makelist(2*(n-1)+i,i,1,4),n,1,Ne);

Declaration of the element matrices with structural size:

for n: 1 thru Ne do(
    K[n] : zeromatrix(Dim,Dim),
    M_lu[n] : zeromatrix(Dim,Dim),
    KG[n] : zeromatrix(Dim,Dim),
    KL[n] : zeromatrix(Dim,Dim));

Calculation of the element matrices with structural size:

for n: 1 thru Ne do
for i: 1 thru 4 do
for j: 1 thru 4 do (
        K[n]     [ dof[n][i] , dof[n][j]] : Ke[i,j],
        M_lu[n]  [ dof[n][i] , dof[n][j]] : Me_lu[i,j],
        KG[n]    [ dof[n][i] , dof[n][j]] : KGe[i,j],
        KL[n]    [ dof[n][i] , dof[n][j]] : KLe[i,j] );

KG[1]$ KG[2]$ KG[3]$ KG[4]$

Assembly of the structural matrices by summing the element matrices with
structural size:

KS: sum(K[n],n,1,Ne)$
MS_lu: sum(M_lu[n],n,1,Ne)$
KGS: sum(KG[n],n,1,Ne)$
KLS: sum(KL[n],n,1,Ne)$

Application of the boundary conditions (clamped-free) - condensed matrices:

K_cond:zeromatrix(Dim-2,Dim-2)$
M_lu_cond:zeromatrix(Dim-2,Dim-2)$
KG_cond:zeromatrix(Dim-2,Dim-2)$
KL_cond:zeromatrix(Dim-2,Dim-2)$

for i: 1 thru Dim-2 do
for j: 1 thru Dim-2 do (
          K_cond[i,j] : KS[i+2,j+2], 
          M_lu_cond[i,j] : MS_lu[i+2,j+2], 
          KG_cond[i,j] : KGS[i+2,j+2], 
          KL_cond[i,j] : KLS[i+2,j+2] );

K_cond$M_lu_cond$KG_cond$KL_cond$

## Equation of motion, characteristic equation, root am-
plitude plots

Frequency parameter:

p:sqrt(ν*Iz*E11/ρ/A/(Ne*L)^4);

Force parameter:

F:λ*Iz*E11/(L*Ne)^2;

The characteristic equation of the system is:

S_lu:p^2*M_lu_cond+K_cond+KG_cond+KL_cond,[p=p,F=F]$

Calculation of the determinant:

f_lu:determinant(S_lu)/Iz^(2*Ne)/E11^(2*Ne)*L^(4*Ne),data$

Load the implicit plot module:

load(implicit_plot)$
ratprint:false$

Display the system behavior in the λ−ν plane! Flutter is expected around λ = 20 and
110, however, the curves, where complex roots are expected cannot be displayed
by implicit plot...

implicit_plot([ev(f_lu,λ=λ,ν=ν)=0],[λ,0,120],[ν,-3000,0],grid2d,[legend,"lumped"]);

Calculation of the characteristic roots in the function of load control parameters
Step size for the numerical solution:

dλ:1$

Number of points in the λ axis:

Nmax:150$

Calculation of solutions (four roots) if λ = 0..150:

for n: 0 thru Nmax do(
P1_lu[n]:flatten(allroots(ev(f_lu,λ=dλ*n)=0)),
RP_lu:map(lambda([val],ev(ν,val)),P1_lu[n]),
RP1_lu[n]:RP_lu[1],
RP2_lu[n]:RP_lu[2],
RP3_lu[n]:RP_lu[3],
RP4_lu[n]:RP_lu[4] );

Plot the solution in the λ - sqrt(|ν|) parameter plane:

wxplot2d([
    [discrete,makelist([dλ*n,sqrt(abs(RP1_lu[n]))],n,0,Nmax)],
    [discrete,makelist([dλ*n,sqrt(abs(RP2_lu[n]))],n,0,Nmax)],
    [discrete,makelist([dλ*n,sqrt(abs(RP3_lu[n]))],n,0,Nmax)],
    [discrete,makelist([dλ*n,sqrt(abs(RP4_lu[n]))],n,0,Nmax)] ],
[x,0,120],grid2d,[xlabel,"λ"],[ylabel,"|ν|"],[legend,"lumped","","","","","",""],
[style,[points,.5]],[point_type,bullet],[color,blue,blue,blue,blue])$

## Calculation of the mode shapes

Define the size of the vector of nodal displacements:

UU:zeromatrix(Dim-2,1);

Fill up the vector with the nodal displacement parameters:

for i: 1 thru Dim-2 do
UU[i,1]:u[i];

We set the last nodal transverse displacement to be equal to 1 mm:

UU[Dim-3,1]:1e-3;

The vector of nodal displacements is:

UU;

The characteristic equation of the system is:

char_poly_lu:ratsimp(S_lu.UU)$

Solution of the system of linear equations for the nodal displacements:

uisol_lu:linsolve(makelist(char_poly_lu[i,1]=0,i,1,Dim-3),
 append(makelist(u[i],i,1,Dim-4),[u[Dim-2]]))$

Declaration of the structural vector of nodal displacements:

U_lu:zeromatrix(Dim,1);

Fill up the vector with the nodal displacement parameters:

for i: 3 thru Dim-1 do
( U_lu[i,1]:ev(u[i-2],uisol_lu[i-2]));
  

U_lu[Dim-1,1]:1E-3$
  U_lu[Dim,1]:ev(u[Dim-2],uisol_lu[Dim-3])$

Definition of the Hermitian interpolation functions:

Nv:[ 1-3*xi^2+2*xi^3 , L*xi-2*L*xi^2+L*xi^3 , 3*xi^2-2*xi^3 , -L*xi^2+L*xi^3 ],data;

Calculation of the deflection function for each element:

for n: 1 thru Ne do
(
    v_lu[n]:first(at(at(makelist(U_lu[i],i,dof[n]).Nv,xi=x/L-(n-1)),data))
)$

Calculate the ν parameter (proportional to the square of char. exponent) at λ = 0:

q0_1:0.0$
R_1_lu:sort(allroots(at(f_lu,λ=q0_1)),ordergreatp);

Substitute the parameters back and calculate the first four mode shapes:

wxdraw2d(
    title = "q0 = 0, Lumped",
    grid = true,
    color=black,
    makelist(parametric(x,at(v_lu[i],[λ=q0_1,R_1_lu[1]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
    color=blue,
    makelist(parametric(x,at(v_lu[i],[λ=q0_1,R_1_lu[2]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
    color=green,
    makelist(parametric(x,at(v_lu[i],[λ=q0_1,R_1_lu[3]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
    color=red,
    makelist(parametric(x,at(v_lu[i],[λ=q0_1,R_1_lu[4]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne)
)$

Calculate the ν parameters close to the first flutter instability point:

q0_2_lu:18.3$
R_2_lu:sort(allroots(at(f_lu,λ=q0_2_lu)),ordergreatp);

Substitute the parameters back and calculate the first four mode shapes:

wxdraw2d(
title="q0 = 18.3, Lumped",
grid = true,
color=black,
makelist(parametric(x,at(v_lu[i],[λ=q0_2_lu,R_2_lu[1]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=blue,
makelist(parametric(x,at(v_lu[i],[λ=q0_2_lu,R_2_lu[2]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=green,
makelist(parametric(x,at(v_lu[i],[λ=q0_2_lu,R_2_lu[3]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=red,
makelist(parametric(x,at(v_lu[i],[λ=q0_2_lu,R_2_lu[4]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne)
)$

Calculate the ν parameters close to the second flutter instability point:

q0_3_lu:100$
R_3_lu:sort(allroots(at(f_lu,λ=q0_3_lu)),ordergreatp);

Substitute the parameters back and calculate the first four mode shapes:

wxdraw2d(
title="q0 = 100, Lumped",
grid = true,
color=black,
makelist(parametric(x,at(v_lu[i],[λ=q0_3_lu,R_3_lu[1]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=blue,
makelist(parametric(x,at(v_lu[i],[λ=q0_3_lu,R_3_lu[2]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=green,
makelist(parametric(x,at(v_lu[i],[λ=q0_3_lu,R_3_lu[3]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=red,
makelist(parametric(x,at(v_lu[i],[λ=q0_3_lu,R_3_lu[4]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne)
)$

Calculate the frequency parameters after the second flutter instability point:

q0_4_lu:110$
R_4_lu:sort(allroots(at(f_lu,λ=q0_4_lu)),ordergreatp);

Substitute the parameters back and calculate the first four mode shapes:

wxdraw2d(
title="q0 = 110, Lumped",
grid = true,
color=black,
makelist(parametric(x,at(v_lu[i],[λ=q0_4_lu,R_4_lu[1]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=blue,
makelist(parametric(x,at(v_lu[i],[λ=q0_4_lu,R_4_lu[2]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=green,
makelist(parametric(x,at(v_lu[i],[λ=q0_4_lu,R_4_lu[3]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne),
color=red,
makelist(parametric(x,at(v_lu[i],[λ=q0_4_lu,R_4_lu[4]]),x,at((i-1)*L,data),at(i*L,data)),i,1,Ne)
)$

## Calculation the time response of the system at flutter instability ( λ = 20)

Calculate the frequency parameters after the first flutter instability point:

R_lu:allroots(at(f_lu,λ=20));

Calculate the real and imaginary parts of the characteristic exponent:

eqa_lu:at(Iz*E11/(ρ*A*((Ne*L)^4))*realpart(at(ν,R_lu[1])),data)=(α^2-ω^2);
eqb_lu:at(Iz*E11/(ρ*A*((Ne*L)^4))*imagpart(at(ν,R_lu[1])),data)=2*α*ω;

Solve the system of equations:

R2_lu:solve([eqa_lu,eqb_lu],[α,ω]),numer;

Calculate the time response of the system:

T_lu(t):=exp(t*at(α,R2_lu[2][1]))*cos(t*at(ω,R2_lu[2][2]));

Plot the time function:

wxplot2d([T_lu(t)],[t,0,0.03],grid2d,[ylabel,"amplitude [m]"],[xlabel,"t[s]"],[legend,"lumped"])$

Plot the phase plane portrait (unstable focus):

ft_lu(t):=at(at(v_lu[Ne],[x=Ne*L,λ=20,R_lu[1]]),data)*T_lu(t)$

wxdraw2d(
grid=true,
nticks=500,
color=red,
key="lumped",
parametric(ft_lu(t),diff(ft_lu(t),t),t,0,10*%pi/abs(at(ω,R2_lu[2][2]))))$