# Degenerate BEAM element (von Kármán type nonlinearity, moderately large disp. & rotation) and Newton-Raphson iteration

Cantilever beam subjected to end force

## Calculation of stiffness matrices

kill(all);

Start the modules:

load("vect");
ratprint:false$

Give an accessible temporary directory:

maxima_tempdir:"C:\\Users\\PoloskeiTamasHunor\\Pictures\\Maxima"$

Define the vector of nodal displacements  7DOFs(!):

ue:[u_1,v_1,φ_1,u_12,u_2,v_2,φ_2]$

degenerate BEAM element with 7DOFs

Define the interpolation functions for bending (Hermite) and extension (hierarchical):

N1:1-3*xi^2+2*xi^3$
N2:l*(xi-2*xi^2+xi^3)$
N3:3*xi^2-2*xi^3$
N4:l*(-xi^2+xi^3)$
N5:1-xi$
N6:xi$
N56:4*xi-4*xi^2$

Define the vector of interpolation functions for bending:

Nv_t:matrix([0,N1,N2,0,0,N3,N4])$

Calculate the strain-displacement vectors:

Bv_t:1/l*diff(Nv_t,xi);

BBv_t:1/l^2*diff(Nv_t,xi,2);

...ind the transposed vectors:

Bv:transpose(Bv_t);

BBv:transpose(BBv_t);

Define the vector of interpolation functions for extension:

Nu_t:matrix([N5,0,0,N56,N6,0,0]);

Bu_t:1/l*diff(Nu_t,xi);

Bu:transpose(Bu_t);

Calculation of the stress resultants, N - normal force, M - bending moment:

N:A*E11*(Bu_t.ue+1/2*(Bv_t.ue)^2);
M:-Iz*E11*(BBv_t.ue);

Calculation of the stiffness matrices.
Linear stiffness matrix for tension (LINK1D with quadratic interpolation)

K11:A*E11*l*integrate(Bu.Bu_t,xi,0,1);

Nonlinear coupling stiffness matrix (between bending and tension):

K12:A*E11*l*integrate((Bv_t.ue)*(Bu.Bv_t + Bv.Bu_t),xi,0,1);

Linear stiffness matrix for bending (BEAM1D - cubic interpolation):

K22a:Iz*E11*l*integrate(BBv.BBv_t,xi,0,1);

Nonlinear stiffness matrices, K = K(U) - quadratic dependence on the DOFs:

K22b:A*E11*l*integrate((Bv_t.ue)^2*(Bv.Bv_t),xi,0,1);

K22c:l*integrate(N*Bv.Bv_t,xi,0,1);

Total stiffness matrix for bending:

K22:K22a+K22b+K22c$

Creation of the total nonlinear element stiffness matrix (tension+bending):

Ke_NL:K11+K12+K22$

Creation of the linear element stiffness matrix (tension+bending):

Ke_L: K11+K22a;

## Calculation of the internal force vector

Calculation of the vector of internal forces, axial forces:

Ru_int:l*integrate(N*Bu,xi,0,1),ratsimp;

Calculation of the vector of transverse internal forces:

Rv_int:l*integrate( N*(Bv_t.ue)*Bv - M*BBv ,xi,0,1);

## The assembly process, structural matrices and vectors

Definition of the number of elements:

Ne:3$

Define the data of the problem:

data:[h=10e-3, b=40e-3, E11=60e9, L=1, A=200e-5, Iz=40e-3*(10e-3)^3/12, F=120, l=1/Ne]$

Evaluate the stiffness matrices and vector of forces:

Ke_L: ev(K11+K22a),data$

Ke_NL:ev(K11+K12+K22),data$

Re_int:ev(Ru_int+Rv_int),data$

Definition of the structural matrix dimension:

elemlist:makelist(makelist(4*(n-1)+k,k,1,7,1),n,1,Ne,1);

Dim:4*Ne+3$

Definition of the structural vector of external forces:

F_cond:zeromatrix(Dim-3,1)$
F_cond[Dim-3-1,1]:-ev(F),data$
F_cond;

load("mmpack.wxm")$

Calculation of the nonlinear element matrices with structural size:

for n: 1 thru Ne do (
    uList:makelist(u[i],i,elemlist[n]),
    Ke_NL_subs:subst(map("=",ue,uList),   Ke_NL),
    KNL[n]:ExtMatrix(Ke_NL_subs,[Dim,Dim],elemlist[n],elemlist[n])
);

Assembly of the nonlinear structural matrices by using the nonlinear element matrices:

KSNL:sum(KNL [n], n,1,Ne)$

Application of the boundary conditions (clamped-free - see Figure 1) - condensed stiffness matrix:

freedof:makelist(i,i,4,Dim);

KNL_cond:SubMatrix(KSNL,freedof,freedof)$

## Algorithm for the classical Newton-Raphson scheme

sketch of Newton-Raphson iteration for a degressive system

Number of load steps:

N_ls:6$

Number of iterations within a single load step:

N_it:8$

Newton-Raphson iteration procedure:
    - creation of the system of equations for each load step and iteration step
    - solution of the system of equations with the updated tangent stiffness matrix and internal force
    vector
    - recalculation of the internal force vector with the actual displacements for each element
    - assembly of the element force vectors
    - initialize the displacement and internal force vectors for the next step

load("lapack")$

U[1,0]:zeromatrix(Dim,1)$
R[1,0]:zeromatrix(Dim,1)$

http://maxima.sourceforge.net/docs/manual/maxima_85.html

for p: 1 thru N_ls do(
    
    print("LS: ", p),
    F_cond_ss:p*F_cond/N_ls,
    
    for k: 1 thru N_it do(
        
        R_cond:SubMatrix(R[p,k-1], freedof,[1]),
        KNL_cond_ss:subst(makelist(u[s] = U[p,k-1][s,1] ,s,1,Dim), ev(KNL_cond)),
        residual:(F_cond_ss-R_cond),
        dU[p,k]:dgesv(KNL_cond_ss,residual),
        
        printf(true, "IT:~t~03d~t |~tF conv:~t~12,3e~t|~tU conv:~t~12,3e~%", k,sqrt(residual.residual),sqrt(dU[p,k].dU[p,k])),
        
        U[p,k]:U[p,k-1]+ExtMatrix(dU[p,k],[Dim,1],freedof,[1]),
        
        for n: 1 thru Ne do  (
            
            uList:transpose(SubMatrix(U[p,k],elemlist[n],[1]))[1],
            RNL_subs:subst(map("=",ue,uList), Re_int),
            RNL [n]:ExtMatrix(RNL_subs ,[Dim,1],elemlist[n],[1])
            
        ),
        
        R[p,k]: sum(RNL [ n], n, 1, Ne)
        
    ),
    /*print(U[p,N_it],R[p,N_it]),*/
    print("======="),
    U[p+1,0]:U[p,N_it],
    R[p+1,0]:R[p,N_it]
);

dU[2,0]$
R[2,0]$

## Plot of load-displacement response

Create the plot of the nonlinear FE solution, where "dd" sets the position of DOF in the
displacement vector:

dd:1$

Nonlinear FE solution for the beam end displacement in [m]:

Response: makelist([abs(U[q,0][Dim-dd,1]),abs(R[q,0][Dim-1,1])],q,1,N_ls+1)$

Analytical linear elastic solution and its plot:

v_an:ev(-F*(3*L*x^2-x^3)/(6*Iz*E11)),data;

Response_lin:[[0,0],[  at(subst(x=L, abs(v_an)),data)   , at(F,data)   ]]$

Display the load displacement response of the endpoint of the beam:

wxplot2d([ 
[discrete,Response_lin]   ,
[discrete,Response]], 
[x,0,0.3],[y,0,120],
grid2d,
[xlabel,"Deflection at beam end [m]"],
[ylabel,"External load [N]"],
[legend,"Linear solution","Nonlinear solution","","","",""],
[color,blue,red,black,magenta])$

Display the load rotation response of the endpoint of the beam:

dd:0$

Response_phi: makelist([abs(U[q,0][Dim-dd,1]),abs(R[q,0][Dim-1,1])],q,1,N_ls+1)$

Analytical linear elastic solution for the rotation and its plot:

Response_lin_phi:[[0,0],[  at(subst(x=L, abs(diff(v_an,x))),data)   , at(F,data)   ]]$

Display the load rotation response of the endpoint of the beam:

wxplot2d([ 
[discrete,  Response_lin_phi]   ,
[discrete,Response_phi]], 
[x,0,0.3],[y,0,120],
grid2d,
[xlabel,"Rotation at beam end [rad]"],
[ylabel,"External load [N]"],
[legend,"Linear solution","Nonlinear solution","","","",""],
[color,blue,red])$

## Plots of deflection and axial displacement functions

Calculation of the axial displacement function for each element:

for n: 1 thru Ne do(
    u0[n]:subst(xi=x/l-(n-1), Nu_t.SubMatrix(U[N_ls+1,0],elemlist[n],[1]))
);

Display the axial displacement function

wxdraw2d(
    title = "u0(x) - axial displacement function",
    grid = true,
    color=red,
    makelist(parametric( x, at(u0[n],data) , x , (n-1)*at(l,data) , n*at(l,data) ) , n , 1 , Ne ) );

Calculation of the deflection function for each element:

for i: 1 thru Ne do(
    v[i]:subst(xi=x/l-(i-1), Nv_t.SubMatrix(U[N_ls+1,0],elemlist[i],[1]))
);

Display the deflection functions from linear and nonlinear solutions:

wxdraw2d(
    title = "v(x) - transverse deflection",
    grid = true,
    key   = "Nonlinear solution",
    color=red,
    makelist(parametric( x, at(v[i],data) , x , (i-1)*at(l,data) , i*at(l,data) ) , i , 1 , Ne ),
    key   = "Linear solution",
    color=blue,
    parametric (x, v_an , x, 0, at(L,data)) );