# Eigenvalues and eigenvectors of stiffness matrices



In [1]:

from sympy.utilities.codegen import codegen
from sympy import *
from sympy import init_printing




In [2]:

init_printing()




In [3]:

r, s, t, x, y, z = symbols('r s t x y z')
k, m, n = symbols('k m n', integer=True)
rho, nu, E = symbols('rho, nu, E')



## Predefinition

The constitutive model tensor in Voigt notation (plane stress) is $$C = \frac{E}{(1 - \nu^2)} \begin{pmatrix} 1 & \nu & 0\\ \nu & 1 & 0\\ 0 & 0 & \frac{1 - \nu}{2)} \end{pmatrix}$$



In [4]:

K_factor = E/(1 - nu**2)
C = K_factor * Matrix([
[1, nu, 0],
[nu, 1, 0],
[0, 0, (1 - nu)/2]])

C




Out[4]:

$\displaystyle \left[\begin{matrix}\frac{E}{1 - \nu^{2}} & \frac{E \nu}{1 - \nu^{2}} & 0\\\frac{E \nu}{1 - \nu^{2}} & \frac{E}{1 - \nu^{2}} & 0\\0 & 0 & \frac{E \left(\frac{1}{2} - \frac{\nu}{2}\right)}{1 - \nu^{2}}\end{matrix}\right]$



## Interpolation functions

The shape functions are



In [5]:

N = S(1)/4*Matrix([(1 + r)*(1 + s),
(1 - r)*(1 + s),
(1 - r)*(1 - s),
(1 + r)*(1 - s)])
N




Out[5]:

$\displaystyle \left[\begin{matrix}\frac{\left(r + 1\right) \left(s + 1\right)}{4}\\\frac{\left(1 - r\right) \left(s + 1\right)}{4}\\\frac{\left(1 - r\right) \left(1 - s\right)}{4}\\\frac{\left(1 - s\right) \left(r + 1\right)}{4}\end{matrix}\right]$



Thus, the interpolation matrix renders

## Derivatives interpolation matrix



In [6]:

dHdr = zeros(2,4)
for i in range(4):
dHdr[0,i] = diff(N[i],r)
dHdr[1,i] = diff(N[i],s)

jaco = eye(2)  # Jacobian matrix, identity for now
dHdx = jaco*dHdr

B = zeros(3,8)
for i in range(4):
B[0, 2*i] = dHdx[0, i]
B[1, 2*i+1] = dHdx[1, i]
B[2, 2*i] = dHdx[1, i]
B[2, 2*i+1] = dHdx[0, i]

B




Out[6]:

$\displaystyle \left[\begin{matrix}\frac{s}{4} + \frac{1}{4} & 0 & - \frac{s}{4} - \frac{1}{4} & 0 & \frac{s}{4} - \frac{1}{4} & 0 & \frac{1}{4} - \frac{s}{4} & 0\\0 & \frac{r}{4} + \frac{1}{4} & 0 & \frac{1}{4} - \frac{r}{4} & 0 & \frac{r}{4} - \frac{1}{4} & 0 & - \frac{r}{4} - \frac{1}{4}\\\frac{r}{4} + \frac{1}{4} & \frac{s}{4} + \frac{1}{4} & \frac{1}{4} - \frac{r}{4} & - \frac{s}{4} - \frac{1}{4} & \frac{r}{4} - \frac{1}{4} & \frac{s}{4} - \frac{1}{4} & - \frac{r}{4} - \frac{1}{4} & \frac{1}{4} - \frac{s}{4}\end{matrix}\right]$



Being the stiffness matrix integrand $$K_\text{int} = B^T C B$$



In [7]:

K_int = B.T*C*B



## Analytic integration

The stiffness matrix is obtained integrating the product of the interpolator-derivatives (displacement-to-strains) matrix with the constitutive tensor and itself, i.e. \begin{align*} K &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} K_\text{int} dr\, ds\\ &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} B^T C\, B\, dr\, ds \enspace . \end{align*}



In [8]:

K = zeros(8,8)
for i in range(8):
for j in range(8):
K[i,j] = integrate(K_int[i,j], (r,-1,1), (s,-1,1))

simplify(K/K_factor)




Out[8]:

$\displaystyle \left[\begin{matrix}\frac{1}{2} - \frac{\nu}{6} & \frac{\nu}{8} + \frac{1}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{3 \nu}{8} - \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & - \frac{\nu}{8} - \frac{1}{8} & \frac{\nu}{6} & \frac{1}{8} - \frac{3 \nu}{8}\\\frac{\nu}{8} + \frac{1}{8} & \frac{1}{2} - \frac{\nu}{6} & \frac{1}{8} - \frac{3 \nu}{8} & \frac{\nu}{6} & - \frac{\nu}{8} - \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & \frac{3 \nu}{8} - \frac{1}{8} & - \frac{\nu}{12} - \frac{1}{4}\\- \frac{\nu}{12} - \frac{1}{4} & \frac{1}{8} - \frac{3 \nu}{8} & \frac{1}{2} - \frac{\nu}{6} & - \frac{\nu}{8} - \frac{1}{8} & \frac{\nu}{6} & \frac{3 \nu}{8} - \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & \frac{\nu}{8} + \frac{1}{8}\\\frac{3 \nu}{8} - \frac{1}{8} & \frac{\nu}{6} & - \frac{\nu}{8} - \frac{1}{8} & \frac{1}{2} - \frac{\nu}{6} & \frac{1}{8} - \frac{3 \nu}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{\nu}{8} + \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4}\\\frac{\nu}{12} - \frac{1}{4} & - \frac{\nu}{8} - \frac{1}{8} & \frac{\nu}{6} & \frac{1}{8} - \frac{3 \nu}{8} & \frac{1}{2} - \frac{\nu}{6} & \frac{\nu}{8} + \frac{1}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{3 \nu}{8} - \frac{1}{8}\\- \frac{\nu}{8} - \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & \frac{3 \nu}{8} - \frac{1}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{\nu}{8} + \frac{1}{8} & \frac{1}{2} - \frac{\nu}{6} & \frac{1}{8} - \frac{3 \nu}{8} & \frac{\nu}{6}\\\frac{\nu}{6} & \frac{3 \nu}{8} - \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & \frac{\nu}{8} + \frac{1}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{1}{8} - \frac{3 \nu}{8} & \frac{1}{2} - \frac{\nu}{6} & - \frac{\nu}{8} - \frac{1}{8}\\\frac{1}{8} - \frac{3 \nu}{8} & - \frac{\nu}{12} - \frac{1}{4} & \frac{\nu}{8} + \frac{1}{8} & \frac{\nu}{12} - \frac{1}{4} & \frac{3 \nu}{8} - \frac{1}{8} & \frac{\nu}{6} & - \frac{\nu}{8} - \frac{1}{8} & \frac{1}{2} - \frac{\nu}{6}\end{matrix}\right]$



We can check some numerical vales for $E=1$ Pa and $\nu=1/3$



In [9]:

K_num = K.subs([(E, 1), (nu, S(1)/3)])




In [10]:

K_num.eigenvects()




Out[10]:

$\displaystyle \left[ \left( 0, \ 3, \ \left[ \left[\begin{matrix}1\\0\\1\\1\\0\\1\\0\\0\end{matrix}\right], \ \left[\begin{matrix}1\\0\\1\\0\\1\\0\\1\\0\end{matrix}\right], \ \left[\begin{matrix}-1\\1\\-1\\0\\0\\0\\0\\1\end{matrix}\right]\right]\right), \ \left( \frac{1}{2}, \ 2, \ \left[ \left[\begin{matrix}-1\\0\\1\\0\\-1\\0\\1\\0\end{matrix}\right], \ \left[\begin{matrix}0\\-1\\0\\1\\0\\-1\\0\\1\end{matrix}\right]\right]\right), \ \left( \frac{3}{4}, \ 2, \ \left[ \left[\begin{matrix}0\\-1\\-1\\0\\0\\1\\1\\0\end{matrix}\right], \ \left[\begin{matrix}1\\0\\0\\-1\\-1\\0\\0\\1\end{matrix}\right]\right]\right), \ \left( \frac{3}{2}, \ 1, \ \left[ \left[\begin{matrix}-1\\-1\\1\\-1\\1\\1\\-1\\1\end{matrix}\right]\right]\right)\right]$




In [ ]: