Local stiffness matrix Euler-Bernoulli element

Local stiffness matrix Euler-Bernoulli element#

In Force-displacement relations single extension element and Combine elements using matrix formulation you’ve seen how to derive the local stiffness matrix for a simple extension element. But can you do the same for other elements?

Previously, the local stiffness matrix for a simple extension element was found as:

\[\begin{split}\mathbf{K}^{(e)} = \cfrac{EA}{\ell}\begin{bmatrix} 1&-1\\-1&1 \end{bmatrix}\end{split}\]

The same procedure can be followed for other element, like a combined extension and Euler-Bernoulli element:

../_images/elemtypes.svg

Fig. 12 Combined extension and Euler-Bernoulli element#

The amount of degrees of freedom increases, as both ends of the element can translate in two directions and rotate. However, the approach is exactly the same, leading to the following element stiffness matrix:

\[\begin{split} \mathbf{K}^{(e)} = \begin{bmatrix} \cfrac{EA}{\ell} & 0 & 0 & -\cfrac{EA}{\ell} & 0 & 0\\ 0 & \cfrac{12EI}{\ell^3} & -\cfrac{6EI}{\ell^2} & 0 & -\cfrac{12EI}{\ell^3} & -\cfrac{6EI}{\ell^2}\\ 0 & -\cfrac{6EI}{\ell^2} & \cfrac{4EI}{\ell} & 0 & \cfrac{6EI}{\ell^2} & \cfrac{2EI}{\ell}\\ -\cfrac{EA}{\ell} & 0 & 0 & \cfrac{EA}{\ell} & 0 & 0\\ 0 & -\cfrac{12EI}{\ell^3} & \cfrac{6EI}{\ell^2} & 0 & \cfrac{12EI}{\ell^3} & \cfrac{6EI}{\ell^2}\\ 0 & -\cfrac{6EI}{\ell^2} & \cfrac{2EI}{\ell} & 0 & \cfrac{6EI}{\ell^2} & \cfrac{4EI}{\ell}\\ \end{bmatrix}\end{split}\]

for

\[\begin{split} \mathbf{u}^{(e)} =\begin{bmatrix} u_1\\ w_1 \\ \varphi_1 \\ u_2 \\ w_2\\ \varphi_2\\ \end{bmatrix} \end{split}\]

Derivation using SymPy#

We can make use of software like SymPy, as we did before in Solving differential equation of one field using SymPy to do the calculations in this derivation:

import sympy as sym
sym.init_printing()

import sympy as sym

EI, x, L = sym.symbols('EI, x, L')
w = sym.Function('w')

ODE_bending = sym.Eq(w(x).diff(x, 4) * EI, 0)
display(ODE_bending)
\[\displaystyle EI \frac{d^{4}}{d x^{4}} w{\left(x \right)} = 0\]
w = sym.dsolve(ODE_bending, w(x)).rhs
display(w)
\[\displaystyle C_{1} + C_{2} x + C_{3} x^{2} + C_{4} x^{3}\]
phi = -w.diff(x)
kappa = phi.diff(x)
M = EI * kappa
V = M.diff(x)

w_1, w_2, phi_1, phi_2 = sym.symbols('w_1, w_2, phi_1, phi_2')

eq1 = sym.Eq(w.subs(x,0),w_1)
eq2 = sym.Eq(w.subs(x,L),w_2)
eq3 = sym.Eq(phi.subs(x,0),phi_1)
eq4 = sym.Eq(phi.subs(x,L),phi_2)

sol = sym.solve([eq1, eq2, eq3, eq4 ], sym.symbols('C1, C2, C3, C4'))
for key in sol:
    display(sym.Eq(key, sol[key]))
\[\displaystyle C_{1} = w_{1}\]
\[\displaystyle C_{2} = - \phi_{1}\]
\[\displaystyle C_{3} = \frac{2 L \phi_{1} + L \phi_{2} - 3 w_{1} + 3 w_{2}}{L^{2}}\]
\[\displaystyle C_{4} = \frac{- L \phi_{1} - L \phi_{2} + 2 w_{1} - 2 w_{2}}{L^{3}}\]
F_1_z, F_2_z, T_1_y, T_2_y = sym.symbols('F_1_z, F_2_z, T_1_y, T_2_y')

eq5 = sym.Eq(-V.subs(sol).subs(x,0), F_1_z)
eq6 = sym.Eq(V.subs(sol).subs(x,L), F_2_z)
eq7 = sym.Eq(-M.subs(sol).subs(x,0), T_1_y)
eq8 = sym.Eq(M.subs(sol).subs(x,L), T_2_y)

K_e, f_e = sym.linear_eq_to_matrix([eq5,eq7, eq6, eq8], [w_1, phi_1, w_2, phi_2])
display(K_e)
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{12 EI}{L^{3}} & - \frac{6 EI}{L^{2}} & - \frac{12 EI}{L^{3}} & - \frac{6 EI}{L^{2}}\\- \frac{6 EI}{L^{2}} & \frac{4 EI}{L} & \frac{6 EI}{L^{2}} & \frac{2 EI}{L}\\- \frac{12 EI}{L^{3}} & \frac{6 EI}{L^{2}} & \frac{12 EI}{L^{3}} & \frac{6 EI}{L^{2}}\\- \frac{6 EI}{L^{2}} & \frac{2 EI}{L} & \frac{6 EI}{L^{2}} & \frac{4 EI}{L}\end{matrix}\right]\end{split}\]

To use the stiffness matrix without manually copying it over, you can make use of the lambdify which converts a symbolic SymPy object in a python function. This allows you to evaluate it for specific numerical values and continue using it in the numerical framework of the matrix method.

K = sym.lambdify((L, EI), K_e)
print(K.__doc__)

Created with lambdify. Signature:

func(L, EI)

Expression:

Matrix([[12*EI/L**3, -6*EI/L**2, -12*EI/L**3, -6*EI/L**2], [-6*EI/L**2,...

Source code:

def _lambdifygenerated(L, EI):
    return array([[12*EI/L**3, -6*EI/L**2, -12*EI/L**3, -6*EI/L**2], [-6*EI/L**2, 4*EI/L, 6*EI/L**2, 2*EI/L], [-12*EI/L**3, 6*EI/L**2, 12*EI/L**3, 6*EI/L**2], [-6*EI/L**2, 2*EI/L, 6*EI/L**2, 4*EI/L]])


Imported modules:

print('Example of K with L=5 and EI=1000:\n',K(5,1000))

Example of K with L=5 and EI=1000:
 [[  96. -240.  -96. -240.]
 [-240.  800.  240.  400.]
 [ -96.  240.   96.  240.]
 [-240.  400.  240.  800.]]