Sunday, October 22, 2023

Simulation #3: Bending of 2D beam under static force (part 1 - theory of mechanics for beam)


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 2D beam with length L and applied force $\vec{P}$. What is typical for beam - external force is applied perpendicularly to beam axis $x$. Beam here is in the state before simulation (before it is bent).

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:
Fig 1.2 2D beam from fig. 1.1, state after 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.
We know that first case is not what relevant to our beam as there is force $\vec{P}$ acting on it. Regarding second bullet, it is met when:
  1. sum of projections of forces onto x-axis is zero: $$\sum F_x=0\qquad\qquad (1)$$
  2. sum of projections of forces onto y-axis is zero: $$\sum F_y=0\qquad\qquad (2)$$
  3. sum of moments at any point is zero: $$\sum M_i=0\qquad\qquad (3)$$
Above equations are called external equilibrium equations (for 2D structure). Our 2D beam case, is the system, in which we have 3 degrees of freedom: movement along x-axis, movement along y-axis and rotation. They are related respectively to equations: 1st (moving along x-axis), 2nd (moving along y-axis) and 3rd (rotation).


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).
    Fig 3.1 Displacement of point A as effect of acting force.
    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.2 Deflection function for bending beam.
    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.
  • 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}}$$
Stress expresses the internal forces that neighbouring particles of a continuous material exert on each other, while strain is the measure of the relative deformation of the material. But what is relation between them?


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):
Fig 4.1 Stress–strain curve for steel (picture from Wikipedia)
We would like to apply enough small force that beam can behave as linear-elastic solid (in range from 0 up to point 2 in above fig.). In that range formula (4.1) is applicable.


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).
Fig 5.1 External and internal forces acting on bending beam.
When we cut beam into two parts at x – we have left and right part that reacts on each other by inner forces in the crosssection. Those forces are equal but acts in opposite direction. Therefore we can consider only left part of the beam – fig 5.1 b). In its crosssection we can define small section called $\Delta h$, then we can say that $\Delta P_w$ is total force that acts on that section as reaction of external forces acting on beam. $\Delta P_w$ is internal force. It is related to stress value. Let's figure out this relation.

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.
Fig 6.1 Moments acting on bending beam.
Analogically as previously, lets cut the beam into two parts left and right (figure 6.1 b)). We can see, that in its crossection, there is another moment called bending moment Mb. In the right part there is moment with same lenght but acting in oposite direction to balance the one from left part in given crossection. There is relation between $M_b$ and moments of external forces like for example moment with relation to point $M_A$: $$M_b+\sum\limits_{i}M_{li}=0,$$ $$M_b+\sum\limits_{i}M_{ri}=0,$$ where $M_{li}$ are all moments of external forces acting on left part of the beam, $M_{ri}$ are all moments of external forces acting on right part of the beam. Therefore: $$M_b=-\sum\limits_{i}M_{li},$$ $$M_b=-\sum\limits_{i}M_{ri}.$$ Notice that it is sum of vectors.


7. Sign convention.


Below figure describes sign convention that are valid in Statics for quantities like: $M_b$, T, N.
Fig 7.1 Sign convention for normal force, shear force and bending moment.
Sign of moment which is defined with respect to a fixed reference point, for example $M_A$ (with respect to point A) is determined differently than sign of bending moment $M_b$:

Moment M is positive if it is rotating clockwise: Moment M is negative if it is rotating counterclockwise:
Fig 7.2 a) Sign of moment in relation to point is positive when it is rotating clockwise.
Fig 7.2 b) Sign of moment in relation to point is negative when 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.
Fig 8.1 Forces, moments and its reactions acting on bending beam.
From equilibrium equations we have:

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:
Fig 8.2 Bending moment and sheer force plots for bending beam.
As we were able to calculate reaction forces and moments using equilibrium equations, therefore we can say that beam is statically determined.


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


Fig 10.1 Beam fragment for normal stress formula deriviation.
In above figure it is presented small part of bending beam. We can see that upper edge of that beam fragment ($y=h/2$) gets expanded, whereas bottom edge ($y=-h/2$) gets compressed. Because one edge extended and opposite one compressed, therefore there should be a fiber between them that does not change its length when beam is bending. With blue line I marked, so called neutral axis, which is that particular fiber. In case of beam with regular rectangular shape (our beam), neutral axis is at the center of crossection ($y=0$). In further part of this article, we will see why it is like that. We assume that face of the cut $|BD|$ remains plane (pure bending), and it remains perpendicular to the deformed axis of the beam (neutral axis).

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.
Fig 10.2 Beam fragment, bent due to external force with marked normal stress values on the right crossection.
Because part of the beam is in static equilibrium, we can use static equilibrium equations: $$\sum F=0,$$ because normal stress $\sigma$ is related to distributed force over the crossection, therefore: $$\int\limits_{-h/2}^{h/2}\sigma dy=0,$$ taking out constants: $$\frac{E}{R}\int\limits_{-h/2}^{h/2}y dy=0.$$ It can be seen, that existence of blue line (neutral axis) which we assume before is now proofed because it is determined by integral from above equation of static equilibrium.

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:
Fig 11.1 Shear force and bending moment acting on the beam cut with corresponting shear and normal stresses.
We have seen that the normal stresses due to bending moment M are linearly distributed over the cross section, with maximum magnitudes of normal stress occuring on the outer fibers of the beam and with zero normal stress at the neutral axis (the neutral axis passing through the centroid of the cross section). In our case, shear force T acts along with the bending moment $M_b$ and a component of shear stress will exist. Shear strains correspond to a change in angle of the stress element. This angle change is somewhat in contradiction with the pure bending assumption of the cross section remaining perpendicular to the deformed beam axis. For our normal stress formula, we will assume that the shear strain effects will be slight and that, even in the presence of shear stress, the distribution of flexural stress on a given cross section is unaffected by the deformation due to shear and flexure formula equation is still valid for computing the normal stresses on the cross section.

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).
Fig 11.2 Beam cut element with marked all types of stresses acting on it.
Consider the aribitrarily-loaded beam shown below:
Fig 11.3 Aribitrarily-loaded beam with selected cut.
As can be seen, we isolate a section of the beam between locations $x$ and $x + \Delta x$. Zoom of this part of beam is presented in below figure:
Fig 11.4 Zoom of selected cut from fig. 11.3 with marked forces, moments and stresses.
The the resultant of shear forces and bending moments acting on the beam section is shown in 11.4 a). The resultant bending moments $M(x)$ and $M(x + \Delta x)$ produce normal stresses of $\sigma(x)$ and $\sigma(x + \Delta x)$ on the left and right faces of the beam section, respectively (fig. 11.4 b)).

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$.
After substitution of constant values and transformations, this yields to final formula of displacement function: $$v(x)=\frac{6Px^3}{Eh^3}\left(L-\frac{h}{3}\right).$$

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):
Fig 13.1 Sheer force plot for bending beam
Fig 13.2 Bending moment plot for bending beam
Fig 13.3 Displacement function plot for bending beam
Fig 13.4 Normal stress plot for bending beam
Fig 13.4 Sheer stress plot for bending beam

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.