Finite Elements Lecture 2



In [1]:

    
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())









    Out[1]:

Georges Limbert

Aims

Present the general procedure to convert a BVP into its discretised version

Outcomes

Derive variational form of BVP
Discretise the associated field variable using shape functions and derive discrete version of the BVP
Write a procedure to solve a simple 2D elasticity problem

FE in a nutshell

Transform PDE into a variational problem
Introduce piecewise approximation to variables
Discreise the physical domain into elements; write approximate equations for each elements (meshing). Gives matrix equations.
Assemble each into global system to solve

Weak form

Multiply by test function; convert to integral form. Key that test function vanishes on the boundary of the domain, and that it's arbitrary.

Field interpolation

In direct stiffness approach, derived the stiffness matrix of each element using elementary mechanics of solids (Hooke's law). General approach:

Approx displcement in each elemnt
Approx strain and stress from displacement
Derive stiffness matrix
Derive global matrix
Solve.

Trial solution

Guess a trial form for the solution (finite dimensional)

$$ u(x) = \sum a_i \phi_i(x) $$

where $\phi_i$ are the basis functions.

Simplest: basis elements are linear polynomials.

Obtain $a_i$ coefficients from nodal displacements; general form is

$$ w(x) = N^T {\bf d} $$

where ${\bf d}$ are the displacements and $N$ is the matrix of shape functions.

Shape functions

Properties:

Kronecker delta property: Has value 1 at one node and zero at all other nodes.
Compatibility: displacement approximation is continuous across element boundaries
Completeness: Sum of shape functions is $1$ so that $N_1(x) x_1 + N_2(x) x_2 = x$ for all $x$.

Derive strain and stress

Continuous strain

$$ \varepsilon = \frac{d w}{dx} $$

becomes in the discrete version

$$ \varepsilon = \frac{d N^T {\bf d}}{dx} = \frac{dN_1}{dx} \frac{dN_2}{dx} {\bf d} = B {\bf d} $$

The strain-displacement matrix $B$ is

$$ B = \frac{d N^T}{dx} $$

In 2D case have simple result that $B$ is independent of $x$ (depends only on the nodal locations) so the strain is constant within a linearly interpolated displacement-based element (this follows for stress as well). This highlights the importance of the interpolation scheme, and shows why in some cases it's important to split the parts of the response into volumetric and deviatric pieces.

1D Hooke's law is

$$ \sigma = E \varepsilon = E \frac{dw}{dx} $$

which generalizes to give the stress

$$ \sigma = E B {\bf d}. $$

Generally stress and strain approximations (esp. for Lagrange polynomials) are discontinuous across element boundaries.

Small 3d detour

Weak form links internal and external virtual work. The constitutive equations (eg linear elasticity) link strain to stress via

$$ \sigma = \mathbb{C} : \varepsilon $$

The 2D case converts the tensor form to a matrix equation; so the symmetry of the stress and strain can be used to convert the 2-tensor into a vector containing only the symmetric terms, and the 4-tensor into a simple matrix.

Using that the strain is the symmetric part of the gradient of the virtual displacement can relate the displacement, virtual displacement and elasticity tensor $\mathbb{C}$.

Following this through to the end we'll get a weak form applied to the shape functions.

Stiffness matrix

The matrix form of this final weak form can be written (at the element level)

$$ K_{AiBk} u_i^A - F^A_i $$

where $K$ is the stiffness matrix, $F$ the nodal foce vector, $A, B, \dots$ are node indices, and $i, j, \dots$ are coordinate indices.

In tensor form

$$ K = \int B \mathbb{C} B^T, \qquad F = - \left( \int b N^a \, dv + \int_{\partial \Omega} \bar{t} N^A \right). $$

In 1d this reduces to the direct stiffness method of the previous lecture.

Example for the lab

Go from the FE mesh to the actual element (in the physical domain), then to the parent element (in the parametric domain) by changing to parametric coordinates. In the parent element the shape functions have very simple forms.

Next step is to interpolate the position vector; ${\bf x}$ will just be a linear combination of the shape functions.

After this you can interpolate the displacement vector in the same way (this is the isoparametric approach; same shape function for both displacement and position).

This leads to the shape function vector written in terms of the parametric coordinates. When taking derivatives you have to link to the physical coordinates. So

$$ \frac{\partial N_I}{\partial x} = \frac{\partial N}{\partial \xi} \left( \frac{\partial x}{\partial \xi} \right)^{-1} $$

where the final term is the (inverse of the) Jacobian. The Jacobian of the geometric transformation from the partent to the actual element must have positive determinant (non-zero or the transformation fails; positive to retain orientation).

As the shape functions have simple, fixed forms in parametric coordinates the $\partial N / \partial \xi$ also have simple, fixed forms.

$B$ matrix and 2D strain tensor

Multiple ways of writing the matrix out; careful consideration of order important as it makes the matrix solution more/less efficient.

Integrals

Standard approach: Gauss integration. Take outer product of 1D stencils, use standard Gauss rules to compute integrals inside the elements.

The transformation to the physical domain uses the same transformation (Jacobian) matrix; the Gauss integration formula will pick up a contribution from the Jacobian.

Assembly

Have node coordinates array and the IEN (element nodes) array. The latter is a matrix connecting the nodes to the elements (connectiveity matrix) Rows represent the elements and columns represent the nodes that support the element.

Displacement prescribed at specific nodes (some of the boundary).

The global system matrix formed by assembling the individual element siffness matrices according to the element topology. By iterating over each element, then iterating over each Gauss point.

Call Gauss quad points
Call shape functions and derivatives
Call J matrix and physical derivs
Form strain-disp. matrix $B$
Form element stiffness matrix $k$

At the end, assemble the global $K$ matrix by placing the stiffness contributions at the DOF numbers according to the IEN array.

Global matrix

Can split out those DOFs which are fee (active) and prescribed. This gives a smaller matrix.