The displacement based finite elment stiffness matrix can be written like:
$${K^{QP}} = \int\limits_V {B_{ij}^Q{C_{ijkl}}B_{kl}^PdV} $$where:
$$C = \frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}}\left[ {\begin{array}{*{20}{c}} 1&{\frac{\nu }{{1 - \nu }}}&0\\ {\frac{\nu }{{1 - \nu }}}&1&0\\ 0&0&{\frac{{1 - 2\nu }}{{2(1 - \nu )}}} \end{array}} \right]$$is the elstic tensor and $B_{ij}^Q$ is the contribution to the strain-displacement interpolator from the $Q$ degree of freedom.
The 4-noded perfectly square element of side $2h$ is shown in the figure below:
with each node having 2 degrees of freedom corresponding to the rectangular components of the displacement vector.
In [4]:
%matplotlib notebook
from __future__ import division
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt
from IPython.display import Image
In the displacement based finite element method we assume that the displacemnts vector at any point $x_i$ over the element are expressed, via interpolation, in terms of the nodal displacements like:
$${u_i} = N_i^Q(r){{\hat u}^Q}$$and where $N_i^K(r)$ is the shape function associated to the $k$-th degree of freedom.
For the computation of the stiffness matrix ${K^{QP}}$ we require the term $B_{ij}$ relating nodal displacements to strains. This interpolator can be written in terms of drivativs of the shape functions as follows:
In the implementation this operator is computed by the subroutine strain-displacement matrix for a 4-noded element stdm4() as listed blow:
In [5]:
def shape4(x , y , h):
"""Shape functions for the bi-lineal element.
Parameters
----------
x , y: Space variables.
h : Element halwidth.
Returns
-------
N : Array
"""
N=sym.zeros(4)
N = 1/(4*h**2)*sym.Matrix([(h + x)*(h + y),
(h - x)*(h + y),
(h - x)*(h - y),
(h + x)*(h - y)])
return N
In [6]:
def stdm4(x , y , h ):
"""Strain-displacement interpolator for the bi-lineal element.
Parameters
----------
x , y: Space variables.
h : Element halwidth.
Returns
-------
B : Array
"""
B = sym.zeros(3,8)
N = shape4(x , y , h)
dhdx=sym.zeros(2,4)
for i in range(4):
dhdx[0,i]=sym.diff(N[i],x)
dhdx[1,i]=sym.diff(N[i],y)
#
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]
#
return B
We also need to evaluate the elastic constitutive matrix:
$$C = \frac{{E(1 - \nu )}}{{(1 + \nu )(1 - 2\nu )}}\left[ {\begin{array}{*{20}{c}} 1&{\frac{\nu }{{1 - \nu }}}&0\\ {\frac{\nu }{{1 - \nu }}}&1&0\\ 0&0&{\frac{{1 - 2\nu }}{{2(1 - \nu )}}} \end{array}} \right]$$by the subroutine umat().
In [7]:
def umat(nu,E):
"""Plane stress constitutive tensor.
Parameters
----------
nu : Posisson ratio.
E : Young modulus.
Returns
-------
C : Array. Constitutive matrix.
"""
#
C=sym.zeros(3,3)
G=E/(1-nu**2)
mnu=(1-nu)/2.0
C[0,0]=G
C[0,1]=nu*G
C[1,0]=C[0,1]
C[1,1]=G
C[2,2]=G*mnu
#
return C
In [8]:
C = sym.zeros(3,3)
B = sym.zeros(3,8)
K = sym.zeros(8,8)
In [9]:
x , y = sym.symbols('x y')
nu, E = sym.symbols('nu, E')
h = sym.symbols('h')
In [10]:
C = umat(nu,E)
B = stdm4(x , y , h)
K_int = B.T*C*B
In [11]:
nuu = 1.0/3.0
EE = 8.0/3.0
In [12]:
for i in range(8):
for j in range(8):
K[i,j] = sym.integrate(K_int[i,j], (x,-h,h), (y,-h,h))
kk=K.subs([(E, EE), (nu, nuu), (h, 1.00)])
In [18]:
print(sym.N(kk , 3))
In [19]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom_barba.css', 'r').read()
return HTML(styles)
css_styling()
Out[19]:
In [ ]: