9.1. Decoupled beam

9.1.1. Single straight beam

9.1.1.1. Axial dynamics

\[\begin{aligned} m \ddot{s}_z - \left( EA s'_z \right)' = f_z \end{aligned}\]

Piecewise first-degree polynomial approximation of the displacement. The base functions on the reference element \(\xi \in [-1,1]\) are

\[ \phi_1(\xi) = \frac{1}{2} \left( 1 - \xi \right) \quad , \quad \phi_2(\xi) = \frac{1}{2} \left( 1 + \xi \right) \]

and their first order derivatives read

\[ \phi'_1(\xi) = - \frac{1}{2} \quad , \quad \phi'_2(\xi) = \frac{1}{2} \]

Weak formulation.

\[\int_{\ell} w m \ddot{s} + \int_{\ell} w' EA s' = \int_{\ell} w f + \left.\left[ w EA s' \right]\right|_{\partial \ell}\]

Discrete problem. Exploiting additivity of the integrals over partition of domains,

\[\int_{\ell} = \sum_{k} \int_{\ell_k}\]

introducing the approximation \(s(x,t) = \sum_{j} \phi_j(x) s_j(t) = \boldsymbol\phi(x)^T \mathbf{s}(t)\), and testing over all the base functions \(\phi_i(x)\), the discrete problem reads

\[\int_{\ell} \phi_i(x) m \phi_j(x) \, dx \, \ddot{s}_j + \int_{\ell} \phi'_i(x) EA \phi'_j(x) \, dx \, s_j = \int_{\ell} \phi_i(x) f(x) \, dx + \left.\left[ \phi_i EA s' \right]\right|_{\partial \ell} \ .\]

In the \(i^{th}\) equation, non-zero contributions only come from elements \(k\) where both \(\phi_i\), \(\phi_j\) are simultaneously non-zero.

9.1.1.1.1. Example: clamped-free beam

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

#> Parameters
#> Physical parameters
l, m, EA = 1., 1., 1.
f = .0

#> Discretization parameters, for building mesh
nelems = 3
nnodes = nelems + 1

#> Boundary conditions
u0, F1 = .0, 1.   # Dirichlet (displ) in P0, Neumann (force) in P1
i_dir = [ 0 ]
u_dir = [ u0 ]
i_neu = [ nnodes-1 ]
f_neu = [ F1 ]

#> Mesh
#> Node coordinates
rr = np.linspace(0, 1, nnodes) * l
#> Node-element connectivity matrix
ee = np.array([ [ iel, iel+1 ] for iel in np.arange(nelems) ])

print(rr)
print(ee)

[0.         0.33333333 0.66666667 1.        ]
[[0 1]
 [1 2]
 [2 3]]

\[\begin{split}\begin{aligned} M_{11}^{loc} = M_{22}^{loc} & = \int_{\xi=-1}^{1} \left( \dfrac{1}{2} \left( 1- \xi \right)\right)^2 d \xi = \frac{1}{4} \left( 2 + \dfrac{2}{3} \right) = \dfrac{2}{3} \\ M_{12}^{loc} & = \int_{\xi=-1}^{1} \dfrac{1}{4} \left( 1 - \xi \right) \left( 1 + \xi \right) d \xi = \dfrac{1}{3} \end{aligned}\end{split}\]

#> Local mass and stiffness matrices
#> Length of the elements
el_len = np.array([ np.abs(rr[ee[iel,0]] - rr[ee[iel,1]]) for iel in np.arange(nelems) ])
el_jac = el_len * .5
# print(el_len)

#> Local mass matrix, stripped by row
me_loc = np.array([ 2., 1., 1., 2. ]) / 3.

#> Local stiffness matrix, stripped by row
ke_loc = np.array([ .5, -.5, -.5, .5 ])

#> Assembling global matrices
mi = np.concatenate(np.array([ [ e[0], e[0], e[1], e[1] ] for e in ee]))
mj = np.concatenate(np.array([ [ e[0], e[1], e[0], e[1] ] for e in ee]))

me = np.concatenate([ me_loc * m  * el_jac[iel] for iel in np.arange(nelems) ])
ke = np.concatenate([ ke_loc * EA / el_jac[iel] for iel in np.arange(nelems) ])

M = sp.sparse.coo_array((me, (mi, mj)),)
K = sp.sparse.coo_array((ke, (mi, mj)),)

# print(M.todense())
# print(K.todense())

9.1.1.1.1.1. Static problem

#> Prescribe boundary conditions
# 1) Slicing or 2) augmented system
#> 1) Slicing
# [ Kuu Kud ][ uu ] = [ fu ] + [ Fu ]
# [ Kdu Kdd ][ ud ]   [ fd ]   [ Fd ]
iD = i_dir.copy()
iU = list(set(np.arange(nnodes))-set(iD))

Kuu = (K.tocsc()[:,iU]).tocsr()[iU,:]
Kud = (K.tocsc()[:,iD]).tocsr()[iU,:]
Kdu = (K.tocsc()[:,iU]).tocsr()[iD,:]
Kdd = (K.tocsc()[:,iD]).tocsr()[iD,:]

F_vol = np.zeros(nnodes)
f = F_vol.copy()
f[i_neu] += F1
fu = f[iU]
fd = f[iD]

ud = u_dir.copy()

print(fu)
print(Kuu.todense())

#> Solve the linear problem, and retrieve the solution
uu = sp.sparse.linalg.spsolve(Kuu, fu - Kud @ ud )
fd = Kdu @ uu + Kdd @ ud

print(uu)

u = np.zeros(nnodes);  u[iU] = uu;  u[iD] = ud
f = np.zeros(nnodes);  f[iU] = fu;  f[iD] = fd

fig, ax = plt.subplots(1,1, figsize=(5,5))
ax.plot(rr, u)
ax.set_xlabel('x')
ax.set_ylabel('u')

[0. 0. 1.]
[[ 6. -3.  0.]
 [-3.  6. -3.]
 [ 0. -3.  3.]]
[0.33333333 0.66666667 1.        ]

Text(0, 0.5, 'u')

9.1.1.2. Torsional dynamics

\[\begin{aligned} I_z \ddot{\theta}_z - \left( GJ_z \theta'_z \right)' = m_z \end{aligned}\]

Same as axial dynamics…

9.1.1.3. Shear-Bending dynamics

\[\begin{split}\begin{cases} & m \ddot{s}_x - \left( \chi_x^{-1}GA ( s'_x - \theta_y ) \right)' = f_x \\ & I_y \ddot{\theta}_y - \left( EJ \theta_y' \right)' - \chi_x^{-1}GA \left( s_x' - \theta_y \right) = m_x \end{cases}\end{split}\]

\[\begin{split}\begin{cases} & m \ddot{s}_y - \dots = f_y \\ & I_x \ddot{\theta}_x - \dots = m_x \end{cases}\end{split}\]

Finite-dimensional approximation

\[\begin{split}\begin{aligned} s(x,t) & = \boldsymbol\phi_s^T(x) \mathbf{s}(t) \\ \theta(x,t) & = \boldsymbol\phi_\theta^T(x) \boldsymbol{\theta}(t) \\ \end{aligned}\end{split}\]

\[\begin{split}\begin{cases} \int_{\ell} \boldsymbol\phi_s m \boldsymbol\phi_s^T \, dx \, \ddot{\mathbf{s}} + \int_{\ell} \boldsymbol\phi'_s \chi^{-1} GA \left( \boldsymbol\phi_s^{' \ T} \mathbf{s} - \boldsymbol\phi^T_\theta \boldsymbol\theta \right) \, dx = \int_{\ell} \boldsymbol\phi^T_s \, f + \left.\left[ \boldsymbol\phi_s \underbrace{ \chi^{-1}GA( \boldsymbol\phi^{' \ T}_s \mathbf{s} - \boldsymbol\phi^T_\theta \boldsymbol\theta) }_{T} \right]\right|_{\partial \ell} \\ \int_{\ell} \boldsymbol\phi_\theta I \boldsymbol\phi_\theta^T \, dx \, \ddot{\boldsymbol{\theta}} + \int_{\ell} \boldsymbol\phi'_\theta EJ \boldsymbol\phi_\theta^{' \ T} \boldsymbol{\theta} - \int_{\ell} \boldsymbol\phi_\theta \chi^{-1} GA \left( \boldsymbol\phi^{' \ T}_s \mathbf{s} - \boldsymbol\phi^T_\theta \boldsymbol\theta \right)\, dx = m + \left.\left[ \boldsymbol\phi_\theta \underbrace{ EJ \boldsymbol\phi^{' \ T}_\theta \boldsymbol{\theta} }_{M} \right]\right|_{\partial \ell} \end{cases}\end{split}\]

\[\begin{split}\begin{bmatrix} \int_{\ell} \boldsymbol\phi_s m \boldsymbol\phi_s^T & 0 \\ 0 & \int_{\ell} \boldsymbol\phi_\theta I \boldsymbol\phi_\theta^T \end{bmatrix}\begin{bmatrix} \ddot{\mathbf{s}} \\ \ddot{\boldsymbol{\theta}} \end{bmatrix} + \begin{bmatrix} \int_{\ell} \boldsymbol\phi'_s \chi^{-1} GA \boldsymbol\phi_s^{' \ T} & -\int_{\ell} \boldsymbol\phi'_s \chi^{-1} GA \boldsymbol\phi_\theta^{T} \\ -\int_{\ell} \boldsymbol\phi_\theta \chi^{-1} GA \boldsymbol\phi_s^{' \ T} & \int_{\ell} \boldsymbol\phi'_\theta EJ \boldsymbol\phi_\theta^{' \ T} + \int_{\ell} \boldsymbol\phi_\theta \chi GA \boldsymbol\phi_\theta^{T} \end{bmatrix}\begin{bmatrix} \mathbf{s} \\ \boldsymbol{\theta} \end{bmatrix} = \begin{bmatrix} \dots \\ \dots \end{bmatrix}\end{split}\]

Different approximations for displacement and rotation

\[\begin{split}\begin{aligned} \phi^s_1(\xi) & = \dfrac{1}{2} \xi ( \xi - 1 ) && , \quad \phi^{s \ '}_1(\xi) = \xi - \dfrac{1}{2} \\ \phi^s_2(\xi) & = 1 - \xi^2 && , \quad \phi^{s \ '}_2(\xi) = - 2 \xi \\ \phi^s_3(\xi) & = \dfrac{1}{2} \xi ( \xi + 1 ) && , \quad \phi^{s \ '}_3(\xi) = \xi + \dfrac{1}{2} \\ \end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} \phi^\theta_1(\xi) & = \dfrac{1}{2} ( 1 - \xi ) && , \quad \phi^{\theta \ '}_1(\xi) = - \dfrac{1}{2} \\ \phi^\theta_2(\xi) & = \dfrac{1}{2} ( 1 + \xi ) && , \quad \phi^{\theta \ '}_2(\xi) = \dfrac{1}{2} \\ \end{aligned}\end{split}\]

\[\begin{split}\mathbf{M}_{ss} = \begin{bmatrix} \frac{4}{15} & \frac{2 }{15} & -\frac{1}{15} \\ \cdot & \frac{16}{15} & \frac{2}{15} \\ \cdot & \cdot & \frac{4}{15} \end{bmatrix}\end{split}\]

\[\begin{split}\mathbf{K}_{ss} = \begin{bmatrix} \frac{7}{6} & -\frac{4}{3} & \frac{1}{6} \\ \cdot & \frac{8}{3} & -\frac{4}{3} \\ \cdot & \cdot & \frac{7}{6} \end{bmatrix}\end{split}\]

\[\begin{split}\mathbf{A}_{s\theta} = \begin{bmatrix} -\frac{5}{6} & -\frac{1}{6} \\ \frac{2}{3} & -\frac{2}{3} \\ \frac{1}{6} & \frac{5}{6} \end{bmatrix}\end{split}\]

as

\[\begin{split}\begin{aligned} A_{11} & = \int_{-1}^{1} \phi^{s \ '}_1(\xi) \phi^{\theta}_1(\xi) \, d \xi = \int_{-1}^{1} \left( \xi - \frac{1}{2} \right) \frac{1}{2} \left( 1 - \xi \right) d \xi = \int_{-1}^{1} \left( -\frac{1}{4} - \frac{1}{2} \xi^2 \right) d \xi = - \dfrac{1}{2} - \dfrac{1}{3} = - \dfrac{5}{6} \\ A_{12} & = \int_{-1}^{1} \phi^{s \ '}_1(\xi) \phi^{\theta}_2(\xi) \, d \xi = \int_{-1}^{1} \left( \xi - \frac{1}{2} \right) \frac{1}{2} \left( 1 + \xi \right) d \xi = \int_{-1}^{1} \left( -\frac{1}{4} + \frac{1}{2} \xi^2 \right) d \xi = - \dfrac{1}{2} + \dfrac{1}{3} = - \dfrac{1}{6} \\ A_{21} & = \int_{-1}^{1} \phi^{s \ '}_2(\xi) \phi^{\theta}_1(\xi) \, d \xi = \int_{-1}^{1} \quad - 2 \xi \quad \ \frac{1}{2} \left( 1 - \xi \right) d \xi = \int_{-1}^{1} \qquad \quad \ \ \ \ \ \ \xi^2 \quad d \xi = \ \ \ \dfrac{2}{3} \\ A_{22} & = \int_{-1}^{1} \phi^{s \ '}_2(\xi) \phi^{\theta}_2(\xi) \, d \xi = \int_{-1}^{1} \quad - 2 \xi \quad \ \frac{1}{2} \left( 1 + \xi \right) d \xi = \int_{-1}^{1} \qquad \quad \ \ - \xi^2 \quad d \xi = - \dfrac{2}{3} \\ A_{31} & = \int_{-1}^{1} \phi^{s \ '}_3(\xi) \phi^{\theta}_1(\xi) \, d \xi = \int_{-1}^{1} \left( \xi + \frac{1}{2} \right) \frac{1}{2} \left( 1 - \xi \right) d \xi = \int_{-1}^{1} \left( \ \ \ \frac{1}{4} - \frac{1}{2} \xi^2 \right) d \xi = \ \ \ \dfrac{1}{2} - \dfrac{1}{3} = \dfrac{1}{6} \\ A_{32} & = \int_{-1}^{1} \phi^{s \ '}_3(\xi) \phi^{\theta}_2(\xi) \, d \xi = \int_{-1}^{1} \left( \xi + \frac{1}{2} \right) \frac{1}{2} \left( 1 + \xi \right) d \xi = \int_{-1}^{1} \left( \ \ \ \frac{1}{4} + \frac{1}{2} \xi^2 \right) d \xi = \ \ \ \dfrac{1}{2} + \dfrac{1}{3} = \dfrac{5}{6} \\ \end{aligned}\end{split}\]

%reset -f

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

#> Parameters
#> Physical parameters
l = 1.                  # length
m, I = 1., 1.           # mass and inertia density
EJ = 1.                 # bending stiffness
xGA = 10000. * EJ        # shear stiffness (relative to bending stiffness)
qf, qm = .0, .0         # distributed loads, force and moments

#> Discretization parameters, for building mesh
nelems = 10
nnodes_p1 = nelems + 1
nnodes_p2 = nelems
nnodes = nnodes_p1 + nnodes_p2

#> Boundary conditions
u0, theta0 = .0, .0    # clamp in P0
F1, M1 = 1., .0        # only force in P1

#> Mesh
# Not only P1 elements, but P1(rotation)-P2(displacement) elements
#> Node coordinates (nnodes,nd), with space dimension nd=1
rr = np.linspace(0, 1, nnodes_p1) * l                 # P1 nodes (equi-spaced)
rr = np.concatenate([rr, .5*(rr[1:] + rr[:-1])])      # concatentate P1 with P2 nodes
#> Node-element connectivity matrix (nelems,3)
# First two ids: P1-nodes, third id: P2-node
ee = np.array([ [ iel, iel+1, nnodes_p1+iel ] for iel in np.arange(nelems) ])

#> Global indices
# Some nodes have 2 variables (P1 have both displacement and rotation). In order to 
# assemble global matrices, global indices of the variables are required
i_displ = ee.copy()
i_theta = ee[:,:2].copy() + nnodes
i_var = np.block([i_displ, i_theta])

n_vars = nnodes + nnodes_p1

print(f"rr:\n {rr}")
print(f"ee:\n {ee}")
print(f"i_var: \n {i_var}")

rr:
 [0.   0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.   0.05 0.15 0.25
 0.35 0.45 0.55 0.65 0.75 0.85 0.95]
ee:
 [[ 0  1 11]
 [ 1  2 12]
 [ 2  3 13]
 [ 3  4 14]
 [ 4  5 15]
 [ 5  6 16]
 [ 6  7 17]
 [ 7  8 18]
 [ 8  9 19]
 [ 9 10 20]]
i_var: 
 [[ 0  1 11 21 22]
 [ 1  2 12 22 23]
 [ 2  3 13 23 24]
 [ 3  4 14 24 25]
 [ 4  5 15 25 26]
 [ 5  6 16 26 27]
 [ 6  7 17 27 28]
 [ 7  8 18 28 29]
 [ 8  9 19 29 30]
 [ 9 10 20 30 31]]

#> Local mass and stiffness matrices
#> Length of the elements
el_len = np.array([ np.abs(rr[ee[iel,0]] - rr[ee[iel,1]]) for iel in np.arange(nelems) ])
el_jac = el_len * .5
# print(el_len)

#> Local mass matrix, stripped by row
me_tt_loc = np.array([ 2., 1., 1., 2. ]) / 3.
me_ss_loc = np.array([ 4., -1., 2., -1., 4., 2., 2., 2., 16.]) / 15.

#> Local stiffness matrix, stripped by row
ke_tt_loc = np.array([ .5, -.5, -.5, .5 ])
ke_ss_loc = np.array([ 7., 1., -8., 1., 7., -8., -8., -8., 16.]) / 6.

#> Ast
# [ -5, -1 ]
# [  1,  5 ]
# [  4, -4 ] / 6
ae_st_loc = np.array([-5., -1., 1., 5., 4., -4.]) / 6.
ae_ts_loc = np.array([-5., 1., 4., -1., 5., -4.]) / 6.

#> Assembling global matrices
#> Mass
mi_ss = np.concatenate(np.array([ np.repeat(i_var[iel,:3], [3,3,3]) for iel in np.arange(nelems)]))
mj_ss = np.concatenate(np.array([ np.tile(i_var[iel,:3],3) for iel in np.arange(nelems)]))
me_ss = np.concatenate([ me_ss_loc * m * el_jac[iel] for iel in np.arange(nelems)])

mi_tt = np.concatenate(np.array([ np.repeat(i_var[iel,3:], [2,2]) for iel in np.arange(nelems)]))
mj_tt = np.concatenate(np.array([ np.tile(i_var[iel,3:],2) for iel in np.arange(nelems)]))
me_tt = np.concatenate([me_tt_loc * I * el_jac[iel] for iel in np.arange(nelems)])

mi = np.concatenate([mi_ss, mi_tt])
mj = np.concatenate([mj_ss, mj_tt])
me = np.concatenate([me_ss, me_tt])

M = sp.sparse.coo_array((me, (mi,mj)),)

#> Stiffness
ke_ss    = np.concatenate([ ke_ss_loc * xGA / el_jac[iel] for iel in np.arange(nelems)])
ke_tt_EJ = np.concatenate([ ke_tt_loc * EJ  / el_jac[iel] for iel in np.arange(nelems)])
ke_tt_GA = np.concatenate([ me_tt_loc * xGA * el_jac[iel] for iel in np.arange(nelems)])
ke_tt    = ke_tt_EJ + ke_tt_GA

ai_st = np.concatenate(np.array([ np.repeat(i_var[iel,:3], [2,2,2]) for iel in np.arange(nelems)]))
aj_st = np.concatenate(np.array([ np.tile(i_var[iel,3:],3) for iel in np.arange(nelems)]))

ai_ts = np.concatenate(np.array([ np.repeat(i_var[iel,3:], [3,3]) for iel in np.arange(nelems)]))
aj_ts = np.concatenate(np.array([ np.tile(i_var[iel,:3],2) for iel in np.arange(nelems)]))

ke_st    = np.concatenate([ - ae_st_loc * xGA for iel in np.arange(nelems)])
ke_ts    = np.concatenate([ - ae_ts_loc * xGA for iel in np.arange(nelems)])

ki = np.concatenate([ mi_ss, ai_st, ai_ts, mi_tt ])
kj = np.concatenate([ mj_ss, aj_st, aj_ts, mj_tt ])
ke = np.concatenate([ ke_ss, ke_st, ke_ts, ke_tt ])

K = sp.sparse.coo_array((ke, (ki,kj)),)

# Some check/output for low-dimensional problems
if ( nelems < 4 ):
  print("\n> Mass matrix:")
  print(M.todense())
  print(f"Sum of elements of mass matrix: {np.sum(M)}")
  print("Needs to be equal to the total mass of the system")

  print("\n> Stiffness matrix:")
  print(K.todense())

  if ( nelems < 3 ):
    print("Null space:")
    print("Comment: 2-dimensional null space (rigid translation and rotation, or their linear combination)")
    print(sp.linalg.null_space(K.todense()))

9.1.1.3.1. Example: clamped beam

# HERE, a very rude implementation:
# a good design needs a cleaner impelemntation of boundary conditions and loads
# and can be easily (with some work) achieved relying on OO-programming, relying
# on classes

# TODO: could it be a good idea to provide an example of a FEM code implementation,
# - clearly separating pre-processing, analyses, and post-processing
# - relying on OO-programming
# - factor method for maintenability and expansion with the min amount of pain
# - ...
# ???

#> Essential boundary conditions
iD = np.array([ 0, nnodes ])   # var index of displ and rot of the first node
uD = np.array([ .0, .0 ])      # value of the displ and rot of the first node

iN = np.array([ nnodes_p1 - 1 ])    # var index of the displ of the last node
fN = np.array([ 1. ])

iU = list(set(np.arange(n_vars)) - set(iD))

#> Slicing
Kuu = (K.tocsc()[:,iU]).tocsr()[iU,:]
Kud = (K.tocsc()[:,iD]).tocsr()[iU,:]
Kdu = (K.tocsc()[:,iU]).tocsr()[iD,:]
Kdd = (K.tocsc()[:,iD]).tocsr()[iD,:]

F_vol = np.zeros(n_vars)
f = F_vol.copy()
f[iN] += fN
fu = f[iU]
fd = f[iD]

#> Solve the linear problem, and retrieve the solution
uu = sp.sparse.linalg.spsolve(Kuu, fu - Kud @ uD )
fd = Kdu @ uu + Kdd @ uD

nU = len(iU)
u = np.zeros(n_vars);  u[iU] = uu;  u[iD] = uD

u_displ = u[:nnodes]
u_rot   = u[nnodes:]

fig, ax = plt.subplots(2,1, figsize=(5,5))
ax[0].plot(rr            , u_displ, 'o')
ax[1].plot(rr[:nnodes_p1], u_rot  , 'o')
ax[0].set_ylabel('Displ.')
ax[1].set_ylabel('Rot.')
ax[1].set_xlabel('x')
ax[0].grid()
ax[1].grid()

print(f"u_displ: {u_displ}")
print(f"u_rot  : {u_rot  }")

u_displ: [0.        0.00476   0.01852   0.04028   0.06904   0.1038    0.14356
 0.18732   0.23408   0.28284   0.3326    0.0011925 0.0105775 0.0284625
 0.0538475 0.0857325 0.1231175 0.1650025 0.2103875 0.2582725 0.3076575]
u_rot  : [0.    0.095 0.18  0.255 0.32  0.375 0.42  0.455 0.48  0.495 0.5  ]