# Beam finite element matrices
Based on: \
[1] Neto, M. A., Amaro, A., Roseiro, L., Cirne, J., & Leal, R. (2015). Finite Element Method for Beams. In Engineering Computation of Structures: The Finite Element Method (pp. 115–156). Springer International Publishing. http://doi.org/10.1007/978-3-319-17710-6_4


The beam finite element, its nodes and degrees-of-freedom (DOFs) can be seen below:



```python
from sympy import *
init_printing()


def symb(x, y):
    return symbols('{0}_{1}'.format(x, y), type = float)


E, A, L, G, I_1, I_2, I_3, rho = symbols('E A L G I_1 I_2 I_3 rho', type = float)
```

The finite element matrices should have order 12, accounting for each of the element's DOFs, shown below:


```python
u = Matrix(12, 1, [symb('u', v + 1) for v in range(12)])
transpose(u)
```




$\displaystyle \left[\begin{array}{cccccccccccc}u_{1} & u_{2} & u_{3} & u_{4} & u_{5} & u_{6} & u_{7} & u_{8} & u_{9} & u_{10} & u_{11} & u_{12}\end{array}\right]$



## Axial deformation along $x_1$
In terms of generic coordinates $v_i$:


```python
v_a = Matrix(2, 1, [symb('v', v + 1) for v in range(2)])
transpose(v_a)
```




$\displaystyle \left[\begin{matrix}v_{1} & v_{2}\end{matrix}\right]$



which are equivalent to the $u_i$ coordinates in the following way:
$$
    \mathbf{v} = \mathbf{R} \mathbf{u},
$$
where:
$$
    v_1 = u_1, \\
    v_2 = u_7,
$$
with the following coordinate transformation matrix:


```python
R = zeros(12)
R[1 - 1, 1 - 1] = 1
R[2 - 1, 7 - 1] = 1
R[:2, :]
```




$\displaystyle \left[\begin{array}{cccccccccccc}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\end{array}\right]$



### Stiffness matrix
Eq. (3.15) of [1]:


```python
K_a = (E * A / L) * Matrix([[1, -1],
                            [-1, 1]])
K_a
```




$\displaystyle \left[\begin{matrix}\frac{A E}{L} & - \frac{A E}{L}\\- \frac{A E}{L} & \frac{A E}{L}\end{matrix}\right]$



### Inertia matrix
Eq. (3.16) of [1]:


```python
M_a = (rho * A * L / 6) * Matrix([[2, 1],
                                  [1, 2]])
M_a
```




$\displaystyle \left[\begin{matrix}\frac{A L \rho}{3} & \frac{A L \rho}{6}\\\frac{A L \rho}{6} & \frac{A L \rho}{3}\end{matrix}\right]$



## Torsional deformation around $x_1$
According to [1], one can obtain the matrices for the torsional case from the axial case by replacing the elasticity modulus $E$ and the cross-sectional area $A$ by the shear modulus $G$ and the polar area moment of inertia $I_1$.

In terms of generic coordinates $v_i$:


```python
v_t = Matrix(2, 1, [symb('v', v + 3) for v in range(2)])
transpose(v_t)
```




$\displaystyle \left[\begin{matrix}v_{3} & v_{4}\end{matrix}\right]$



which are equivalent to the $u_i$ coordinates in the following way:
$$
    v_3 = u_4, \\
    v_4 = u_{10},
$$
with the following coordinate transformation matrix:


```python
R[3 - 1,  4 - 1] = 1
R[4 - 1, 10 - 1] = 1
R[0:4, :]
```




$\displaystyle \left[\begin{array}{cccccccccccc}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\end{array}\right]$



### Stiffness matrix


```python
K_t = K_a.subs([(E, G), (A, I_1)])
K_t
```




$\displaystyle \left[\begin{matrix}\frac{G I_{1}}{L} & - \frac{G I_{1}}{L}\\- \frac{G I_{1}}{L} & \frac{G I_{1}}{L}\end{matrix}\right]$



### Inertia matrix


```python
M_t = M_a.subs([(E, G), (A, I_1)])
M_t
```




$\displaystyle \left[\begin{matrix}\frac{I_{1} L \rho}{3} & \frac{I_{1} L \rho}{6}\\\frac{I_{1} L \rho}{6} & \frac{I_{1} L \rho}{3}\end{matrix}\right]$



## Bending on the plane $x_1-x_3$
In this case the bending torsion occurs around the $x_2$ axis. In terms of generic coordinates $v_i$:


