Tuesday, May 10, 2011

Video Tutorial : How To Implement Finite Element Method 3

From this post we saw how to formulate one row of the stiffness matrix. But this is not how you want the computer to calculate it. In this post we are going to explore, through some animations, how to calculate the entire global stiffness matrix by assembling local or elementary matrices. Please make sure to read the instructions at the end to view the animations full screen and to have play and pause control.

Think Globally Act Locally
Ultimately you would want to program your approximate solution in some computing language. The first step is to convert the global system (often a large system) into a local framework. This is exactly, think globally act locally kind of a thing. To further elucidate, let me give you another quote:

the longest journey begins with a single step
Let me further explain the transformation with the help of the following animation:
Photobucket FEM Video Tutorial

Element Stiffness Matrix Calculation
Now that we know how to isolate a single element and what are the elemental conventions, we could calculate the element stiffness matrix. The animation below shows us how:
Photobucket FEM Video Tutorial

You must have noticed that in order to simplify the calculations, I assumed that $EA(x)=1$. The intuition behind this simplification is the fact that $EA$ is a material property and if that is equal to a constant then the stiffness matrix entries are only dependent upon the shape functions. This makes it easy to think about other problems, without knowing a lot about their physics. For example for the one dimensional electrostatic boundary value problem, the element stiffness matrix is given as:
$$k_{ij}=\int_{x_{1}^{l}}^{x_{2}^{l}} \varepsilon^{l}(x){N^{l}_{j}}^{\prime}{N^{l}_{i}}^{\prime} dx $$

Does it not look similar. It sure does. The only difference between the above expression and the one we saw in the animation is that $EA(x)$ is now replaced by the permittivity of the medium $\varepsilon (x)$, where the former is the product of the Young's modulus of elasticity $E$ and the area of the bar $A$ and the latter is the product of dielectric constant of the medium $\varepsilon_{r}$ and permittivity of the free space $\varepsilon_{o}$ ($\varepsilon=\varepsilon_{r}\varepsilon_{o}$). So for a case where $\varepsilon(x)=1$, there is no difference between the element stiffness matrices, provided the element sizes $h^{l}$ are the same.

Assembling the Element Matrices
The next step is to assemble all the local matrices in the global system. The first thing we need for this is the element connectivity table. An element connectivity table tells us how different elements are connected and which nodes connect which elements. For a one dimensional case, every element is linear and is connected by two nodes on either side. In turn each node connects two elements on either side, except for the first and the last node. To further explain this, the first node connects the first element only, the first element is connected by the first and the second node, and the second node connects the first and the second element. This means that the second node is common to both the first and the second element. Same is true for all the elements until the last element arrives. The last node only connects the last element and that is it.

Once we have the connectivity table, properly identifying the inner nodes and the boundary nodes and all the elements therein, we can move on to assemble the local matrices into a global one. What better way to do that than to show you another small animation.
Photobucket FEM Video Tutorial
Each row has three non zero entries which is coherent with what was said about the global stiffness matrix for a simple one dimensional case with linear elements. However, the first and the last rows contain only two non-zero entries. The reason is that the first and the last nodes connect to one element only.

Instructions for Viewing the Animations Full Screen
The instructions to watch the animation are repeated here from this post. Here is how to play this animation full screen. Right click and save file as a gif file. If you do not already have storm codec Media Player Classic then download it using this link. Open the saved gif file with this media player and play and pause as you like. If you do not play this gif file with media player classic then you can not pause the frames and some frames might be too fast to read properly. This animation could also be played with Quicktime Player. You will find it here.

See you next post.

4 comments:

  1. Change Notification

    I just changed the last animation. The new animation is clearer and more understandable.

    ReplyDelete
  2. Respected Sir, I am unable to save the files.
    jeetender

    ReplyDelete
  3. Just right click on the gifs and select save image as, that should solve the problem

    ReplyDelete
  4. Hey, it is really a nice tutorial. Could you help to explain how to set boundary condition? Especially the first row below.
    F = F - K(:,1) * u_o;
    F (1) = u_o;
    F (end) = F (end) + P;

    K(1,:) = 0;
    K(:,1) = 0;
    K(1,1) = 1;

    ReplyDelete

Follow Husnain_Inayat on Twitter View Husnain Inayat Hussain's LinkedIn profile