I. Introduction
Previous simulation #2 was related to phenomenon of ball in gravitational field about which we were teached at school. Formulas describing the ball movement are known and simple (they do not contain derivative equations in their basic form). Any quantity like x(t), y(t), V(t) etc. can be expressed by its own formula, which mathematicaly is same simple.
This blog is about numerical methods. Numerical methods are used to solve formulas that does not have analytic solution or they are very difficult to find. For instance formulas describing physical phenomenon occurring on the top of a body with irregular geometry. If it is possible to find analytically, it would take much time.
This article is related to simulation from area of structural engineering. We will consider the case of 2 dimentional beam under static force. Formulas describing it are complex and analytical solution may be difficult to find. Luckily, numerical methods are our helpful friend. We will use Finite Element Method (FEM) to find approximate solution - values of different parameters that describe the beam.
II. 2D beam - theory
This article is theoretical and provide required understanding about theory of the beam, from the field of mechanical engineering. The next part will be introduction to FEM method. Theory will be limited to particular case that we deal with as topic is too large to provide lecture for any beam case in one article.
1. What is beam?
First question you may ask is: what is beam? According to Wikipedia:
a beam is a structural element that primarily resists loads applied laterally to the beam's axis (an element designed to carry primarily axial load would be a strut or column).What is specific for beam and different from other structural elements is that it has one dimention much longer than than other two. For 2 dimentional case, its lenght is much longer than its high. Below picture presents beam that will be modeled. That case has its own name - cantilever beam. It is a beam fixed on one left edge (for example bolted to the wall) and the right edge is free. Static force is acting on the right edge downwards.
Fig 1.1 presents the initial state (before simulation). After simulation beam bends under acting force $\vec{P}$ or in language of engineering - displacement occurs. Fig. 1.2 presents that state of simulation:
2. Static equilibrium
To start beam analysis we need to introduce fundamental term of static equilibrium. Statics (or Engineering Statics) is part of Engineering Mechanics that deals with rigid bodies (such as beams) that are in equilibrium state. Equilibrium means that they are not moving. According to Newton's Law, body is not moving when:
- there are no forces acting on body,
- all forces acting on the body are balanced.
- sum of projections of forces onto x-axis is zero: $$\sum F_x=0\qquad\qquad (1)$$
- sum of projections of forces onto y-axis is zero: $$\sum F_y=0\qquad\qquad (2)$$
- sum of moments at any point is zero: $$\sum M_i=0\qquad\qquad (3)$$
3. Quantities describing beam
There are multiple quantities describing beam. Three especially important are:
- displacement - (in Polish: przemieszczenie) movement of a point that belongs to a beam due to its bend caused by external force. Below figure presents displacement of point A (point is marked as A' after displacement).
It can be seen that force $\vec{P}$ causes mainly displacement of y component of the point denotes by v. Displacement of x component denoted by u is relatively much smaller. It makes sense as force $\vec{P}$ acts perpendicular to y axis. We can introduce function v(x) which represents displacement function in y-axis direction (u deflection is much smaller, therefore we can skip it). Below picture presents v(x), we can call it deflection/displacement curve:
Fig 3.1 Displacement of point A as effect of acting force. Based on above figure true is that: $$v(x)=\frac{d\theta}{dx}$$ At the end of this article we will derive formula for v(x) for our beam.
Fig 3.2 Deflection function for bending beam. - stress - (in Polish: naprężenie) denoted as $\vec{p}$ it is vector expressing the intensity of internal force. It can be describe as force per unit area or intensity of forces distributed over given section. Internal forces in the beam are effect of external force (force $\vec{P}$ in our case). Stress vector $\vec{p}$ is sum of two vectors - normal stress vector ($\vec{\sigma}$) and sheer stress vector ($\vec{\tau}$): $$\vec{p}=\vec{\sigma}+\vec{\tau}$$
- strain - (in Polish: odkształcenie) denoted as ($\epsilon$) it is the change of displacement per unit length: $$\epsilon = \frac{\text{final_length} – \text{initial_length}}{\text{initial_length}}$$
4. Stress strain relation
Beam behaviour under the force to some extent is similar to spring behaviour under the force. Commonly known is law, which says that spring deflection from equilibrium (state without force) is proportional to force applied. It is called Hooke's law. For the beam it is analogical equation, which is applicable for small deflections: $$\sigma = E\cdot \epsilon,\qquad\qquad (4.1)$$ where:
E – is Young modulus, which is measured by experience. It is constant value specific for material that beam is made from.
Below figure presents relation between strain and stress (stress – strain curve):
5. What forces are acting on the beam?
We mentioned in previous section that there are 2 types of forces - external and internal. Below picture (5.1) presents types of forces that are acting on the beam. In part a) there is external static force $P$ that we already know. Another external force is force $R_A$, which is reaction force caused by fixed support (wall).
When we limit size of $\Delta h$ to zero: $$\lim_{\Delta h \to 0} \frac {\vec{\Delta P_w}}{\Delta h} = \vec{p},$$ where p is called stress vector. It is intensity of internal force at specific point. In 2D case crossection is just line therefore p depends only on y: p(y). Stress vector $\vec{p}$ is sum of two vectors: normal stress vector ($\vec{\sigma}$) and sheer stress vector ($\vec{\tau}$) as we state before. Those 2 types of stress corresponds respectively to 2 types of internal forces: normal force and sheer force.
We can calculate average vector of p for the whole crosssection when we calculate it in its centroid (center of mass/center of an area, in our 2D case it is center of crossection y=0). Total $p$ is marked as $P_w$ and presented at fig. 5.1 c). When we cast it on x and y axis, we get respectively total normal force N and total sheer force T for that crossection. Their values are: $$N=\int \limits_{h} \sigma dy,$$ $$T=\int \limits_{h} \tau dy,$$ where integral is calculated on crossection which in this case is (2D) line segment with lenght h. There is relation between $P_w$ and external forces acting on beam: $$P_w+\sum \limits_{i}P_{li}=0,$$ $$P_w+\sum \limits_{i}P_{ri}=0,$$ where $P_{li}$ are vectors of all external forces acting on left part of the beam, $P_{ri}$ are all external forces vectors acting on right part of the beam. It make sense that, we can calculate internal force $P_w$ considering either left or right part of the beam, because both of them stays in static equilibrium.
Above equation can be rewrite as: $$P_w=-\sum \limits_{i}P_{li},$$ $$P_w=-\sum \limits_{i}P_{ri}.$$
6. What moments are acting on the beam?
Now it is time to consider moment of forces that are acting on our beam. On the left we have wall into which beam is attached. In the language of Statics it is fixed support (cantilever). This kind of support produce both reaction force and reaction moment $M_A$ as well – figure 6.1 a). How to understand reaction moment? Imagine you are holding a long board by one end and lifting the other end. Where you hold it, you have to rotate your wrist to keep the board balanced. The same happens with the moment in a cantilever. Why does this happen? This is due to the balance of forces and moments. Introducing a reaction moment helps the beam maintain balance by counteracting the forces and moments acting on it. Therefore $M_A$ has orientation that makes it counteracting of bending of the beam.
7. Sign convention.
Below figure describes sign convention that are valid in Statics for quantities like: $M_b$, T, N.
| Moment M is positive if it is rotating clockwise: | Moment M is negative if it is rotating counterclockwise: |
|
|
8. Finding reaction force $R_A$, reaction moment $M_A$, sheer force $T$ and bending moment $M_b$
If we ignore wall on the left, but take into account moment and reaction that it is causing, we can try to calculate forces $R_A$, $M_A$ using equilibrium equations.
ad. 1) $\sum F_x=0 \quad \implies 0=0$ (no external forces along x-array)
ad. 2) $\sum F_y=0 \quad \implies -P+R_A=0 \implies R_A=P$
ad. 3) lets consider point x=0; y=0. Moment in it should be zero: $$\sum M_{(0;0)}=0 \quad \implies P \cdot L–M_A=0 \implies M_A=P\cdot L$$ $M_A$ and $R_A$ values we can use now to calculate $M_b(x)$ and $T(x)$. They can be determined when we split beam into 2 parts – left and right and considerint moments and forces only from left part. As mention before, we are allowed to consider left part only due to fact that it is in static equilibrium. It is possible to analyze right part of course and it should give the same results. For the crossection of cut (at x in above fig.), bending moment is equal to sum of all moments from left part of the beam in respect to center point of that crossection ($y=0$): $$M_b(x)=R_A\cdot x-M_A = P\cdot x - P\cdot L= P(x-L)$$ According to sign convetion $R_A$ force makes beam bent upwards so its moment's ($R_A\cdot x$) sign is positive, moment $M_A$ whereas makes beam bent downwards, so its sign is negative.
Sheer force T at x is equal to sum of forces acting along y-axis on the left from crossection: $$T(x) = R_A = P$$ According to sign convention figure, $R_A$ force is positive, because it is on the left from crossection and it is pointing upwards.
$T(x)$ can be confirmed using relation $T(x) = \frac{dM_b(x)}{dx}$ (Schwedler's Theorem). It is equal to P, therefore value is right.
Now we can draw graphs for both cases:
9. Beam bending classification
When beam bending happened due to moment only (when sheer force T is zero) we call it pure bending. Then crossection shape remains flat during bending. When sheer force exists also then bending is called non-uniform. Sheer force makes bending more complicated, when it is strong comparing to beam strenght against lateral force and crossection shape may be distorted during bending instead of remaining flat. We can see, that in our case bending moment $M_b$ is non-zero for all beam lenght (except from point $x=L$) and also sheer force T is non-zero. It is very difficult to calculate in such case analitically beam parameters like stresses, displacement function. We are able to do this numerically using finite element method. However for now, let's perform analytical calculations applying common simplification, used when lenght of the beam is much greater than its width, and then impact of sheer force can be skipped. We will use this simplification when calculating formula for normal stress. So we will assume that bending is pure: $M_b\neq 0$, $T=0$.
In such model, we can imagine the beam to be made up of longitudinal fibers parallel to the longitudinal axis of the beam. Each fiber has crossection $\Delta h$ and lenght L. Fibers does not acts on each other along lateral direction (y-axis).
10. Normal stress formula deriviation
Above beam fragment we will use to derive formula for bending stress $\sigma$. Firstly, let's write what is the length of upper edge of beam fragment: $$\epsilon = \frac{\text{final_length} – \text{initial_length}}{\text{initial_length}}$$ We assume that blue line (neutral axis) does not change its length so it has initial length after bend. Its length is: $$\text{inital_length} = |EF| = \theta \cdot R$$ $$\text{final_length} = |AB| = \theta \cdot (R+ \frac{h}{2})$$ Generally, final_lenght for any fiber is: $$\text{final_length} = \theta \cdot (R+y),$$ therefore strain formula is equal to: $$\epsilon = \frac{\theta\cdot (R+y)-\theta R}{\theta \cdot R}= \frac{y}{R}$$ Note: often in books, in above formula there is minus on RHS. It depends on how we set XY coordinate system.
We are talking about linear – elastic solid, therefore relation between stress and strain is: $$\sigma = E\cdot \epsilon$$ When we substitute $\epsilon$ to above formula, we get: $$\sigma = E \cdot \frac{y}{R},$$ we can see that the highest value of stress we get at the bottom and at the top of the beam, however sign in both cases are different. For extension we have positive, for compression negative. Above formula for sigma depends on radius R. Let's do more transformations and get rid of it. Below figure presents again our beam part with marked $\sigma$ at the right edge of the beam part.
For second equilibrium equation:
$\sum M =0$, therefore $$\int\limits_{-h/2}^{h/2}\sigma y dy =M.$$ When we substitute $\sigma$ value: $$\frac{E}{R}\int\limits_{-h/2}^{h/2}y^2 dy=M.$$ Integral in above equation has special name – it is Moment of Intertia I. Therefore: $$M=\frac{EI}{R}.$$ We derived 3 relations so far: $$\epsilon=\frac{y}{R}; \quad \sigma=\frac{Ey}{R}; \quad M=\frac{EI}{R}.$$ We can combine 2nd and 3rd relation to get rid of R parameter: $$\sigma=\frac{My}{I},$$ the above is formula is important relation in Statics and it is called flexure formula. It says, when we know the moment that is applied to the beam and its moment of inertia for particular crossection, then we can calculate maximum stress in it.
For 2D case, moment of inertia I we can called $I_y$ and calculate as: $$I_y=\int\limits_h y^2 dy=\left[\frac{y^3}{3}\right]^{h/2}_{-h/2}=\frac{h^3}{12}.$$ Now formula for normal stress is: $$\sigma=\frac{12My}{h^3},$$ we can finally substitute known value of M calculated before: $$\sigma=\frac{12(Px-PL)y}{h^3}=\frac{12(x-L)Py}{h^3}.$$
11. Sheer stress formula deriviation
Now we are going to derive formula for shear moment $\tau$. Below picture shows both shear force and bending moment acting on beam cut:
Suppose we consider a stress element on the side of a beam with a non-zero shear force resultant on the face of the cut (Fig. 11.2). Our goal here is to determine the transverse shear stress component $\tau_{xy}$ that corresponds to the shear force resultant T. Note, however, that since $\tau_{yx} = \tau_{xy}$ (based on the rule that sheer stress in mutually perpendicular planes are equal*), the transverse shear stress component $\tau_{xy}$ is the same as the longitudinal shear stress component $\tau_{yx}$. Stated in different words, we can determine the transverse shear stress by calculating the longitudinal shear stress.
(*) rule is known in Statics, it is easy to proof – try to calculate moment in any point of the element and equal it to zero (according to equlibrium equation).
Normal stress on beam section value is: $$\sigma(x,y)=\frac{-M(x)y}{I}, \\ \sigma(x+\Delta x, y) = \frac{-M(x+\Delta x)y}{I}.$$ Total sheer force resultants on beam section (marked in figure 11.4 a)) is: $$T(x)=\frac{dM}{dx}.$$ Suppose we further isolate a slice of this beam section, found below a given value of y – it is marked as $A'$ in fig. 11.4 c) also present in figure 11.4 d). Sheer stress $\tau$ is related to sheer force. It is acting on both left and right crosssections of the beam part and as well in longitudinal direction (fig. 11.4 d)).
We can define forces at $x$ and $x + \Delta$, acting on slice A' that are corresponding to normal stress $\sigma$: $$F(x)=\int \sigma(x,y)dh', \\ F(x+\Delta x)=\int \sigma (x+\Delta x,y)dh'.$$ Thre is one more force factor along x-axis which corresponds to longitudinal sheer stress. Longitudinal sheer stress acts on the upper surface of the slice at y. Corresponding force to it is marked as $\Delta H$ in fig. 11.4 d).
Because $A'$ element is in equilibrium state, total force along x-axis should be zero: $$\sum F_x=F(x)-F(x+\Delta x)+\Delta H=0 \implies \Delta H=F(x+\Delta x)-F(x).$$ The shear stress corresponding to this resultant shear force is found from the usual definition of stress in terms of the force resultant as: $$\tau=\lim_{\Delta x \to 0} \frac{\Delta H}{\Delta x} =\lim_{\Delta x \to 0} \frac{F(x+\Delta x)-F(x)}{\Delta x} =\frac{dF}{dx}.$$ We can calculate F(x) using below relation: $$F(x)=\int \sigma (x,y) dh'=\frac{M(x)}{I}\int y dh'=\frac{M(x)}{I}\int \limits_{-h/2}^{-y}ydy=\frac{M(x)}{I}\left[ \frac{y^2}{2} \right]_{-h/2}^{-y}= \\ =\frac{M(x)}{I}\frac{1}{2}\left(\frac{h}{2}-y \right)\left(\frac{h}{2}+y\right).$$ Combining above equations gives: $$\tau=\frac{\left(\frac{h}{2}-y \right)\left(\frac{h}{2}+y\right)}{2I}\frac{dM(x)}{dx}=\frac{\left(\frac{h}{2}-y \right)\left(\frac{h}{2}+y\right)}{2I}T,$$ where we used relation (Schwedler's Theorem) $T = dM(x)/dx$.
Moment of inertia $I = h^3/12$ for our beam, therefore final formula for sheer stress is: $$\tau= \frac{\left(\frac{h}{2}-y \right)\left(\frac{h}{2}+y\right)12}{2h^3}T=\frac{\left(\frac{h^2}{4}-y^2\right)6}{h^3}T.$$
12. Displacement function formula
Now it is time to find value of displacement function v(x). Below function, which can be find in books for Statics, is differential equation that describes relationship between displacement function of beam along y-axis (marked as v(x)) with other parameters that we introduced before: $$EI\frac{d^2 v(x)}{dx^2}=-M_b.$$ Substituting value of $M_b$, we have: $$EI\frac{d^2 v(x)}{dx^2}=PL-Px.$$ To find v(x), we need to integrate equation twice.
1st integration gives: $$EI\frac{dv(x)}{dx}=EI\theta=C+PLx-\frac{1}{2}Px^2. \qquad\qquad (12.1)$$ 2nd integration: $$EIv(x)=D+Cx+\frac{1}{2}PLx^2-\frac{1}{6}Px^3. \qquad\qquad (12.2)$$ From border condition we can calculate constant values C and D:
- (from eq. 12.1) for $x=0$, $\theta=0 \implies C=0$,
- (from eq. 12.2) for $x=0$, $v(0)=0 \implies D=0$.
13. Example
We will solve our beam case for below values:
L=6m
E=20 000 MPa
h=800mm
P=150 kN
Using SymPy Python's library I programmed equations derived in previous chapters for sheer force, bending moment, deflection function, normal stress and sheer stress. Substituting above values, I created graph for each parameter (using Matplotlib library):
III. Summary
This article was introduction for Mechanics of the beam. In the next part, I will describe basics for FEM method based on our beam case - so further, we can proceed with the actual simulation.
