Finite Elements Notebook
by Arthur B. Soprano - arthursoprano@gmail.com
Class notes from Professor Giuseppe Gambolati and Professor Carlo Janna: Solving Geomechanical Problems with the Finite Element Method. Special thanks to Jaime Ambrus and Edson Tadeu Manoel for helping with the class notes.
This notebook intends to give a brief review and a short example of the Finite Element Method - FEM using Python programming language. First, let's import the necessary packages. For our analysis, we'll need numpy and matplotlib. Also, we'll define the images to be showed in this notebook.
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Image
image1 = Image(filename='images/square_summable.jpg')
image2 = Image(filename='images/triangle_fe.jpg')
# Remove this line for a pop-up (qt) plot
%matplotlib inline
Indroduction
Variational Principles
We say that a variational principle exists if the solution to a (generally well behaved) problem is obtained through the minimization of
an integral expression = functional = variational principle
Tipically the wanted solution is subject to boundary conditions, so we are in fact looking for a constrained minimum.
Functional Minimization Methods
Minimization of the variational principle $\Omega$ is performed with the aid of a functional linear space. The items of a linear space are functions. The norm of a function $f$ is defined as (in a one-dimensional space):
$$ \left\lVert f \right\rVert = \sqrt{\int^{b}_{a}{f^{2}(x)}\, \mathrm{d}x} = \text{Euclidean "norm"} $$
If the above norm exists for any space fucntion the space is said to be measurable and the norm of $f$ is assumed to be the measure of $f$. The scalar product between $f_1$ and $f_2$ is defined as:
$$ \int^{b}_{a}{f_1(x)f_2(x)}\, \mathrm{d}x \neq \infty $$
if the space is measurable and if
$$ \int^{b}_{a}{f_1(x)f_2(x)}\, \mathrm{d}x = 0 $$
the functions $f_1$ and $f_2$ are said to be orthogonal.
The functional space made of all functions which have finite Euclidean norm is denoted by $L_2$. The functional space $L_2$ is made of all the functions which are square summable.
image1 # Examples of square summable functions
If we prescribe restrictions on the continuity of the derivatives of the space we get higher order spaces called Sobolev spaces. The Sobolev spaces are subspaces of $L_2$. For example, $W^{(1)}_2$ is the space made by all the functions whose first derivatives are also square summable
$$ \sqrt{\int^{b}_{a}{f^{2} + \left(\frac{\partial f}{\partial x}\right)^2}\, \mathrm{d}x} \neq \infty $$
Similarly, $W^{(2)}_2$ is the space made by all the functions whose first and second order derivatives are square summable
$$ \sqrt{\int^{b}_{a}{f^{2} + \left(\frac{\partial f}{\partial x}\right)^2 + \left(\frac{\partial^2 f}{\partial x^2}\right)^2}\, \mathrm{d}x} \neq \infty $$
Higher order spaces are subsets of lower order spaces, i.e., they are subspaces of lower spaces. The lowest possible space is $L_2$, which contains all the Sobolev spaces $W^{(1)}_2,W^{(2)}_2,\ldots,W^{(k)}_2$.
Approximation of a Function in $L_2$
We first identify a (generally) infinite set of functions $\xi_1, \xi_2, \ldots, \xi_n$ belonging to $L_2$ that defines a basis. For $\xi_1, \xi_2, \ldots, \xi_n$ to represent a basis, it is mandatory that for any given function $f$ of $L_2$ and any given (small) number $\epsilon$ the exists an $n$ (depending on $\epsilon$) and a set of coefficient $a_1, a_2, \ldots, a_n$ such that the linear combination $f_n$ given by
$$ f_n = a_1 \xi_1 + a_2 \xi_2 + \ldots + a_n \xi_n $$
differs from $f$ less than $\epsilon$, namely
$$ \left\lVert f - f_n \right\rVert < \epsilon $$
Notice that the set $\xi_1, \xi_2, \ldots, \xi_n$ defines a subspace $S_n$ of the full space $L_2$. It is intuitive that as $n$ increases $S_n$ approximates $L_2$ better and better with $S_n \rightarrow L_2$ as $n \rightarrow \infty$.
Ritz FEM Method
The Ritz FEM method looks for an approximation $\bar{u}_n$ to the solution of the variational problem (i.e., the minimizing function $\bar{u}$) with the aid of an $L_2$ basis:
$$ \bar{u}_n = a_1 \xi_1 + a_2 \xi_2 + \ldots + a_n \xi_n $$
where $\xi_1, \xi_2, \ldots, \xi_n$ are the basis functions. We consider the functional $\Omega(\bar{u}_n)$ in the space $S_n$ and minimize it in $S_n$:
$$ \frac{\partial \Omega}{\partial a_i}(\bar{u}_n) = 0 $$
for $i=1,2,\ldots,n$, giving rise to $n$ algebraic equations that are linear if $\Omega$ is a quadratic functional.
Example of FEM
In order to understand the Ritz FEM method, let's consider the following example
$$ \Omega(\bar{u}_n) = \int^{1}_{0}\left[\frac{1}{2}(u')^2 + \alpha u\right]\, \mathrm{d}x $$
subject to $u(0)=u(1)=0$ boundary conditions. Select a polinomial basis:
$$ 1, x, x^2, x^3, \ldots $$
Now let's look for the approximate solution
$$ \bar{u}_n = x(1-x) \sum^{n}_{i=1}a_i x^{i-1} $$
which satisfies a priori the boundary conditions. The lowest order approximation ($n = 1$) is
$$ u_1 = a_1 x (1-x) $$
Replace $u_1$ in $\Omega(\bar{u}_n)$ and differentiate with respect to $a_1$:
$$ \frac{\partial \Omega}{\partial a_i} = \frac{\partial}{\partial a_i} \int^1_0\left[\frac{1}{2} a^{2}_{1} (1-2x)^2 + \alpha a_1 x (1-x) \right] \, \mathrm{d}x = 0 $$ $$ \int^1_0\left[a_1 (1-2x)^2 + \alpha x (1-x) \right] \, \mathrm{d}x = 0 $$ $$ a_1 = -\frac{\alpha}{2} $$ $$ u_1 = -\frac{\alpha}{2}x(1-x) $$
Is is easily shown that with $n>1$ we would get $a_2=a_3=\ldots=a_n=0$, i.e.,
$$ u_1 = -\frac{\alpha}{2}x(1-x)=\text{analytical solution} $$
Select this other basis:
$$ a_n = \sum^n_{i=1} a_i sin(i \pi x) $$
satisfying the boundary conditions. Let's choose $u_1 = a_1 sin(\pi x)$ and minimize $\Omega(a_1)$:
$$ \frac{\partial \Omega}{\partial a_i} = \frac{\partial}{\partial a_i} \int^1_0\left[\frac{1}{2} a^{2}_{1} \pi^2 cos^2(\pi x) + \alpha a_1 sin(\pi x) \right] \, \mathrm{d}x = 0 $$ $$ \int^1_0\left[a_{1} \pi^2 cos^2(\pi x) + \alpha sin(\pi x) \right] \, \mathrm{d}x = 0 $$ $$ a_1 = -\frac{4 \alpha}{\pi^3} $$ $$ u_1 = -\frac{4 \alpha}{\pi^3}sin(\pi x) $$
We can compare the two approximations and check that they have similar values:
# Numerical example
x = np.array([0.0, 0.25, 0.5, 1.0])
alpha = 1.0
u1 = -0.5 * alpha * x * (1 - x)
u2 = -(4.0 * alpha / (np.pi ** 3)) * np.sin(np.pi*x)
print " x | u1 | u1 (sin())"
for c1, c2, c3 in zip(x, u1, u2):
print "%.2f | %.3f | %.3f" % (c1, c2, c3)
Major concepts to address when implementing FEM
- Selected basis $\rightarrow$ convergence speed $\rightarrow$ numerical efficiency of FEM.
Completeness of the selected basis.
If we had selected $$ u_n = \sum^n_{i=1} a_i sin(2 \pi i x)$$ we would have obtained $$a_1 = a_2 = \ldots=a_n=0 $$ The selected basis is incomplete and totally unable to describe the wanted solution.
Functional complexity from complex basis functions. Hence special care to be devoted to simple basis functions (in other words to simple finite elements, e.g., triangles in 2D and tetrahedrals in 3D).
Galerkin's Variational approach
FEM with Galerkin does not require the existence of a variational principle associated to a corresponding PDE of the type
$$ \mathbf{A} u(x) = f(x) $$
where $x=[x_1\ x_2\ \ldots\ x_n]$ and $\mathbf{A}$ is a differential operator, not necessarily linear. As with Ritz, Galerkin looks for an approximate solution
$$ u_n = \sum^n_{i=1} a_i \xi_i $$
where $\xi_1, \xi_2, \ldots, \xi_n$ are the basis functions. Replacing $u_n$ into the PDE provides the residual $\mathbf{E}$
$$ \mathbf{E} = \mathbf{A} u_n(x) - f(x) $$
Galerkin's coefficients $a_i$ for $i=1,2,\ldots,n$ are found by prescribing the orthogonality between $\mathbf{E}$ and $\xi_i(x)$:
$$ \int_R \mathbf{E} \xi_i(x) \, \mathrm{d}R = \int_R \left[\mathbf{A} u_n(x) - f(x)\right] \xi_i(x) \, \mathrm{d}R = 0 $$
for $i=1,2,\ldots,n$. This is the i-th algebraic Galerkin's equation.
Galerkin is more general than Ritz (as it does not require a variational principle) and can be shown to give the same FEM equations as Ritz if the PDE (i.e., the operator $\mathbf{A}$) is linear and symmetric positive definite, in which case the PDE is the Euler equation of an appropriate variational principle.
Galerkin's method may be easily generalized into the
Method of Weighted Residuals
$$ \int_R \underbrace{\left(\mathbf{A} u_n - f \right)}_\text{residual} \underbrace{w_i}_{\substack{\text{weight or}\\\text{test function}}} \, \mathrm{d}R = 0 $$
The method of Weighted Residuals will give rise to different methods depending on the choice of the weight function $w_i$:
Subdomain method
$$ w_i= \begin{cases} 1 & \text{if } x \in R_i \\ 0 & \text{if } x \notin R_i \\ \end{cases} $$
Collocation: $$w_i = \delta(x-x_i)$$
$$ \int_R \left(\mathbf{A} u_n - f\right) \delta(x-x_i) \, \mathrm{d}R = \left(\mathbf{A} u_n - f\right)_{x=x_i} = 0 $$
Least Squares: $$w_i = \frac{\partial \left(\mathbf{A} u_n - f\right)}{\partial a_i}$$ namely, $$ \int_R \left(\mathbf{A} u_n - f\right) \mathbf{A} \xi_i \, \mathrm{d}R = 0 $$
that is minimum of
$$ I = \int_R \left(\mathbf{A} u_n - f\right)^2 \, \mathrm{d}R = \text{minimum} $$
Solving Equilibrium Equations by Finite Element
Triangular Finite Elements
The usual basis functions for the triangular element is given by
$$ \xi_i = \frac{a_i + b_i x + c_i y}{2\Delta^e} $$
where $\Delta^e$ is the triangle area
$$ A = \left( \begin{matrix} 1 & 1 & 1 \\ x_i & x_j & x_k \\ y_i & y_j & y_k \end{matrix} \right) \left| det(A) \right| = 2 \Delta^e $$
and
$$ \begin{align} a_i &= x_j y_k - x_k y_j \\ b_i &= y_j - y_k \\ c_i &= x_k - x_j \end{align} $$
image2 # Example a triangular finite element
Note that for triangular elements, integrals can be performed in closed form (without using quadrature).
Cauchy Equilibrium Equations
$$ \begin{align} \frac{\partial \sigma_x}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} +\frac{\partial \tau_{xz}}{\partial z} + X &= 0 \\ \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \sigma_y}{\partial y} +\frac{\partial \tau_{yz}}{\partial z} + Y &= 0 \\ \frac{\partial \tau_{xz}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} +\frac{\partial \sigma_z}{\partial z} + Z &= 0 \\ \end{align} $$
Theorem of Virtual Work
A body is in equilibrium if and only if the virtual work of internal forces is equal to the virtual work of the external forces for any virtual displacement.
$$ \underbrace{\delta L^i}_{\substack{\text{Work of}\\\text{internal forces}}} = \underbrace{\delta L^e}_{\substack{\text{Work of}\\\text{external forces}}} $$
$$ \delta L^i = \iint_\Omega \underline{\sigma}^T \cdot \underline{\delta \epsilon}\, \mathrm{d}\Omega $$
$$ \delta L^e = \iint_\Omega \underline{X}^T \cdot \underline{\delta u}\, \mathrm{d}\Omega + \int_\Gamma \underline{t}^T \cdot \underline{\delta u}\, \mathrm{d}S $$
$$ \underline{u}^e = \left( \begin{matrix} u_e \\ v_e \end{matrix} \right) = \underbrace{ \left( \begin{matrix} \xi^e_i & 0 & \xi^e_j & 0 & \xi^e_k & 0 \\ 0 & \xi^e_i & 0 & \xi^e_j & 0 & \xi^e_k \end{matrix} \right)}_{\mathbf{N}^e} \underbrace{\left( \begin{matrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \end{matrix} \right)}_{\underline{u}^e_N} $$
$$ \underline{\epsilon}^e = \left( \begin{matrix} \frac{\partial}{\partial x} & 0 \\ 0 & \frac{\partial}{\partial y} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} \end{matrix} \right) \left( \begin{matrix} \xi^e_i & 0 & \xi^e_j & 0 & \xi^e_k & 0 \\ 0 & \xi^e_i & 0 & \xi^e_j & 0 & \xi^e_k \end{matrix} \right) \left( \begin{matrix} u_i \\ v_i \\ u_j \\ v_j \\ u_k \\ v_k \end{matrix} \right) $$
$$ \underline{\epsilon}^e = \underbrace{\frac{1}{2 \Delta^e} \left( \begin{matrix} b_i & 0 & b_j & 0 & b_k & 0 \\ 0 & c_i & 0 & c_j & 0 & c_k \\ b_i & c_i & b_j & c_j & b_k & c_k \end{matrix} \right)}_{\mathbf{B}^e} \underline{u}^e_N $$
$$ \underline{\epsilon}^e = \mathbf{B}^e \cdot \underline{u}^e_N $$
We need to link deformation to stresses. Using linear elasticity and focusing on the isotropic case the constitutive equation, in matrix notation, is given by
$$ \underline{\sigma} = \left( \begin{matrix} \sigma_x \\ \sigma_y \\ \tau_{xy} \end{matrix} \right) = \underbrace{\frac{E}{1 - \nu^2} \left( \begin{matrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac{1-\nu}{2} \end{matrix} \right)}_{\mathbf{D} = \text{Elastic Constitutive Matrix}} \left( \begin{matrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{matrix} \right) $$
$$ \underline{\sigma} = \mathbf{D}^e \cdot \underline{\epsilon} $$
$$ \delta L^i = \iint_{\Delta^e} \underline{\sigma}^T \cdot \underline{\delta \epsilon}\, \mathrm{d}A = \iint_{\Delta^e} \underline{\delta \epsilon}^T \cdot \underline{\sigma} \, \mathrm{d}A $$
$$ \underline{\delta \epsilon} = \mathbf{B}^e \cdot \delta \underline{u}^e_N $$
$$ \underline{\sigma} = \mathbf{D}^e \cdot \underline{\epsilon} = \mathbf{D}^e \mathbf{B}^e \underline{u}^e_N $$
$$ \delta L^i = (\underline{u}^e_N)^T \iint_{\Delta^e} (\mathbf{B}^e)^T \mathbf{D}^e \mathbf{B}^e \underline{u}^e_N \, \mathrm{d}A = (\underline{u}^e_N)^T \left[ \Delta^e (\mathbf{B}^e)^T \mathbf{D}^e \mathbf{B}^e \right] \underline{u}^e_N $$
Now we need to compute $\delta L^e$:
$$ \begin{align} \delta L^e &= \iint_{\Delta^e} \underline{X}^T \mathbf{N}^e \delta \underline{u}^e_N\, \mathrm{d}A + \int_{\partial \Delta^e} \underline{t}^T \mathbf{N}^e \delta \underline{u}^e_N \, \mathrm{d}S \\ &= (\delta \underline{u}^e_N)^T \underbrace{\left[ \iint_{\Delta^e} (\mathbf{N}^e)^T \underline{X} \, \mathrm{d}A + \int_{\partial \Delta^e} (\mathbf{N}^e)^T \underline{t}\, \mathrm{d}S \right]}_{\underline{f}^e_N} \end{align} $$
By making the works equal, we get
$$ (\delta \underline{u}^e_N)^T \cdot \Delta^e \left[ (\mathbf{B}^e)^T \mathbf{D}^e \mathbf{B}^e \right] \underline{u}^e_N = (\delta \underline{u}^e_N)^T \cdot \underline{f}^e_N $$
The linear system above is for each finite element and it results in a global linear system as follows
$$ \mathbf{K} \vec{U} = \vec{F} $$
where $\mathbf{K}$ is the stiffness matrix, $\vec{U}$ the unknowns and $\vec{F}$ the external forces.
Numerical Implementation
Let's load our data. We'll work with two example cases: a small case and a large case. The input data is rather large, hence it will be configured in an external Python script (fem_data.py). From this script, we'll import the class Setup and the methods GetSetupLargeCase and GetSetupSmallCase, as they will configure our two examples. Of course that, from their own names, we can infere that one corresponds to a large problem and the other one to a small one.
The setup includes constants and variables such as the number of degrees of freedom, the coordinates of each node, the number of elements and the conectivity array, as well as the boundary conditions and properties of our problem.
Some hipotesis that will be assumed are:
- All elements have the same number of nodes (triangles)
- Properties are the same for all elements
from fem_data import Setup, GetSetupLargeCase, GetSetupSmallCase
setup = Setup()
#GetSetupSmallCase(setup)
GetSetupLargeCase(setup)
number_of_dof = setup.NumGL # number of defrees of freedom
number_of_nodes_el = setup.NumNosEl
number_of_nodes = setup.NumNos
coordinates = setup.Coord
number_of_elements = setup.NumElem
connectivity = setup.Incid
material_properties = setup.PropMat
geometric_properties = setup.PropGeo
dirichlet_bc = setup.CCDirichlet
newmann_bc = setup.CCNewmann
body_force = setup.Fcorpo
The method below is an auxiliary method that will plot the results that we'll obtain from the FEM simulation, i.e, it will allow us to view the initial and final conditions of our nodal variables.
def ViewMesh(
coord,
connec,
hold=False,
show=True,
label='unknown',
node_color='r',
edge_color='k'
):
import networkx as nx
G = nx.Graph()
number_of_nodes = coord.shape[0]
for i,j,k in connec:
G.add_edge(i,j)
G.add_edge(i,k)
G.add_edge(j,k)
pos = {}
for n in xrange(number_of_nodes):
pos[n] = coord[n]
node_size = 20 + (1000 / number_of_nodes)
nx.draw(G, pos=pos, hold=hold, node_color=node_color,
label=label, edge_color=edge_color, node_size=node_size)
if show:
pyplot.show()
Element Stiffness Matrix
The following method returns the Stiffness Matrix for the linear triangle element in plane stress:
def KeT3_EPT(coord, mat_prop, geo_prop, body_force):
'''
Method that returns the Stiffness Matrix
for the linear triangle element in plane
stress.
'''
Ex = mat_prop[0] # Young's modulus
Nu = mat_prop[1] # Poisson's ratio
t = geo_prop
# Auxiliary functions
ci = - coord[1,0] + coord[2,0]
cj = - coord[2,0] + coord[0,0]
ck = - coord[0,0] + coord[1,0]
bi = + coord[1,1] - coord[2,1]
bj = + coord[2,1] - coord[0,1]
bk = + coord[0,1] - coord[1,1]
d = (ck * bj - cj * bk)
area = d / 2
B = (1 / d) * np.matrix([
[bi, 0, bj, 0, bk, 0],
[ 0, ci, 0, cj, 0, ck],
[ci, bi, cj, bj, ck, bk]
])
D = Ex / (1 - Nu ** 2) * np.matrix([
[ 1, Nu, 0],
[Nu, 1, 0],
[ 0, 0, (1 - Nu) / 2]
])
Ke = B.T * D * B * (area * t)
Fe1 = (area * body_force[0,1] * t
* np.array([1/3, 0.0, 1/3, 0.0, 1/3, 0.0]))
Fe2 = (area * body_force[1,1] * t
* np.array([0.0, 1/3, 0.0, 1/3, 0.0, 1/3]))
Fe = Fe1 + Fe2
return Ke, Fe
The code below will generate the matrix ID that enumerates each equation (for $x$ and $y$ direction) to be solved.
count = 0
ID = np.zeros((number_of_nodes, number_of_dof))
for n in xrange(number_of_nodes):
for dof in xrange(number_of_dof):
ID[n,dof] = count
count += 1
assert count == number_of_dof * number_of_nodes
Create the Linear System
With the algorithm below, we'll build the global matrix $\mathbf{K}$ and global vector $\vec{F}$.
global_size = number_of_dof * number_of_nodes
K = np.zeros((global_size,global_size))
F = np.zeros((global_size,1))
size_el = number_of_dof * number_of_nodes_el
for e in xrange(number_of_elements):
# Builds element LM vector and element coordinates
# LM has the equations associated to element e
LM = np.zeros(size_el)
coord_el = np.zeros((size_el,2))
z = 0
for j in xrange(number_of_nodes_el):
J = connectivity[e,j] # global node of element
for dof in xrange(number_of_dof):
LM[z] = ID[J,dof]
z += 1
for k in xrange(2):
coord_el[j,k] = coordinates[J,k]
# Calls method that calculates element
# matrix Ke and element load vector Fe
Ke, Fe = KeT3_EPT(
coord_el,
material_properties,
geometric_properties,
body_force
)
for i in xrange(size_el):
I = LM[i]
for j in xrange(size_el):
J = LM[j]
# Adds element matrix Ke to global matrix K
K[I,J] += Ke[i,j]
# Adds load vector Fe in F
F[I] = F[I] + Fe[i]
Boundary Conditions
Now we'll introduce the Newmann boundary condition to the load vector $\vec{F}$:
for node, dof, value in newmann_bc:
eq = ID[node,dof]
F[eq] = F[eq] + value
For the Dirichlet boundary conditions, the unkown has a prescribed value. In this case, we set the line of the matrix associated to that unknown to zero with the value one in the diagonal. Moreover, we subtract the product of the column of the matrix that corresponds to the precribed unknown time the prescribed value from the independent vector and make that column equal to zero (the size of the final linear system could be reduced).
for node, dof, value in dirichlet_bc:
eq = ID[node,dof]
F[:, 0] = F[:, 0] - value * K[:, eq]
# zeros the whole line and column
# associated to eq, except the
# diagonal that is equal to 1.0
K[eq,:] = 0.0
K[:,eq] = 0.0
K[eq,eq] = 1.0
F[eq] = value
Solve Linear System
With the global linear system complete with the boudnary condtiions, we can move on to solve the resulting linear system $\mathbf{K} \vec{U} = \vec{F}$ for the static problem:
U = np.linalg.solve(K, F)
With the resulting deformation vector $\vec{U}$ we can generate the deformed coordinates for the elements:
coord_def = np.zeros((number_of_nodes, 2))
Uf = np.zeros((number_of_nodes, number_of_dof))
b = 0
for n in xrange(number_of_nodes):
for dof in xrange(number_of_dof):
Uf[n,dof] = U[b]
b += 1
# Amplifying factor for vizualization purposes
factor = 100
for i in xrange(number_of_nodes):
coord_def[i,0] = coordinates[i,0] + factor * Uf[i,0]
coord_def[i,1] = coordinates[i,1] + factor * Uf[i,1]
View results
With our ViewMesh method defined previously, we can view the results from the linear system:
f = plt.figure()
ViewMesh(coord_def, connectivity, hold=True, show=False,
label='Deformed', node_color='r', edge_color='k')
ViewMesh(coordinates, connectivity, hold=True, show=False,
label='Original', node_color='b', edge_color='0.5')
ax = plt.gca()
handles, labels = ax.get_legend_handles_labels()
plt.legend(handles=[handles[0], handles[2]],bbox_to_anchor=(0.8, 1),
bbox_transform=plt.gcf().transFigure)
f.set_size_inches(10,8)
plt.show()
Element Stress
The following method returns the stress for the linear triangle element in plane stress. The principal stress are obtained by calculating the eigenvalues of the element stress matrix.
def StressT3_EPT(coord, UEl, mat_prop, geo_prop):
'''
Method that returns the Stress for the
linear triangle element in plane stress.
'''
Ex = mat_prop[0] # Young's modulus
Nu = mat_prop[1] # Poisson's ratio
t = geo_prop
# Auxiliary functions
ci = - coord[1,0] + coord[2,0]
cj = - coord[2,0] + coord[0,0]
ck = - coord[0,0] + coord[1,0]
bi = + coord[1,1] - coord[2,1]
bj = + coord[2,1] - coord[0,1]
bk = + coord[0,1] - coord[1,1]
d = (ck * bj - cj * bk)
area = d / 2
B = (1 / d) * np.matrix([
[bi, 0, bj, 0, bk, 0],
[ 0, ci, 0, cj, 0, ck],
[ci, bi, cj, bj, ck, bk]
])
D = Ex / (1 - Nu ** 2) * np.matrix([
[ 1, Nu, 0],
[Nu, 1, 0],
[ 0, 0, (1 - Nu) / 2]
])
StressE = D * B * UEl.T
# Create 2D stress matrix
S = np.zeros((number_of_dof, number_of_dof))
S[0,0] = StressE[0]
S[1,1] = StressE[1]
S[1,0] = S[0,1] = StressE[2]
# Stress in z direction is ZERO
# under plane stress conditions
StressE = np.append(StressE, [[0]], axis=0)
# Calculate principal stresses by
# finding the eigenvalues of S
w, v = np.linalg.eig(S)
s1, s2 = w
# Von Mises Stress
SVME = np.sqrt(s1 ** 2 - s1*s2 + s2 ** 2)
# Tresca Stress
STRE = 0.5 * np.abs(s1 - s2)
return StressE, SVME, STRE
Stress calculation
Calulates the stress and performs averaging.
count = np.zeros(number_of_nodes)
Sx = np.zeros(number_of_nodes)
Sy = np.zeros(number_of_nodes)
Sxy = np.zeros(number_of_nodes)
Sz = np.zeros(number_of_nodes)
SVM = np.zeros(number_of_nodes)
STR = np.zeros(number_of_nodes)
for e in xrange(number_of_elements):
# Builds element LM vector and element coordinates
# and element displacement
LM = np.zeros(size_el)
UEl = np.zeros(size_el)
coord_el = np.zeros((size_el,2))
z = 0
for j in xrange(number_of_nodes_el):
J = connectivity[e,j]
for i in xrange(number_of_dof):
LM[z] = ID[J,i]
UEl[z] = Uf[J,i]
z += 1
for k in xrange(2):
coord_el[j,k] = coordinates[J,k]
# Calls method that calculates element
# stress StressE, von Mises stress SVME
# and Tresca stress STRE
StressE, SVME, STRE = StressT3_EPT(
coord_el,
np.matrix(UEl),
material_properties,
geometric_properties
)
i, j, k = connectivity[e]
Sx[i] = Sx[i] + StressE[0]
Sx[j] = Sx[j] + StressE[0]
Sx[k] = Sx[k] + StressE[0]
Sy[i] = Sy[i] + StressE[1]
Sy[j] = Sy[j] + StressE[1]
Sy[k] = Sy[k] + StressE[1]
Sxy[i] = Sxy[i] + StressE[2]
Sxy[j] = Sxy[j] + StressE[2]
Sxy[k] = Sxy[k] + StressE[2]
Sz[i] = Sz[i] + StressE[3]
Sz[j] = Sz[j] + StressE[3]
Sz[k] = Sz[k] + StressE[3]
SVM[i] = SVM[i] + SVME
SVM[j] = SVM[j] + SVME
SVM[k] = SVM[k] + SVME
STR[i] = STR[i] + STRE
STR[j] = STR[j] + STRE
STR[k] = STR[k] + STRE
count[i] += 1
count[j] += 1
count[k] += 1
# Divide by counters
for i in xrange(number_of_nodes):
Sx[i] = Sx[i] /count[i]
Sy[i] = Sy[i] /count[i]
Sxy[i] = Sxy[i]/count[i]
Sz[i] = Sz[i] /count[i]
SVM[i] = SVM[i]/count[i]
STR[i] = STR[i]/count[i]
Stress Field
Let's view the different stress fields, such as, the $\sigma_x$ and $\sigma_y$ as well as the shear stress $\tau_{xy}$, the von Mises and Tresca stress.
# Only works for the LARGE CASE
x = coordinates[:,0]
y = coordinates[:,1]
# Choose with stress you want to vizualize:
# Sx, Sy, Sxy, SVM or STR
f, axes = plt.subplots(2, 2, sharex='col', sharey='row')
((ax1, ax2), (ax3, ax4)) = axes
vmin_Sx = min(Sx)
vmin_Sy = min(Sy)
vmin_Sxy = min(Sxy)
vmin_SVM = min(SVM)
vmin_STR = min(STR)
vmin = min([vmin_Sx, vmin_Sy, vmin_Sxy, vmin_SVM, vmin_STR])
vmax_Sx = max(Sx)
vmax_Sy = max(Sy)
vmax_Sxy = max(Sxy)
vmax_SVM = max(SVM)
vmax_STR = max(STR)
vmax = max([vmax_Sx, vmax_Sy, vmax_Sxy, vmax_SVM, vmax_STR])
linewidths = 0.2
n_layers = 100
ax1.tricontour(x, y, connectivity, Sx, 15, linewidths=linewidths, colors='k')
im1 = ax1.tricontourf(x, y, connectivity, Sx, n_layers, cmap=plt.cm.rainbow,
vmin=vmin, vmax=vmax)
ax1.set_title(r'$\sigma_x$' + ' Stress Field')
ax2.tricontour(x, y, connectivity, Sy, 15, linewidths=linewidths, colors='k')
im2 = ax2.tricontourf(x, y, connectivity, Sy, n_layers, cmap=plt.cm.rainbow,
vmin=vmin, vmax=vmax)
ax2.set_title(r'$\sigma_y$' + ' Stress Field')
ax3.tricontour(x, y, connectivity, Sxy, 15, linewidths=linewidths, colors='k')
im3 = ax3.tricontourf(x, y, connectivity, Sxy, n_layers, cmap=plt.cm.rainbow,
vmin=vmin, vmax=vmax)
ax3.set_title(r'$\tau_{xy}$' + ' Stress Field')
ax4.tricontour(x, y, connectivity, SVM, 15, linewidths=linewidths, colors='k')
im4 = ax4.tricontourf(x, y, connectivity, SVM, n_layers, cmap=plt.cm.rainbow,
vmin=vmin, vmax=vmax)
ax4.set_title('Von Mises Stress Field')
#ax4.tricontour(x, y, connectivity, STR, 15, linewidths=linewidths, colors='k')
#im4 = ax4.tricontourf(x, y, connectivity, STR, n_layers, cmap=plt.cm.rainbow,
# vmin=vmin, vmax=vmax)
#ax4.set_title('Tresca Stress Field')
f.subplots_adjust(right=0.8)
cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])
f.colorbar(im2, cax=cbar_ax)
f.set_size_inches(10,8)
plt.show()
The source code of this notebook can be found on github.