In [1]:
from __future__ import division
from sympy.utilities.codegen import codegen
from sympy import *
from sympy import init_printing
from IPython.display import Image
init_printing()
In [2]:
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')
The constitutive model tensor in Voigt notation (plane strain) is $$C = \frac{(1 - \nu) E}{(1 - 2\nu) (1 + \nu) } \begin{pmatrix} 1 & \frac{\nu}{1-\nu} & 0\\ \frac{\nu}{1-\nu} & 1 & 0\\ 0 & 0 & \frac{1 - 2\nu}{2(1 - \nu)} \end{pmatrix}$$
But for simplicity we are going to use $$\hat{C} = \frac{C (1 - 2\nu) (1 + \nu)}{E} = \begin{pmatrix} 1-\nu & \nu & 0\\ \nu & 1-\nu & 0\\ 0 & 0 & \frac{1 - 2\nu}{2} \end{pmatrix} \enspace ,$$ since we can always multiply by that factor afterwards to obtain the correct stiffness matrices.
In [3]:
C = Matrix([[1 - nu, nu, 0],
[nu, 1 - nu, 0],
[0, 0, (1 - 2*nu)/2]])
C_factor = E/(1-2*nu)/(1 + nu)
C
Out[3]:
The enumeration that we are using for the elements is shown below
In [4]:
Image(filename='../img_src/4node_element_enumeration.png', width=300)
Out[4]:
What leads to the following shape functions
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]:
Thus, the interpolation matrix renders
In [6]:
H = zeros(2,8)
for i in range(4):
H[0, 2*i] = N[i]
H[1, 2*i + 1] = N[i]
H.T # Transpose of the interpolation matrix
Out[6]:
And the mass matrix integrand is $$M_\text{int}=H^TH$$
In [7]:
M_int = H.T*H
In [8]:
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[8]:
Being the stiffness matrix integrand $$K_\text{int} = B^T C B$$
In [9]:
K_int = B.T*C*B
The mass matrix is obtained integrating the product of the interpolator matrix with itself, i.e. $$\begin{align*} M &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} M_\text{int} dr\, ds\\ &= \int\limits_{-1}^{1}\int\limits_{-1}^{1} H^T H\, dr\, ds \enspace . \end{align*}$$
In [10]:
M = zeros(8,8)
for i in range(8):
for j in range(8):
M[i,j] = rho*integrate(M_int[i,j],(r,-1,1), (s,-1,1))
M
Out[10]:
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 [11]:
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))
K
Out[11]:
We can generate automatically code for Fortran, C or Octave/Matlab, although it will be useful just for non-distorted elements.
In [12]:
K_local = MatrixSymbol('K_local', 8, 8)
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "f95")
print code[0][1]
In [13]:
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "C")
print code[0][1]
In [14]:
code = codegen(("local_stiff", Eq(K_local, simplify(K))), "Octave")
print code[0][1]
We can check some numerical vales for $E=8/3$ Pa, $\nu=1/3$ and $\rho=1$ kg\m$^3$, where we can multiply by the factor that we took away from the stiffness tensor
In [15]:
(C_factor*K).subs([(E, S(8)/3), (nu, S(1)/3)])
Out[15]:
In [16]:
M.subs(rho, 1)
Out[16]:
As stated before, the analytic expressions for the mass and stiffness matrices is useful for non-distorted elements. In the general case, a mapping between distorted elements and these canonical elements is used to simplify the integration domain. When this transformation is done, the functions to be integrated are more convoluted and we should use numerical integration like Gauss-Legendre quadrature.
The Gauss-Legendre quadrature approximates the integral: $$ \int_{-1}^1 f(x)\,dx \approx \sum_{i=1}^n w_i f(x_i)$$
The nodes $x_i$ of an order $n$ quadrature rule are the roots of $P_n$ and the weights $w_i$ are given by: $$w_i = \frac{2}{\left(1-x_i^2\right) \left(P'_n(x_i)\right)^2}$$
For the first four orders, the weights and nodes are
In [17]:
wts = [[2], [1,1], [S(5)/9, S(8)/9, S(5)/9],
[(18+sqrt(30))/36,(18+sqrt(30))/36, (18-sqrt(30))/36, (18-sqrt(30))/36]
]
pts = [[0], [-sqrt(S(1)/3), sqrt(S(1)/3)],
[-sqrt(S(3)/5), 0, sqrt(S(3)/5)],
[-sqrt(S(3)/7 - S(2)/7*sqrt(S(6)/5)), sqrt(S(3)/7 - S(2)/7*sqrt(S(6)/5)),
-sqrt(S(3)/7 + S(2)/7*sqrt(S(6)/5)), sqrt(S(3)/7 + S(2)/7*sqrt(S(6)/5))]]
And the numerical integral is computed as
In [18]:
def stiff_num(n):
"""Compute the stiffness matrix using Gauss quadrature
Parameters
----------
n : int
Order of the polynomial.
"""
if n>4:
raise Exception("Number of points not valid")
K_num = zeros(8,8)
for x_i, w_i in zip(pts[n-1], wts[n-1]):
for y_j, w_j in zip(pts[n-1], wts[n-1]):
K_num = K_num + w_i*w_j*K_int.subs([(r,x_i), (s,y_j)])
return simplify(K_num)
In [19]:
K_num = stiff_num(3)
K_num - K
Out[19]:
A best approach is to use Python built-in functions for computing the Gauss-Legendre nodes and weights
In [20]:
from sympy.integrals.quadrature import gauss_legendre
In [21]:
x, w = gauss_legendre(5, 15)
print x
print w
In [22]:
from IPython.core.display import HTML
def css_styling():
styles = open('../styles/custom_barba.css', 'r').read()
return HTML(styles)
css_styling()
Out[22]: