Benutzer:Dirk Huenniger/fem

A simple example

Consider the spring on the right. Let the deformation of the spring x be zero for a mass m of zero. If you know try various masses you will find that the deformation of the spring depends linearly on the mass. An since the gravitational force depends linearly on mass it becomes clear that there is a linear dependence between force and deformation. This is can be modeled using a single constant called the spring constant $k$. We can write:

$F=k\cdot x$

For small deformations this stays true for complex systems, like welded frames made out of steel, or railway bridges. Since you can look at the deformation $\mathbf{x}$ at many points in the system and since there are many forces $\mathbf{F}$acting on such a system on constant is not enough to model the system correctly, you rather need a matrix of constants, this is called the stiffness matrix $\mathbf{A}$. The linear relation in this case can be written as:

$\mathbf{F}=\mathbf{A}\mathbf{x}$

An this text is basically about finding $\mathbf{A}$

Introduction

The equations of statics often allow to find all forces in a static system. Sometimes it is not possible to find a unique solution and additional assumptions are needed in order to find one. Often one assumes that the internal forces are zero if all external loads are zero. The finite element method is one way of finding a solution under this assumption. It is further assumes that the internal forces depend linearly on the external loads and that the deformations of the system depend linearly on the internal forces.

Calculation of stiffness matrices

In this text I am going to calculate some basic FEM stiffness matrices for members. The results of our calculation can be checked with: mvm.pdf or with frame-mth.pdf

According to Hooks Law. The relationship between stress and strain. Is given by:

$\sigma = E \epsilon$

Since stress is defined as force per area we go:

$F= A \sigma$

Furthermore we got the definition of strain:

$\epsilon=\frac{\triangle L}{L}$

From these equations we conclude:

$F=\frac{AE}{L}\triangle L$

Now we look at a member parallel to the x axis, with two end nodes called 1 and 2. Moving any of the two nodes for an infinitesimal amount of distance in x direction, we will cause an infinitesimal amount of force in x direction. But moving any node in y direction will not cause any change of force since the length of the member will not change if consider only the first order approximation. So we got the following relation between displacement and force.

$\left( \begin{matrix} F_{x1} \\ F_{y1} \\ F_{x2} \\ F_{y2} \end{matrix} \right) = \frac{AE}{L} \left( \begin{matrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \\ \end{matrix} \right) \left( \begin{matrix} \triangle x_1 \\ \triangle y_1 \\ \triangle x_2 \\ \triangle y_2 \end{matrix} \right)$

We will use the symbol $\mathbf{A}$ for the above matrix, $\mathbf{F}$ for the above vector of forces and $\mathbf{x}$ for the above vector of displacements. If we want to calculated the forces or a beam that has got an arbitrary position in space, we can rotate the system to the parallel and then rotate the results back to the original position. Thus we got:

$\mathbf{F}= \mathbf{R} \mathbf{A} \mathbf{R}^{-1} \mathbf{x}$

Where $\mathbf{R}$ is the matrix describing the rotation in our case we got:

$\mathbf{R}= \left( \begin{matrix} c & -s & 0 & 0 \\ s & c & 0 & 0 \\ 0 & 0 & c & -s \\ 0 & 0 & s & c \\ \end{matrix} \right)$ Where $s=\sin(\alpha)$ and $c=\cos(\alpha)$ and $\alpha$ is the angle of rotation. The inverse matrix is the same with the signs of the sines swapped. It is easy to calculate $\mathbf{R} \mathbf{A} \mathbf{R}^{-1}$ by blocks. So we calculate:

$\left( \begin{matrix} c & -s \\ s & c \end{matrix} \right) \left( \begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix} \right) \left( \begin{matrix} c & s \\ -s & c \end{matrix} \right) = \left( \begin{matrix} c & 0 \\ s & 0 \end{matrix} \right) \left( \begin{matrix} c & s \\ -s & c \end{matrix} \right) = \left( \begin{matrix} c^2 & sc \\ sc & s^2 \end{matrix} \right)$

So the overall result is:

$\mathbf{R} \mathbf{A} \mathbf{R}^{-1} = \frac{AE}{L} \left( \begin{matrix} c^2 & sc & -c^2 & -sc \\ sc & s^2 & -sc & -s^2 \\ -c^2 & -sc & c^2 & sc \\ -sc & -s^2 & sc & c^2 \\ \end{matrix} \right)$

Bending

We start with the Euler-Bernoulli beam equation

$\frac{\partial^2}{\partial x^2}u(x)=\frac{M(x)}{EI}$

We assume $u(0)=0$ and $\frac{\partial u}{\partial x}(0)=0$ and start with a beam that has clamped on both sided. We thus got:

$\begin{matrix} M(x)&=&M-Fx \\ EI \frac{\partial u}{\partial x}(x)&=& Mx-\frac{1}{2} F x^2 \\ EI u(x)&=&\frac{1}{2}Mx^2-\frac{1}{6} F x^3 \end{matrix}$

From the static equilibrium we can infer assuming a beam of length $L$ and a momentum $M_2$ and a force $F_2$ at the other clamping:

$\begin{matrix} M+M_2-F\cdot L&=&0 \\ F+F_2&=&0 \end{matrix}$

To find a unique solution we have to fix two more degrees of freedom. Assuming $\frac{\partial u}{\partial x}(L)=0$ and $u(L)=y$ we find:

$\begin{matrix} \frac{1}{2}ML^2-\frac{1}{6} F L^3 &=&EIy \\ ML-\frac{1}{2} F L^2 &=&0 \end{matrix}$

Solving this for $F$ and $M$ we find:

$\begin{matrix} F&=&\frac{12EI}{L^3}y \\ M&=&\overset{\text{ }}{\frac{6EI}{L^2}y} \end{matrix}$

Thus we see that these two quantities are linear in the displacement $y$. Which means that we just calculated two coefficients of the stiffness matrix. We can find two more by using the equations of static equilibrium.

$\begin{matrix} M_2&=& -M+FL=\frac{6EI}{L^2}y\\ F_2&=&-F=\overset{\text{ }}{-\frac{12EI}{L^3}y} \end{matrix}$

The next case we have to consider is $u(L)=0$ and $\frac{\partial u}{\partial x}(L)=-b$ thus we got:

$\begin{matrix} \frac{1}{2}ML^2-\frac{1}{6} F L^3 &=&0 \\ ML-\frac{1}{2} F L^2 &=&-b \end{matrix}$

Solving this for $F$ and $M$ we find:

$\begin{matrix} F&=&\frac{6EI}{L^2}b \\ M&=&\overset{\text{ }}{\frac{2EI}{L}b} \end{matrix}$

The choice of the sign of b is not perfectly clear to me. Well graphically it is clear that the slope has to be negative in order to cause a positive force. Again we can use the static equations to infer two more equations.

$\begin{matrix} M_2&=& -M+FL=\frac{4EI}{L^2}b\\ F_2&=&-F=\overset{\text{ }}{-\frac{6EI}{L^3}b} \end{matrix}$