```python
v_b13 = Matrix(4, 1, [symb('v', v + 9) for v in range(4)])
transpose(v_b13)
```




$\displaystyle \left[\begin{matrix}v_{9} & v_{10} & v_{11} & v_{12}\end{matrix}\right]$



which are equivalent to the $u_i$ coordinates in the following way:
$$
    v_9    = u_3, \\
    v_{10} = u_5, \\
    v_{11} = u_9, \\
    v_{12} = u_{11},
$$
with the following coordinate transformation matrix:


```python
R[9  - 1,  3 - 1] = 1
R[10 - 1,  5 - 1] = 1
R[11 - 1,  9 - 1] = 1
R[12 - 1, 11 - 1] = 1
R
```




$\displaystyle \left[\begin{array}{cccccccccccc}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0\end{array}\right]$



### Stiffness matrix


```python
K_b13 = (E * I_2 / 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]])

if (not K_b13.is_symmetric()):
    print('Error in K_b13.')
K_b13
```




$\displaystyle \left[\begin{matrix}\frac{12 E I_{2}}{L^{3}} & \frac{6 E I_{2}}{L^{2}} & - \frac{12 E I_{2}}{L^{3}} & \frac{6 E I_{2}}{L^{2}}\\\frac{6 E I_{2}}{L^{2}} & \frac{4 E I_{2}}{L} & - \frac{6 E I_{2}}{L^{2}} & \frac{2 E I_{2}}{L}\\- \frac{12 E I_{2}}{L^{3}} & - \frac{6 E I_{2}}{L^{2}} & \frac{12 E I_{2}}{L^{3}} & - \frac{6 E I_{2}}{L^{2}}\\\frac{6 E I_{2}}{L^{2}} & \frac{2 E I_{2}}{L} & - \frac{6 E I_{2}}{L^{2}} & \frac{4 E I_{2}}{L}\end{matrix}\right]$



### Inertia matrix


```python
M_b13 = (rho * L * A / 420) * Matrix([[ 156    ,  22 * L   ,   54    , -13 * L   ],
                                      [  22 * L,   4 * L**2,   13 * L, - 3 * L**2],
                                      [  54    ,  13 * L   ,  156    , -22 * L   ],
                                      [- 13 * L, - 3 * L**2, - 22 * L,   4 * L**2]])

if (not M_b13.is_symmetric()):
    print('Error in M_b13.')
M_b13
```




$\displaystyle \left[\begin{matrix}\frac{13 A L \rho}{35} & \frac{11 A L^{2} \rho}{210} & \frac{9 A L \rho}{70} & - \frac{13 A L^{2} \rho}{420}\\\frac{11 A L^{2} \rho}{210} & \frac{A L^{3} \rho}{105} & \frac{13 A L^{2} \rho}{420} & - \frac{A L^{3} \rho}{140}\\\frac{9 A L \rho}{70} & \frac{13 A L^{2} \rho}{420} & \frac{13 A L \rho}{35} & - \frac{11 A L^{2} \rho}{210}\\- \frac{13 A L^{2} \rho}{420} & - \frac{A L^{3} \rho}{140} & - \frac{11 A L^{2} \rho}{210} & \frac{A L^{3} \rho}{105}\end{matrix}\right]$



## Bending on the plane $x_1-x_2$
In this case the bending torsion occurs around the $x_3$ axis, but in the opposite direction. The matrices are similar to the case of bending on the $x_1-x_3$ plane, needing proper coordinate transformation and replacing the index of the area moment of inertia from 2 to 3.

Written in terms of generic coordinates $v_i$:


```python
v_b12 = Matrix(4, 1, [symb('v', v + 5) for v in range(4)])
transpose(v_b12)
```




$\displaystyle \left[\begin{matrix}v_{5} & v_{6} & v_{7} & v_{8}\end{matrix}\right]$



which are equivalent to the $u_i$ coordinates in the following way:
$$
    v_5 =  u_2, \\
    v_6 = -u_6, \\
    v_7 =  u_8, \\
    v_8 = -u_{12},
$$
with the following coordinate transformation matrix:


```python
R[5 - 1,  2 - 1] =  1
R[6 - 1,  6 - 1] = -1
R[7 - 1,  8 - 1] =  1
R[8 - 1, 12 - 1] = -1
R
```

### Stiffness matrix


```python
K_b12 = K_b13.subs(I_2, I_3)
K_b12
```




$\displaystyle \left[\begin{matrix}\frac{12 E I_{3}}{L^{3}} & \frac{6 E I_{3}}{L^{2}} & - \frac{12 E I_{3}}{L^{3}} & \frac{6 E I_{3}}{L^{2}}\\\frac{6 E I_{3}}{L^{2}} & \frac{4 E I_{3}}{L} & - \frac{6 E I_{3}}{L^{2}} & \frac{2 E I_{3}}{L}\\- \frac{12 E I_{3}}{L^{3}} & - \frac{6 E I_{3}}{L^{2}} & \frac{12 E I_{3}}{L^{3}} & - \frac{6 E I_{3}}{L^{2}}\\\frac{6 E I_{3}}{L^{2}} & \frac{2 E I_{3}}{L} & - \frac{6 E I_{3}}{L^{2}} & \frac{4 E I_{3}}{L}\end{matrix}\right]$



### Inertia matrix


```python
M_b12 = M_b13
if (not M_b12.is_symmetric()):
    print('Error in M_b12.')
M_b12
```




$\displaystyle \left[\begin{matrix}\frac{13 A L \rho}{35} & \frac{11 A L^{2} \rho}{210} & \frac{9 A L \rho}{70} & - \frac{13 A L^{2} \rho}{420}\\\frac{11 A L^{2} \rho}{210} & \frac{A L^{3} \rho}{105} & \frac{13 A L^{2} \rho}{420} & - \frac{A L^{3} \rho}{140}\\\frac{9 A L \rho}{70} & \frac{13 A L^{2} \rho}{420} & \frac{13 A L \rho}{35} & - \frac{11 A L^{2} \rho}{210}\\- \frac{13 A L^{2} \rho}{420} & - \frac{A L^{3} \rho}{140} & - \frac{11 A L^{2} \rho}{210} & \frac{A L^{3} \rho}{105}\end{matrix}\right]$



## Assembly of the full matrices
Accounting for axial loads, torques and bending in both planes.


```python
RAR = lambda A: transpose(R)*A*R
transpose(R**-1*u)
```




$\displaystyle \left[\begin{array}{cccccccccccc}u_{1} & u_{5} & u_{9} & u_{3} & u_{10} & - u_{6} & u_{2} & u_{7} & u_{11} & u_{4} & u_{12} & - u_{8}\end{array}\right]$




```python
K_f = diag(K_a, K_t, K_b12, K_b13)
K = RAR(K_f)
if (not K.is_symmetric()):
    print('Error in K.')
K
```




$\displaystyle \left[\begin{array}{cccccccccccc}\frac{A E}{L} & 0 & 0 & 0 & 0 & 0 & - \frac{A E}{L} & 0 & 0 & 0 & 0 & 0\\0 & \frac{12 E I_{3}}{L^{3}} & 0 & 0 & 0 & - \frac{6 E I_{3}}{L^{2}} & 0 & - \frac{12 E I_{3}}{L^{3}} & 0 & 0 & 0 & - \frac{6 E I_{3}}{L^{2}}\\0 & 0 & \frac{12 E I_{2}}{L^{3}} & 0 & \frac{6 E I_{2}}{L^{2}} & 0 & 0 & 0 & - \frac{12 E I_{2}}{L^{3}} & 0 & \frac{6 E I_{2}}{L^{2}} & 0\\0 & 0 & 0 & \frac{G I_{1}}{L} & 0 & 0 & 0 & 0 & 0 & - \frac{G I_{1}}{L} & 0 & 0\\0 & 0 & \frac{6 E I_{2}}{L^{2}} & 0 & \frac{4 E I_{2}}{L} & 0 & 0 & 0 & - \frac{6 E I_{2}}{L^{2}} & 0 & \frac{2 E I_{2}}{L} & 0\\0 & - \frac{6 E I_{3}}{L^{2}} & 0 & 0 & 0 & \frac{4 E I_{3}}{L} & 0 & \frac{6 E I_{3}}{L^{2}} & 0 & 0 & 0 & \frac{2 E I_{3}}{L}\\- \frac{A E}{L} & 0 & 0 & 0 & 0 & 0 & \frac{A E}{L} & 0 & 0 & 0 & 0 & 0\\0 & - \frac{12 E I_{3}}{L^{3}} & 0 & 0 & 0 & \frac{6 E I_{3}}{L^{2}} & 0 & \frac{12 E I_{3}}{L^{3}} & 0 & 0 & 0 & \frac{6 E I_{3}}{L^{2}}\\0 & 0 & - \frac{12 E I_{2}}{L^{3}} & 0 & - \frac{6 E I_{2}}{L^{2}} & 0 & 0 & 0 & \frac{12 E I_{2}}{L^{3}} & 0 & - \frac{6 E I_{2}}{L^{2}} & 0\\0 & 0 & 0 & - \frac{G I_{1}}{L} & 0 & 0 & 0 & 0 & 0 & \frac{G I_{1}}{L} & 0 & 0\\0 & 0 & \frac{6 E I_{2}}{L^{2}} & 0 & \frac{2 E I_{2}}{L} & 0 & 0 & 0 & - \frac{6 E I_{2}}{L^{2}} & 0 & \frac{4 E I_{2}}{L} & 0\\0 & - \frac{6 E I_{3}}{L^{2}} & 0 & 0 & 0 & \frac{2 E I_{3}}{L} & 0 & \frac{6 E I_{3}}{L^{2}} & 0 & 0 & 0 & \frac{4 E I_{3}}{L}\end{array}\right]$




