In the previous post I discussed some fundamental concepts in implementing the finite element method. In today's post we will move further. The post ends with a small animation on the calculation of the stiffness matrix.
Stiffness Matrix
For those of you who have a background in engineering, probably know Hooke's Law. It states that the restoring force in the direction opposite and equal to the direction of the applied force is directly proportional to the displacement of the material under the applied load. Mathematically this is stated as:
The negative sign simply says that the restoring force is in the opposite direction of the applied force. The constant of proportionality $k$ is called the stiffness of the material. Its magnitude is easily calculated by:
$$k=\frac{f}{u}\ $$
Now why are we discussing it here in the context of finite elements. Recall from this post that the weak formulation is given by:
$$ \int_{0}^{l}EA(x) \frac{du}{dx}\ \frac{dv}{dx}\ dx = \int_{0}^{l}f(x)v(x)dx + \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$
In the previous post we also looked at the formulation of the basis functions and how these basis functions are used in the approximation of finite element method. From the same post let us now write the following:
$$u(x)=\sum_{j=1}^{n}u_{j}\phi_{j}$$
$$v(x)=\phi_{i}$$
This substitution is proposed by Galerkin and hence this kind of finite element method is also known as Galerkin Finite Element Method. Now what is special about this Galerkin method? The special thing about this is that the trial and the test function are defined using the same basis functions. Now remember we were trying to figure out how the term stiffness is involved in finite elements. Let us now substitute the test and the trial functions in the weak form:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx $$
$$+ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$
You may notice that we did not replace $u(x)$ in the term $$ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$. The reason is that this term is the boundary term and is only calculated at $x=0$ and $x=l$. For an overhanging bar (fixed at one end and free at the other) this term is zero at $x=0$ and is equal to $P\phi_{j}$ at $x=l$, where $P$ is the applied load at the free end. Substituting this results in:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$
Writing the above equation in matrix form gives the following:
$$K_{ij}u_{j}=F_{i}\qquad(a)$$
Where $$F_{i}=\int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$ is everything on the right hand side and
$$K_{ij}=\int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx$$.
Now since equation (a) resembles the Hooke's Law, hence the matrix $K$ is called the stiffness matrix. This is the mystery behind calling it the stiffness matrix and it remains as such even if Poisson's Equation is resolved for electric field between two conducting plates.
Now that we have an expression of the stiffness matrix, let us see a small animation which explains how to calculate it.
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.
This is continued until the next post.
$f=-ku$
The negative sign simply says that the restoring force is in the opposite direction of the applied force. The constant of proportionality $k$ is called the stiffness of the material. Its magnitude is easily calculated by:
$$k=\frac{f}{u}\ $$
Now why are we discussing it here in the context of finite elements. Recall from this post that the weak formulation is given by:
$$ \int_{0}^{l}EA(x) \frac{du}{dx}\ \frac{dv}{dx}\ dx = \int_{0}^{l}f(x)v(x)dx + \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$
In the previous post we also looked at the formulation of the basis functions and how these basis functions are used in the approximation of finite element method. From the same post let us now write the following:
$$u(x)=\sum_{j=1}^{n}u_{j}\phi_{j}$$
$$v(x)=\phi_{i}$$
This substitution is proposed by Galerkin and hence this kind of finite element method is also known as Galerkin Finite Element Method. Now what is special about this Galerkin method? The special thing about this is that the trial and the test function are defined using the same basis functions. Now remember we were trying to figure out how the term stiffness is involved in finite elements. Let us now substitute the test and the trial functions in the weak form:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx $$
$$+ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$
You may notice that we did not replace $u(x)$ in the term $$ \left. EA(x) \frac{du}{dx}\ v(x) \right|_{0}^{l}$$. The reason is that this term is the boundary term and is only calculated at $x=0$ and $x=l$. For an overhanging bar (fixed at one end and free at the other) this term is zero at $x=0$ and is equal to $P\phi_{j}$ at $x=l$, where $P$ is the applied load at the free end. Substituting this results in:
$$ \sum_{j=1}^{n}u_{j} \int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx= \int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$
Writing the above equation in matrix form gives the following:
$$K_{ij}u_{j}=F_{i}\qquad(a)$$
Where $$F_{i}=\int_{0}^{l}f(x)\phi_{i}dx + \left. P \phi_{i} \right|_{x=l}$$ is everything on the right hand side and
$$K_{ij}=\int_{0}^{l}EA(x)\frac{d\phi_{j}}{dx}\ \frac{d\phi_{i}}{dx}\ dx$$.
Now since equation (a) resembles the Hooke's Law, hence the matrix $K$ is called the stiffness matrix. This is the mystery behind calling it the stiffness matrix and it remains as such even if Poisson's Equation is resolved for electric field between two conducting plates.
Stiffness Matrix CalculationLet us now see how the stiffness matrix is calculated. This is facilitated with the help of a small animation. The directions to play the animation full screen are given just after.
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.
Q1. Why is the second integral zero in the $i-1^{st}$ column?
A. The second integral as asked in the animation is zero because the function $\phi_{i+1}$ is zero in the range $x_{i-1}\leq x \leq x_{i}$. Pause at the frame when this question is posed in the animation. You could only pause by playing the gif file in media player classic.Q2. Why is the first integral zero in the $i+1^{st}$ column?
A. The first integral as asked in the animation is zero because the function $\phi_{i-1}$ is zero in the range $x_{i}\leq x \leq x_{i+1}$. Pause at the frame when this question is posed in the animation. You could only pause by playing the gif file in media player classic.This is continued until the next post.
No comments:
Post a Comment