Beam Structures: Plane Beam Approximations
Many of the large standing towers, bridges, and buildings owe their existence to some simplified versions of the equilibrium equations applied to slender, linelike structures. The first beam approximation is represented in the Euler Bernoulli’s beam theory, which originated sometime in the eighteenth century. This approximation has been so successful that it is still the basic approximation used in modern structural analysis taught in current structural analysis courses. In this chapter, three different beam approximations are presented. The first is the the Euler Bernoulli’s beam approximation which relies on assuming plane sections perpendicular to the neutral axis of the beam remain plane and perpendicular to the neutral axis after deformation. The second beam approximation is the Timoshenko’s beam approximation which assumes that the cross sections perpendicular to the neutral axis before deformation stay plane, but not necessarily perpendicular to the neutral axis after deformation. These beams are termed: “shear flexible” and because they allow more deformation, the model predicts slightly more deformation than the Euler Bernoulli beam model. The third beam approximation is the beams under axial loading. All these models are presented under the umbrella of small deformations and for isotropic linear elastic material models.
Euler Bernoulli Beam:
There are three basic assumptions for an Euler Bernoulli beam that will be used to derive the equations. These are:
 Plane sections perpendicular to the neutral axis before deformation stay plane and perpendicular to the neutral axis after deformation (Figure 1).
 The deformations are small.
 The beam is linear elastic isotropic and Poisson’s ratio effects are ignored.
The beam is assumed to lie such that its long axis is aligned with . The deformation is controlled by a function which describes the deformation of the neutral axis (Figure 1). The first assumption allows the calculation of the small strain matrix at every point. The coordinates of an arbitrary point before deformation are given by:
Ignoring any effect due to Poisson’s ratio, the coordinates of the point after deformation (Figure 1) are given by:
The second assumption (small deformation) allows the approximations and . Therefore, the displacement function after incorporating the small deformations assumption is given by:
We can drop the superscript since this applies to any arbitrary point. We also note that . Therefore, the displacement function has the form:
we remind the reader that is a function of only and so is its derivative. Therefore, The gradient of the displacement matrix is given by:
Therefore, the small strain matrix has the form:
In essence, the deformation assumptions result in all of the strain components being zero except for . Furthermore, we will look carefully at the variation of on a given cross section. For a given cross section of the beam we have is equal to constant. Therefore, the quantity is constant on the cross section and therefore is linear in the vertical direction (Figure 2). The value of is equal to zero at the neutral axis . The positive convention for and for lead to a positive strain below the neutral axis with the maximum positive value at the bottom fibre of the cross section. The positive convention also leads to a negative strain above the neutral axis with the maximum negative value at the top fibre of the cross section. Using the linear elastic constitutive relationship for the beam material and ignoring Poisson’s ratio lead to the same distribution for (Figure 2). Therefore, the normal stress component distribution on the cross section is given by:
(1)
Internal Forces:
The deformation assumption in beams allows considering only the deformation of the neutral axis. In particular, the deformation of the beam is controlled by the unknown function which is the vertical deformation of the neutral axis. We can also look at the forces acting on the cross section that is perpendicular to the neutral axis. These forces are termed: “Internal Forces” and can be obtained by integrating the various stress components acting on the cross section under consideration (Figure 3). For a beam in plane, only two internal forces are considered, the bending moment , and the shearing force . The positive convention for these internal forces is shown in Figure 3. can be obtained by integrating the stress component on the cross section:
The negative sign is introduced to keep the positive convention for causing compression on the top fibre of the beam. By substituting for using Equation 1, the moment can be calculated as:
(2)
where is the moment of inertia of the beam calculated according to the formula:
Equations 1 and 2 can be used to find the relationship between the internal bending moment and the stress component on the cross section:
(3)
Equilibrium Equations:
It is possible to derive the equilibrium equations for the Euler Bernoulli beam by integrating the general equilibrium equations of a continuum. However, we will follow the classical method as described in many structural analysis books following the sign convention shown in Figure 3. First, two planes that are very close to each other (separated by a small distance and perpendicular to the axis are considered. Then, by assuming that and are smooth functions of , then the first approximation in a Taylor’s series can be utilized to describe the changes in the internal forces between those two planes (Figure 4). The loading on the beam is assumed to be transverse in the direction of upwards and given by with units of load per beam length. Assuming the sum of forces in the vertical direction (direction of ) is equal to zero yields:
Similarly, assuming the sum of moments around is equal to zero yields:
Since , and are only functions of , the partial derivatives can be replaced with total derivatives and the equilibrium equations become:
(4)
We can substitute for using Equation 2 so the equilibrium equation in terms of the displacement function and the loading becomes:
(5)
The solution for has the following form:
Four boundary conditions are needed to fully solve for and its derivatives. These boundary conditions can be given in terms of the displacement , the rotation , the bending moment or the shearing force . Note that if the beam can have a polynomial shape up to the third degree. For this reason, this formulation is sometimes termed: The cubic formulation.
Shear Stress in Euler Bernoulli Beam:
The small strain matrix obtained above for the Euler Bernoulli beam shows that the shear strains are equal to zero. However, this is an approximation that simplifies the beam model. While shear strains are directly proportional to shear stresses in linear elastic materials, for the Euler Bernoulli beam, shear stresses arise from the existence of the shearing force even though the shear strains are equal to zero. The shearing stresses associated with can be calculated by first assuming that the shear stress is constant across the beam width (along the direction). Then, by considering a small slice in the beam with width of the shear stress at the point represented by the black circle in Figure 5 can be calculated by considering the equilibrium of the shaded section shown in the figure. There are three horizontal forces acting on the shaded section. On the left we have caused by integrating on the left cross section of the shaded area. On the right we have caused by integrating on the right cross section of the shaded area. At the bottom we have caused by integrating on the bottom of the shaded area. Using Equation 3 that relates with , the three forces can be calculated as follows:
where is the width of the cross section (the width in the direction of ). Equilibrium dictates:
Therefore,
We will use the equality: . Therefore:
Since , , and are not functions of or , they can be brought outside the integral. Denote to be the moment of the area defined by:
Therefore:
Since is equal to , the shear stress on the cross section at the point in question is equal to:
(6)
A similar proof can be found in the Mechanics ebook here.
Discrepancies of the Beam Theory:
The Euler Bernoulli beam model is a simple approximation. While it is very successful because of it simplicity and practicality, it has a few discrepancies. One discrepancy is that the shear stress is not associated with any shear strains. The second discrepancy is that the calculated stresses, while they satisfy equilibrium in a holistic sense, they don’t necessarily satisfy the equilibrium equation of the continuum at a local sense. In particular, the equation of horizontal equilibrium of the continuum is:
Notice that is assumed to be zero and that we are using the undeformed coordinates since small deformations are assumed. Also, no body forces are assumed to be acting on the beam. By setting and substituting for and using Equations 3 and 6 respectively, we get the following:
Assuming is constant along and since and are functions only of we get:
Therefore:
The left hand side is only equal zero if is constant (i.e., rectangular cross section). Otherwise, the left hand side calculation for any cross section other than a rectangle will give nonzero values, i.e., equilibrium is not necessarily satisfied!
Timoshenko Beam:
The Timoshenko beam formulation is intentionally derived to better describe beams whose shear deformations cannot be ignored. Short beams are a prime example for such beams, and thus, the Timoshenko beam approximation is better suited to describe their behaviour. The basic physical assumptions behind the Timoshenko beam are similar to those described for the Euler Benroulli beam, except that shear deformations are allowed. The following are the three basic assumptions behind the Timoshenko beam theory. (Compare with those described above for the Euler Bernoulli beam)
 Plane sections perpendicular to the neutral axis before deformations remain plane, but not necessarily perpendicular to the neutral axis after deformation (Figure 6).
 The deformations are small.
 The beam is linear elastic isotropic and Poisson’s ratio effects are ignored.
Similar to the Euler Bernoulli beam, the Timoshenko beam is assumed to lie such that its long axis is aligned with . The deformation is controlled by two functions: which describes the deformation of the neutral axis (Figure 6) and which describes the shear rotation of the cross section. In essence, the total rotation of the cross section is equal to . The coordinates of an arbitrary point before deformation are given by:
Ignoring any effect due to Poisson’s ratio, the coordinates of the point after deformation (Figure 6) are given by:
Setting , the small deformations assumption allows the approximations and . Therefore, the displacement function after incorporating the small deformations assumption is given by:
We can drop the superscript since this applies to any arbitrary point. We also note that . Therefore, the displacement function has the form:
we remind the reader that and (and hence ) are functions of only. Therefore, The gradient of the displacement matrix is given by:
Therefore, the small strain matrix has the form:
In essence, the deformation assumptions result in all of the strain components being zero except for and . Furthermore, we will look carefully at the variation of on a given cross section. For a given cross section of the beam we have is equal to constant. Therefore, the quantity is constant on the cross section and therefore is linear in the vertical direction (Figure 2). The value of is equal to zero at the neutral axis . The positive convention for and for lead to a positive strain below the neutral axis with the maximum positive value at the bottom fibre of the cross section. The positive convention also leads to a negative strain above the neutral axis with the maximum negative value at the top fibre of the cross section. On the other hand is constant on the cross section. Using the linear elastic constitutive relationship for the beam material and ignoring Poisson’s ratio lead to the same distribution for (Figure 2). Therefore, the normal stress component distribution on the cross section is given by:
(7)
The internal bending moment on the cross section can be obtained using the integral:
(8)
The shear stress distribution on the cross section is given by:
where is the shear modulus and is a constant introduced to allow for averaging or smearing the shear stress over the full cross section. was introduced by Timoshenko in his book. For more information on you can view this reference. The shear force is the integral of this average over the cross section which yields:
In essence, the term gives the corresponding “effective” shear area such that the shear stress is equal to .
Equilibrium Equations:
The equilibrium Equation 4 developed for the Eulber Bernoulli beam also applies to the Timoshenko beam (Figure 4). By substituting for using Equation 8 in the equilibrium equations, we reach:
(9)
We also have the following differential equation in terms of the beam displacement :
(10)
The solution for and have the following forms:
The above two equations (9 and 10) can be solved if four boundary conditions for the cross section rotation , the moment , the shear , and/or the displacement are available. At least one boundary condition for has to be given in order to find the integration constant . Note that if then the Timoshenko beam solution would approach the Euler Bernoulli beam solution.
Beam under Axial Loading:
The Euler Bernoulli and the Timoshenko beam formulations described account only for the deformations due to bending, without considering any axial deformation in the neutral axis of the beam. It is possible to augment the Euler Bernoulli and the Timoshenko beams so that they also account for axial deformations of the neutral axis, but because of the assumption of small deformations, the resulting equations of axialloading deformations and bending deformations are uncoupled. i.e., the axial loading is only related to the axial deformation, and the lateral loading is only related to the lateral deformation. This uncoupling is only valid for small deformations. However, for large deformations, and normal force in an axially loaded beam can affect the lateral deformation as well. For example, phenomena such as buckling can only be modelled when the interaction between the axial loads and the lateral deformations are coupled.
Deformation Assumption:
Under an axial load, a beam is assumed to deform such that the horizontal displacement is a smooth function and that cross sections remain plane during the deformation (Figure 7). The effect of Poisson’s ratio is usually neglected. Thus, the position vector function is:
In this case, the displacement function is:
Since is a function of only , the only nonzero component of the strain is :
On any cross section, is constant and therefore the values of and its derivatives are constant. Therefore, is constant on any cross section perpendicular to .
Equilibrium Equation:
The beam is assumed to be under a distributed normal load of value of units of force/unit length. Figure 8 shows the applied distributed load along with the equilibrium of a beam slide of length . The equilibrium equation can be written in terms of the normal force which is the integral of the stress component over the cross sectional area:
Since is constant we have:
where is the cross sectional area. Equilibrium can be written by considering the three forces acting on the slice in Figure 8. The first force in acting on the right cross section and directed to the right. The second force is acting on the left cross section and directed to the left. The third is due to the distributed load . Therefore:
Equilibrium dictates:
Therefore:
(11)
Assuming a linear elastic beam, the differential equation can be written in terms of the displacement by converting the stresses into strains:
Therefore, Equation 11 can be written as follows:
(12)
Or equivalently, if is constant:
In the case when the cross sectional area is constant, the equilibrium equation take the simple form:
Visualization Tools:
Difference Between Euler Bernoulli and Timoshenko beams:
The following tool analyzes a cantilever beam of length 10 units and height of 1 unit under an upward unit load. The displaced shape of the beam is drawn for various values of and . The Euler Bernoulli beam has a deflection less than the Timoshenko beam for a small value of kAG=EI. Notice that plane sections remain plane for both beams. For the Euler Bernoulli beam, the sections perpendicular to the neutral axis remain perpendicular after deformation, while this is not the case for the Timoshenko beam. At higher values of (higher shear stiffness) the displacement of both beams are almost the same. Can you comment on what ratio of that would warrant the use of the Timoshenko beam approximation for this cantilever beam?
Euler Bernoulli Cantilever Beam:
The deflection of the centerline of a cantilever Euler Bernoulli beam loaded at the right end with a shearing force and a bending moment (See figure below) follows the following equation:
which can be obtained using the following Mathematica code:
View Mathematica Code
Clear[a,x,V,M,y,L,EI] M=EI*D[y[x],{x,2}]; V=EI*D[y[x],{x,3}]; th=D[y[x],x]; q=a; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; th1=th/.x>0; th2=th/.x>L; s=DSolve[{EI*y''''[x]==0,M2==Mom,V2==PP,y1==0,th1==0},y,x]
The deflected shape can be visualized for units and units with a beam height of 1 unit using the following tool with and
Timoshenko Cantilever Beam:
The deflection and cross section rotation of the centerline of a cantilever Timoshenko beam loaded at the right end with a shearing force and a bending moment (similar to the Euler Bernoulli beam above) follows the following equations:
which can be obtained using the following Mathematica code:
View Mathematica Code
Clear[a,x,V,M,y,L,psi,EI,kAG] M=EI*D[psi[x],x]; V=EI*D[psi[x],{x,2}]; q=a; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; psi1=psi[x]/.x>0; psi2=psi[x]/.x>L; s=DSolve[{EI*psi'''[x]==q,psi[x]==y'[x]+EI/kAG*psi''[x],y1==0,psi1==0,M2==Mom,V2==P},{y,psi},x]
The deflected shape can be visualized for units, units, and units with a beam height of 1 unit using the following tool with and
Accuracy of the Beam Approximation:
The beam approximation is used to analyze and design beams and frames under the effect of combined lateral loads, axial loads, and bending moments. In a coordinate system that includes the beam longitudinal axis, lateral loads lead to the development of “shear” stresses in the beam, and thus, the design codes include formulas for the design of beams under the effect of “bending” and “shear” stresses. When more detailed analysis of the stresses state inside a beam is required, finite element analysis of a beam could be more beneficial since it does not include the “plane sections remain plane” assumption and is able to more accurately characterize the stress state inside the beam. Figure 9 shows the arrows pointing in the direction of the maximum principal stresses obtained from a finite element analysis of a simple beam under lateral loading. The direction of the principal stresses shows that in the middle of the beam span, the maximum principal stresses are primarily horizontal and located at the bottom fibers of the beam. These stresses develop due to the higher bending loads at midspan of beams under lateral loading. The direction of the principal stresses closer to the support, however, are influenced by the “shear” components of the support reactions and thus tend to be oriented at an angle that is dependent on the interaction between the bending and the shear loads. In concrete beams, for instance, cracks develop in locations closer to the support that are oriented perpendicular to the directions of the maximum principal stresses (direction of the maximum tension), and thus, “shear” reinforcement are required closer to the supports. Figure 9 shows that the orientation of the stresses obtained from finite element analysis, which is a better approximation than the beam approximation, are similar to those expected from the beam approximation. The accuracy of the beam approximation is normally guaranteed provided the displacements are small, and the beam cross sectional dimensions are much smaller than the longitudinal direction.
This section presented the derivations for the classical small deformation plane beam theories. The extension to threedimensional deformations is straightforward. For large deformations, the equations presented here are no longer valid. When the deformations are large, the use of the small strains tensor could be erroneous (see the section on small strain). In addition, when the equilibrium equations were derived in this chapter, the horizontal forces and the moment equation were uncoupled. However, for large deformations, the horizontal forces can have large effects on the bending moment induced in a beam. The Euler Bernoulli beam and the Timoshenko beam, as described in this chapter, do not include the effect of axial forces on the bending moment and thus cannot predict the phenomenon of buckling. In addition, no reference was given to any of the traditional or modern methods (the stiffness method in particular)for solution of frame structures composed of different beams connected together.
Examples and Problems:
Example 1:Euler Bernoulli Simple Beam
Find the deflection, rotation, bending moment, and shearing force as a function of the beam length of an Euler Bernoulli simple beam assuming units, units and unit.
Solution:
The DSolve function in Mathematica is used to solve the differential equation of equilibrium with four boundary conditions and on both the left and right ends of the beam.
View Mathematica Code
Clear[a,x,V,M,y,L,EI] M=EI*D[y[x],{x,2}]; V=EI*D[y[x],{x,3}]; th=D[y[x],x]; q=a; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; th1=th/.x>0; th2=th/.x>L; s=DSolve[{EI*y''''[x]==q,M1==0,M2==0,y1==0,y2==0},y,x]; y=y[x]/.s[[1]]; th=D[y,x]; M=EI*D[y,{x,2}]; V=EI*D[y,{x,3}]; y/.x>L/2 EI=1; a=1; L=1; Plot[y,{x,0,L},BaseStyle>Directive[Bold,15],PlotRange>All,AxesLabel>{"x","y"},PlotStyle>{Black,Thickness[0.005]}] Plot[M,{x,0,L},PlotRange>All,BaseStyle>Directive[Bold,15],AxesLabel>{"x","M"},PlotStyle>{Black,Thickness[0.005]}] Plot[V,{x,0,L},BaseStyle>Directive[Bold,15],PlotRange>All,AxesLabel>{"x","V"},PlotStyle>{Black,Thickness[0.005]}] Plot[th,{x,0,L},PlotRange>All,BaseStyle>Directive[Bold,15],AxesLabel>{"x","th"},PlotStyle>{Black,Thickness[0.005]}]
Example 2:Timoshenko Simple Beam
Find the deflection, rotation, bending moment, and shearing force as a function of the beam length of a Timoshenko simple beam assuming units, units, unit and unit.
Solution:
The DSolve function in Mathematica is used to solve the differential equation of equilibrium with four boundary conditions and on both the left and right ends of the beam.
View Mathematica Code
Clear[a,x,V,M,y,L,psi,EI,kAG] M=EI*D[psi[x],x]; V=EI*D[psi[x],{x,2}]; q=a; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; psi1=psi[x]/.x>0; psi2=psi[x]/.x>L; s=DSolve[{EI*psi'''[x]==q,psi[x]==y'[x]+EI/kAG*psi"[x],M1==0,M2==0,y1==0,y2==0},{y,psi},x]; y=y/.s[[1,2]]; psi=psi/.s[[1,1]]; y[L/2] EI=1; a=1; L=1; kAG=1; Plot[y[x],{x,0,L},BaseStyle>Directive[Bold,15],PlotRange>All,AxesLabel>{"x","y"},PlotStyle>{Black,Thickness[0.005]}] Plot[psi[x],{x,0,L},PlotRange>All,BaseStyle>Directive[Bold,15],AxesLabel>{"x","Psi"},PlotStyle>{Black,Thickness[0.005]}] Plot[M,{x,0,L},PlotRange>All,AxesLabel>{"x","M"},PlotStyle>{Black,Thickness[0.005]}BaseStyle>Directive[Bold,15]] Plot[V,{x,0,L},PlotRange>All,AxesLabel>{"x","V"},PlotStyle>{Black,Thickness[0.005]}BaseStyle>Directive[Bold,15]]
Example 3:Bending and Shear Stresses in an Euler Bernoulli Cantilever Beam
Consider a cantilever Euler Bernoulli beam with Young’s modulus GPa. The beam is 5m long and has a rectangular cross section. The height of the beam is 0.5m, while the width is 0.25m.
If a load of 125kN is applied downwards at the free end, draw the contour plot of , , along the beam. Also, sketch the principal stress directions at the top fibers, the bottom fibers, the neutral axis and at quarter of the height of the beam at the support . Also, draw the vector plot of the displacement of the beam.
Solution:
The first step is to solve for the beam deflection by solving the differential equation of equilibrium with :
where is the moment of inertia and and are given. The constants can be solved using the boundary conditions:
Therefore, the displacement function is:
From which the moment and the shear can be obtained as:
The stresses along the beam length and depth can be obtained as functions of and as follows:
the vector of displacement of the beam is given by:
The required plots are shown below.
As shown in the plots, the bending stress is positive at the top (tension) and negative at the bottom, and increases in magnitude towards the support. The shearing stress component is constant at the neutral axis (Why?) throughout the length of the beam and decreases toward zero at the top and bottom fibers of the beam. The von Mises stress distribution follows the bending stress component, but gives a positive value in tension and in compression!
The eigenvalues and the eigenvectors of the stress matrix at the specified locations are obtained using Mathematica and are sketched below. The shear stresses at the top and bottom fibers are zero, and thus, the principal stresses are aligned with the chosen orthonormal coordinate system. At the location of the neutral axis of the beam, the bending stress component is zero, and the only nonzero component is , thus, the principal stresses at the neutral axis are aligned at 45 degrees. The direction of the principal stress with the highest magnitude at is along the vector {0.9993, –0.0374212,0}. At , the principal stress with magnitude –30.042MPa is aligned with the vector {0.9993,0.0374212,0}.
Clear[a,x,V,M,y,L,EI] M=EI*D[y[x],{x,2}]; V=EI*D[y[x],{x,3}]; th=D[y[x],x]; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; th1=th/.x>0; th2=th/.x>L; s=DSolve[{EI*y''''[x]==0,M2==0,V2==125,y1==0,th1==0},y,x] y=y[x]/.s[[1]] EI=Ee*Ii; u = {X2*D[y, x], y} M=FullSimplify[M/.s[[1]]] V=V/.s[[1]] Ee=20000000; b=0.25; t=0.5; L=5; Ii=b*t^3/12; s11=M*X2/Ii; Q=(t/2X2)*b*(t/4+X2/2) s12=V*Q/Ii/b; smatrix={{s11,s12,0},{s12,0,0},{0,0,0}}; FullSimplify[Chop[smatrix]]//MatrixForm VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]sigma[[2,2]])^2+(sigma[[2,2]]sigma[[3,3]])^2+(sigma[[3,3]]sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))]; vonmisesstress=VonMises[smatrix] ContourPlot[s11, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25, ContourLabels > True, PlotLabel > "sigma_11"] ContourPlot[s12, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25, ContourLabels > True, PlotLabel > "sigma_12"] ContourPlot[vonmisesstress, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25, ContourLabels > True, PlotLabel > "sigma_vM"] VectorPlot[u, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25, PlotLabel > "u"] stressmat1=smatrix/.{x>0,X2>0.25} Eigensystem[stressmat1] stressmat2=smatrix/.{x>0,X2>0.25/2} Eigensystem[stressmat2] stressmat3=smatrix/.{x>0,X2>0} Eigensystem[stressmat3] stressmat4=smatrix/.{x>0,X2>0.25/2} Eigensystem[stressmat4] stressmat5=smatrix/.{x>0,X2>0.25} Eigensystem[stressmat5] stressmat1=smatrix/.{x>1,X2>0.25} Eigensystem[stressmat1] stressmat2=smatrix/.{x>1,X2>0.25/2} Eigensystem[stressmat2] stressmat3=smatrix/.{x>1,X2>0} Eigensystem[stressmat3] stressmat4=smatrix/.{x>1,X2>0.25/2} Eigensystem[stressmat4] stressmat5=smatrix/.{x>1,X2>0.25} Eigensystem[stressmat5]
Example 4: Bending and Shear Stresses in an Euler Bernoulli Simple Beam
Consider an Euler Bernoulli simple beam with Young’s modulus . The beam is 8m long and has a rectangular cross section. The height of the beam is 0.5m, while the width is 0.25m. If a load of 125kN/m is distributed over the beam length, draw the contour plot of , , and along the beam. Also, draw the vector plot of the displacement of the beam.
Solution:
The first step is to solve for the beam deflection by solving the differential equation of equilibrium with :
where is the moment of inertia and and are given. The constants can be solved using the boundary conditions:
Therefore, the displacement function is:
From which the moment and the shear can be obtained as:
The stresses along the beam length and depth can be obtained as functions of and as follows:
the vector of displacement of the beam is given by:
The required plots are shown below.
As shown in the contour plots, the contribution of the shear stresses in in this beam is very small compared to the bending stress . This is natural due to the small ratio between the depth of the beam and its length. The longer the beam, the less significant the shear stresses are. The following is the Mathematica code utilized for the above calculations. View Mathematica Code
Clear[a,x,V,M,y,L,EI,Ee,Ii] M=EI*D[y[x],{x,2}]; V=EI*D[y[x],{x,3}]; th=D[y[x],x]; q=125; M1=M/.x>0; M2=M/.x>L; V1=V/.x>0; V2=V/.x>L; y1=y[x]/.x>0; y2=y[x]/.x>L; th1=th/.x>0; th2=th/.x>L; s=DSolve[{EI*y''''[x] ==q,M2==0,y2==0,y1==0,M1==0},y,x] y=y[x]/.s[[1]] EI=Ee*Ii; u={X2*D[y,x],y} L=8; M=FullSimplify[M/.s[[1]]] V=FullSimplify[V/.s[[1]]] Ee=20000000; b=0.25; t=0.5; Ii=b*t^3/12; s11=M*X2/Ii Q=(t/2X2)*b*(t/4+X2/2) s12=V*Q/Ii/b smatrix={{s11,s12,0},{s12,0,0},{0,0,0}}; FullSimplify[Chop[smatrix]]//MatrixForm VonMises[sigma_]:=Sqrt[1/2*((sigma[[1,1]]sigma[[2,2]])^2+(sigma[[2,2]]sigma[[3,3]])^2+(sigma[[3,3]]sigma[[1,1]])^2+6*(sigma[[1,2]]^2+sigma[[1,3]]^2+sigma[[2,3]]^2))]; vonmisesstress=VonMises[smatrix] ContourPlot[s11,{x,0,L},{X2,0.25,0.25},AspectRatio>0.25,ContourLabels>True,PlotLabel>"Sigma_11"] ContourPlot[s12,{x,0,L},{X2,0.25,0.25},AspectRatio>0.25,ContourLabels>True,PlotLabel>"Sigma_12"] ContourPlot[vonmisesstress,{x,0,L},{X2,0.25,0.25},AspectRatio>0.25,ContourLabels>True,PlotLabel>"Sigma_vM"] VectorPlot[u,{x,0,L},{X2,0.25,0.25},AspectRatio>0.25,PlotLabel>"Displacement Vector"]
Example 5:Beam under Axial Load with Displacement Boundary Conditions
A beam fixed at both ends is under a triangular axial load with a maximum value of 125 kN/m at . The beam’s length, width, and height are 5, 0.25, and 0.5 m, respectively. Young’s modulus is 20 GPa. Ignore the effect of Poisson’s ratio. Find the displacement function of this beam along with the stress distribution along the length of the beam. Draw the contour plot of and the vector plot of the displacement along the beam.
Solution:
The displacement vector function of the beam has the form:
The only nonzero strain component is :
Ignoring the effect of Poisson’s ratio, the only nonzero stress component is :
The distributed axial load follows the equation:
The equilibrium equation (Equation 12) for beams under axial loading can be used with constant to write:
The boundary conditions are:
Therefore, the solution to the above ordinary differential equation yields the following forms for , , and the normal force :
The required contour plots and the Mathematica code used are below.
View Mathematica Code
Clear[u,x,L,EA,Ee,A] EA=Ee*A; u1=u[x]/.x>0; u2=u[x]/.x>L; Normalf=EA*D[u[x]]; qnormal=125*x/L; s=DSolve[{EA*u''[x]==qnormal,u1==0,u2==0},u,x] u=u[x]/.s[[1]] s11=Ee*D[u,x] Normalf=Normalf/.s[[1]]; Ee=20000; A=0.25*0.5; L=5; disp={u,0} ContourPlot[s11,{x,0,5},{X2,0.25,0.25},AspectRatio>0.25,ContourLabels>True,PlotLabel>"sigma_11 (MPa)"] VectorPlot[disp,{x,0,5},{X2,0.25,0.25},AspectRatio>0.25,PlotLabel>"Displacement Vector"]
Example 6:Beam under Axial Load with Displacement and Stress Boundary Conditions
The shown beam has one fixed end and one free end. The beam is under a triangular axial load with a maximum value of 125kN/m at . The beam’s length, width, and height are 5, 0.25, and 0.5m, respectively. Young’s modulus is 20 GPa. Ignore the effects of Poisson’s ratio and find the displacement function of this beam along with the stress distribution. Draw the contour plots of and the vector plot of the displacement.
Solution:
The displacement vector function of the beam has the form:
The only nonzero strain component is :
Ignoring the effect of Poisson’s ratio, the only nonzero stress component is :
The distributed axial load follows the equation:
The equilibrium equation (Equation 12) for beams under axial loading can be used with constant to write:
The boundary conditions are:
Therefore, the solution to the above ordinary differential equation yields the following forms for , , and the normal force :
The required contour plots and the Mathematica code used are below. Compare the plots with the previous example.
View Mathematica Code
Clear[u,x,L,EA,Ee,A] EA=Ee*A; u1=u[x]/.x>0; BC2=D[u[x],x]/.x>L; Normalf=EA*D[u[x]]; qnormal=125*x/L; s=DSolve[{EA*u''[x]==qnormal,u1==0,BC2==0},u,x] u=u[x]/.s[[1]] s11=Ee*D[u,x] Normalf=Normalf/.s[[1]]; Ee=20000; A=0.25*0.5; L=5; disp={u,0} ContourPlot[s11,{x,0,5},{X2,0.25,0.25},AspectRatio>0.25,ContourLabels>True,PlotLabel>"sigma_11 (MPa)"] VectorPlot[disp,{x,0,5},{X2,0.25,0.25},AspectRatio>0.25,PlotLabel>"Displacement Vector"]
Example 7:Beam under Axial Load with Varying Cross Sectional Area
The shown beam has a varying circular cross sectional area with a radius units at the top varying linearly to unit at the bottom. The beam is used to carry a concentrated load and has a weight density . Young’s modulus units. Find the vertical displacement field of the beam.
Solution:
Since the beam has a varying cross sectional area, this induces a varying vertical load due to the beam’s own weight. The radius and the area vary with according to the equations:
Similar to the two previous examples, the only nonzero strain component is and ignoring Poisson’s ratio we get:
where is the unknown displacement component along the axis . The equilibrium equation (Equation 12) for beams under axial loading can be used with varying to write:
Substituting for in the above equation we get:
The boundary conditions are:
Mathematica can be used to solve the above differential equation to find the solution for as follows:
View Mathematica Code
Clear[L,A,Al,ro,Ee,u] L=4; A=Pi*(2x/L)^2; Al=A/.x>L; s=DSolve[{D[Ee*A*u'[x],x]+ro*A==0,u[0]==0,u'[L]==P/Ee/Al},u,x]
You can now copy and paste the following code into Mathematica to see the effect of changing the density and the external load on the displacement function produced for and . What is the combination of and that would produce an almost linear displacement field?
View Mathematica Code
Clear[u,x,Ee,P,L,A] Manipulate[Ee=1000; L=4; A=Pi*(2x/L)^2; Al=A/.x>L; s=DSolve[{D[Ee*A*u'[x],x]+ro*A==0,u[0] ==0,u'[L] ==P/Ee/Al},u,x]; s=u/.s[[1]]; Plot[s[x],{x,0,L}],{P,0,500},{ro,0,100}]
Problems:
 Use Mathematica to solve the Euler Bernoulli beam and the Timoshenko beam equations (shear, moment, rotation (slope), and deflection) for the beams shown in the following figure (assume values for the loads and material constants). Also, find the maximum deflection associated with the loads shown.
 Rewrite the final equations presented for the Euler Bernoulli beam and the Timoshenko beam assuming that the crosssectional area parameters and are not constant.
 Conduct a literature search to find the values of the shear coefficient of the Timoshenko beam for different crosssectional shapes. (Square, circle, I beam etc.)
 A beam is 6m long with a rectangular crosssectional area whose width and height are 0.25m and 0.5m, respectively. Young’s modulus is 20GPa. Ignore the effect of Poisson’s ratio. A distributed load of a value 125kN/m is applied vertically downward. Compare the following between a simple beam and a fixedends beam:
 Deflection, bending moment, and shearing force distribution along the beam.
 , , , and the hydrostatic stress (Draw the contour plots and comment on the difference).
 The shown beam has a rectangular crosssectional area with a width of 0.25m. The height of the beam varies linearly along its length as shown in the figure below. Young’s modulus is taken as 20GPa while the effect of Poisson’s ratio is ignored. Assume that the only nonzero component of the stress is and that plane sections perpendicular to the axis remain plane and perpendicular to the axis after deformation. Find the following:
 A generic form for the position and the displacement functions.
 The expression for the infinitesimal strain tensor components.
 The differential equation of equilibrium and its boundary conditions.
 Solve the differential equation of equilibrium to find the displacement of the beam at every point.
 Draw the contour plots of .
 Draw the vector plot of the displacement vector .
 The shown beam structure has a circular cross sectional area with radius of 2m at the bottom varying linearly to 1m at the top. The length of the structure is 4m. A crane is used to carry it from the top. Consider the weight density of the member as 10 N/m^3 and Young’s modulus of 1000 N/m^2. Ignore Poisson’s ratio and assume that the only nonzero component of the stress is the normal stress in the vertical direction and that it is uniformly distributed in the cross section. Find the displacement function of the beam, the total elongation or contraction, and the distribution of the nonzero normal stress component in the following cases:

 When the beam is resting on the ground.
 When the crane is carrying half the weight of the beam but the bottom of the beam is still on the ground.
 When the beam is suspended in the air.

 A cantilever beam with a rectangular cross section is subjected to a uniformly distributed load. The length of the beam is 5 m, its width and height are respectively 400 mm and 500 mm and the distributed load is equal to 15 KN/m. A rigid block is located under the free end of the beam and a gap equal to 10 mm exists between the free end and the rigid block. Consider the Young’s modulus as 20 GPa and find the following:
 the reaction force of the rigid block.
 the deflection of the neutral axis of the beam.
 the displacement field inside the beam.
 the stress component inside the beam. Draw the contour plots of .
Hi! I’ve been reading your blog for a while now and finally
got the courage to go ahead and give you a shout out from Porter Tx!
Just wanted to mention keep up the fantastic job!
Thank you so much!
Thanks Dr. ADEEb, for this helpful presentation.
Thank you!!!
How would the code mathematica for a EULER BERNOULLI CANTILEVER BEAM but with parabolic load?
Clear[a, x, V, M, y, L, EI]
M = EI*D[y[x], {x, 2}];
V = EI*D[y[x], {x, 3}];
th = D[y[x], x];
M1 = M /. x > 0;
M2 = M /. x > L;
V1 = V /. x > 0;
V2 = V /. x > L;
y1 = y[x] /. x > 0;
y2 = y[x] /. x > L;
th1 = th /. x > 0;
th2 = th /. x > L;
(*parabolic load*)
q = x^2;
s = DSolve[{EI*y””[x] == q, M2 == 0, V2 == 0, y1 == 0, th1 == 0}, y,
x]
y = y[x] /. s[[1]]
EI = Ee*Ii;
u = {X2*D[y, x], y}
M = FullSimplify[M /. s[[1]]]
V = V /. s[[1]]
Ee = 20000000;
b = 0.25;
t = 0.5;
Ii = b*t^3/12;
L = 5;
s11 = M*X2/Ii;
Q = (t/2 – X2)*b*(t/4 + X2/2)
s12 = V*Q/Ii/b;
smatrix = {{s11, s12, 0}, {s12, 0, 0}, {0, 0, 0}};
FullSimplify[Chop[smatrix]] // MatrixForm
VonMises[sigma_] :=
Sqrt[1/2*((sigma[[1, 1]] – sigma[[2, 2]])^2 + (sigma[[2, 2]] –
sigma[[3, 3]])^2 + (sigma[[3, 3]] – sigma[[1, 1]])^2 +
6*(sigma[[1, 2]]^2 + sigma[[1, 3]]^2 + sigma[[2, 3]]^2))];
vonmisesstress = VonMises[smatrix]
ContourPlot[s11, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25,
ContourLabels > True, PlotLabel > “sigma_11”]
ContourPlot[s12, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25,
ContourLabels > True, PlotLabel > “sigma_12”]
ContourPlot[vonmisesstress, {x, 0, 5}, {X2, 0.25, 0.25},
AspectRatio > 0.25, ContourLabels > True, PlotLabel > “sigma_vM”]
VectorPlot[u, {x, 0, 5}, {X2, 0.25, 0.25}, AspectRatio > 0.25,
PlotLabel > “u”]
stressmat1 = smatrix /. {x > 0, X2 > 0.25}
Eigensystem[stressmat1]
stressmat2 = smatrix /. {x > 0, X2 > 0.25/2}
Eigensystem[stressmat2]
stressmat3 = smatrix /. {x > 0, X2 > 0}
Eigensystem[stressmat3]
stressmat4 = smatrix /. {x > 0, X2 > 0.25/2}
Eigensystem[stressmat4]
stressmat5 = smatrix /. {x > 0, X2 > 0.25}
Eigensystem[stressmat5]
stressmat1 = smatrix /. {x > 1, X2 > 0.25}
Eigensystem[stressmat1]
stressmat2 = smatrix /. {x > 1, X2 > 0.25/2}
Eigensystem[stressmat2]
stressmat3 = smatrix /. {x > 1, X2 > 0}
Eigensystem[stressmat3]
stressmat4 = smatrix /. {x > 1, X2 > 0.25/2}
Eigensystem[stressmat4]
stressmat5 = smatrix /. {x > 1, X2 > 0.25}
Eigensystem[stressmat5]
For figure 1, shouldn’t grad(u) element (1,2) be zero?
Great work by the way!
Sorry, I don’t understand your question
Thank you for your good lecture.
Could you let me know how you plot the Eigen vector and Eigen values on Example 3????
Hello Samer!
I couldn’t thank you enough for the work you presented here.
I have a question and would be glad if you can help me with that. In the Timoshenko cantilever problem, the deflection that you provided has the first term (3EIPX1). I’m doing it manually and I get everything the same except this term as (6EIPX1). Could you please let me know what boundary conditions you applied and what’s the reason I am getting this term different?
Please email me your calculations and I can take a look.