```python
M_f = diag(M_a, M_t, M_b12, M_b13)
M = RAR(M_f)
if (not M.is_symmetric()):
    print('Error in M.')
M
```




$\displaystyle \left[\begin{array}{cccccccccccc}\frac{A L \rho}{3} & 0 & 0 & 0 & 0 & 0 & \frac{A L \rho}{6} & 0 & 0 & 0 & 0 & 0\\0 & \frac{13 A L \rho}{35} & 0 & 0 & 0 & - \frac{11 A L^{2} \rho}{210} & 0 & \frac{9 A L \rho}{70} & 0 & 0 & 0 & \frac{13 A L^{2} \rho}{420}\\0 & 0 & \frac{13 A L \rho}{35} & 0 & \frac{11 A L^{2} \rho}{210} & 0 & 0 & 0 & \frac{9 A L \rho}{70} & 0 & - \frac{13 A L^{2} \rho}{420} & 0\\0 & 0 & 0 & \frac{I_{1} L \rho}{3} & 0 & 0 & 0 & 0 & 0 & \frac{I_{1} L \rho}{6} & 0 & 0\\0 & 0 & \frac{11 A L^{2} \rho}{210} & 0 & \frac{A L^{3} \rho}{105} & 0 & 0 & 0 & \frac{13 A L^{2} \rho}{420} & 0 & - \frac{A L^{3} \rho}{140} & 0\\0 & - \frac{11 A L^{2} \rho}{210} & 0 & 0 & 0 & \frac{A L^{3} \rho}{105} & 0 & - \frac{13 A L^{2} \rho}{420} & 0 & 0 & 0 & - \frac{A L^{3} \rho}{140}\\\frac{A L \rho}{6} & 0 & 0 & 0 & 0 & 0 & \frac{A L \rho}{3} & 0 & 0 & 0 & 0 & 0\\0 & \frac{9 A L \rho}{70} & 0 & 0 & 0 & - \frac{13 A L^{2} \rho}{420} & 0 & \frac{13 A L \rho}{35} & 0 & 0 & 0 & \frac{11 A L^{2} \rho}{210}\\0 & 0 & \frac{9 A L \rho}{70} & 0 & \frac{13 A L^{2} \rho}{420} & 0 & 0 & 0 & \frac{13 A L \rho}{35} & 0 & - \frac{11 A L^{2} \rho}{210} & 0\\0 & 0 & 0 & \frac{I_{1} L \rho}{6} & 0 & 0 & 0 & 0 & 0 & \frac{I_{1} L \rho}{3} & 0 & 0\\0 & 0 & - \frac{13 A L^{2} \rho}{420} & 0 & - \frac{A L^{3} \rho}{140} & 0 & 0 & 0 & - \frac{11 A L^{2} \rho}{210} & 0 & \frac{A L^{3} \rho}{105} & 0\\0 & \frac{13 A L^{2} \rho}{420} & 0 & 0 & 0 & - \frac{A L^{3} \rho}{140} & 0 & \frac{11 A L^{2} \rho}{210} & 0 & 0 & 0 & \frac{A L^{3} \rho}{105}\end{array}\right]$



d1bf5ee7f5cec8df1dd856ade52e8207d8b3e340
ipynb
Jupyter Notebook
notes/.ipynb_checkpoints/FEM_beam-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
2020-10-17T13:43:01.000Z
2020-10-17T13:43:01.000Z
notes/.ipynb_checkpoints/FEM_beam-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
notes/.ipynb_checkpoints/FEM_beam-checkpoint.ipynb
gfsReboucas/Drivetrain-python
90cc8a0b26fa6dd851a8ddaaf321f5ae9f5cf431
MIT
Qwen/Qwen-72B
1. YES
2. YES