Open Educational Resources

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, line-like 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:

  1. Plane sections perpendicular to the neutral axis before deformation stay plane and perpendicular to the neutral axis after deformation (Figure 1).
  2. The deformations are small.
  3. The beam is linear elastic isotropic and Poisson’s ratio effects are ignored.


Figure 1. Deformation assumption for beams. The black dot represents an arbitrary point before and after deformation. Plane sections perpendicular to the neutral axis remain plane and perpendicular to the neutral axis after deformation. The coordinates of the point after deformation can be obtained accordingly.

Figure 1. Deformation assumption for beams. The black dot represents an arbitrary point before and after deformation. Plane sections perpendicular to the neutral axis remain plane and perpendicular to the neutral axis after deformation. The coordinates of the point after deformation can be obtained accordingly.

The beam is assumed to lie such that its long axis is aligned with X_1. The deformation is controlled by a function y(X_1) 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 p before deformation are given by:

    \[ X=\left(\begin{array}{c}X_1^p\\X_2^p\\X_3^p\end{array}\right) \]

Ignoring any effect due to Poisson’s ratio, the coordinates of the point p after deformation (Figure 1) are given by:

    \[ x=\left(\begin{array}{c}X_1^p-X_2^p\sin{\theta}\\y+X_2^p\cos{\theta}\\X_3^p\end{array}\right) \]

The second assumption (small deformation) allows the approximations \theta\approx \sin{\theta}\approx \tan{\theta} and \cos{\theta}\approx 1. Therefore, the displacement function after incorporating the small deformations assumption is given by:

    \[ u=x-X=\left(\begin{array}{c}-X_2^p\theta\\y\\0\end{array}\right) \]

We can drop the superscript p since this applies to any arbitrary point. We also note that y'=\tan{\theta}\approx \theta. Therefore, the displacement function has the form:

    \[ u=x-X=\left(\begin{array}{c}-X_2\frac{\mathrm{d}y}{\mathrm{d}X_1}\\y\\0\end{array}\right) \]

we remind the reader that y is a function of X_1 only and so is its derivative. Therefore, The gradient of the displacement matrix is given by:

    \[ \nabla u=\left(\begin{matrix}-X_2\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} & -\frac{\mathrm{d}y}{\mathrm{d}X_1} & 0 \\ \frac{\mathrm{d}y}{\mathrm{d}X_1} & 0 & 0\\0&0&0\end{matrix}\right) \]

Therefore, the small strain matrix has the form:

    \[ \varepsilon=\frac{1}{2}\left(\nabla u + \nabla u^T\right)=\left(\begin{matrix}-X_2\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} & 0 & 0 \\0 & 0 & 0\\0&0&0\end{matrix}\right) \]

In essence, the deformation assumptions result in all of the strain components being zero except for \varepsilon_{11}. Furthermore, we will look carefully at the variation of \varepsilon_{11} on a given cross section. For a given cross section of the beam we have X_1 is equal to constant. Therefore, the quantity \frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} is constant on the cross section and therefore \varepsilon_{11} is linear in the vertical direction (Figure 2). The value of \varepsilon_{11} is equal to zero at the neutral axis X_2=0. The positive convention for y and for y' lead to a positive strain below the neutral axis with the maximum positive value \varepsilon_a 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 \varepsilon_b 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 \sigma_{11}=E\varepsilon_{11} (Figure 2). Therefore, the normal stress component distribution on the cross section is given by:

(1)   \begin{equation*} \sigma_{11}=-EX_2\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} \end{equation*}


The linear distribution of the longitudinal strain component \varepsilon_{11} and the longitudinal stress component \sigma{11} on the beam cross section.

Figure 2. The linear distribution of the longitudinal strain component \varepsilon_{11} and the longitudinal stress component \sigma_{11} on the beam cross section.

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 y(X_1) 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 M_3(X_1)=M(X_1), and the shearing force V_{12}(X_1)=V(X_1). The positive convention for these internal forces is shown in Figure 3. M(X_1) can be obtained by integrating the stress component \sigma_{11} on the cross section:

    \[ M=\int \int \! -X_2\sigma_{11} \,\mathrm{d}X_2\mathrm{d}X_3 \]

The negative sign is introduced to keep the positive convention for M causing compression on the top fibre of the beam. By substituting for \sigma_{11} using Equation 1, the moment can be calculated as:

(2)   \begin{equation*} M=\int \int \! EX_2^2\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} \,\mathrm{d}X_2\mathrm{d}X_3=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} \end{equation*}

where I is the moment of inertia of the beam calculated according to the formula:

    \[ I=\int \int \! X_2^2 \,\mathrm{d}X_2\mathrm{d}X_3 \]

Equations 1 and 2 can be used to find the relationship between the internal bending moment M and the stress component \sigma_{11} on the cross section:

(3)   \begin{equation*} \sigma_{11}=-\frac{MX_2}{I} \end{equation*}


Integration of Stress Components on the Beam's Cross Section to Obtain the Internal Forces

Figure 3. Integration of stress components on the beam’s cross section to obtain the internal forces

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 \Delta X_1 and perpendicular to the axis X_1 are considered. Then, by assuming that V and M are smooth functions of X_1, 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 X_2 upwards and given by q with units of load per beam length. Assuming the sum of forces in the vertical direction (direction of X_2) is equal to zero yields:

    \[ \frac{\partial V}{\partial X_1}=q \]

Similarly, assuming the sum of moments around X_3 is equal to zero yields:

    \[ \frac{\partial M}{\partial X_1}=V \]

Since V, q and M are only functions of X_1, the partial derivatives can be replaced with total derivatives and the equilibrium equations become:

(4)   \begin{equation*} \begin{split} \frac{\mathrm{d} V}{\mathrm{d} X_1}& =q\\ \frac{\mathrm{d} M}{\mathrm{d} X_1}& =V\end{split} \end{equation*}

We can substitute for M using Equation 2 so the equilibrium equation in terms of the displacement function y and the loading q becomes:

(5)   \begin{equation*} EI\frac{\mathrm{d}^4 y}{\mathrm{d} X_1^4} =q \end{equation*}

The solution for y has the following form:

    \[ EIy=\int\int\int\int \!q \,\mathrm{d}X_1 \mathrm{d}X_1 \mathrm{d}X_1 \mathrm{d}X_1 + C_1\frac{X_1^3}{6}+C_2\frac{X_1^2}{2}+C_3X_1+C_4 \]

Four boundary conditions are needed to fully solve for y and its derivatives. These boundary conditions can be given in terms of the displacement y, the rotation \theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}, the bending moment M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} or the shearing force V=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}. Note that if q=0 the beam can have a polynomial shape up to the third degree. For this reason, this formulation is sometimes termed: The cubic formulation.


Figure 4. Equilibrium Equations for the Euler Bernoulli and the Timoshenko beams subjected to transverse loading.

Figure 4. Equilibrium Equations for the Euler Bernoulli and the Timoshenko beams subjected to transverse loading.

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 V even though the shear strains are equal to zero. The shearing stresses associated with V can be calculated by first assuming that the shear stress \sigma_{12} is constant across the beam width (along the X_3 direction). Then, by considering a small slice in the beam with width of \Delta X_1 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 F_1 caused by integrating \sigma_{11}+\Delta\sigma_{11} on the left cross section of the shaded area. On the right we have F_2 caused by integrating \sigma_{11} on the right cross section of the shaded area. At the bottom we have F_3 caused by integrating \sigma_{21} on the bottom of the shaded area. Using Equation 3 that relates \sigma_{11} with M, the three forces can be calculated as follows:

    \[\begin{split} F_1 & =\int\int\! \left(\sigma_{11}+\Delta \sigma_{11}\right)\,\mathrm{d}X_2\mathrm{d}X_3=\int\int\! -\frac{M+\Delta M}{I}X_2\,\mathrm{d}X_2\mathrm{d}X_3\\ F_2 & =\int\int\! \sigma_{11} \,\mathrm{d}X_2\mathrm{d}X_3=\int\int\! -\frac{M}{I}X_2\,\mathrm{d}X_2\mathrm{d}X_3\\ F_3 & = \sigma_{21}\Delta X_1 b \end{split} \]

where b is the width of the cross section (the width in the direction of X_3). Equilibrium dictates:

    \[ F_1-F_2-F_3=0 \]

Therefore,

    \[ F_1-F_2-F_3=\int\int\! -\frac{\Delta M}{I}X_2\,\mathrm{d}X_2\mathrm{d}X_3-\sigma_{21}\Delta X_1 b=0 \]

We will use the equality: \Delta M\approx\frac{\mathrm{d}M}{\mathrm{d}X_1}\Delta X_1=V\Delta X_1. Therefore:

    \[ \int\int\! -\frac{V}{I}\Delta X_1 X_2\,\mathrm{d}X_2\mathrm{d}X_3-\sigma_{21}\Delta X_1 b=0 \]

Since V, I, and \Delta X_1 are not functions of X_2 or X_3, they can be brought outside the integral. Denote Q to be the moment of the area defined by:

    \[ Q=\int\int\! X_2\,\mathrm{d}X_2\mathrm{d}X_3 \]

Therefore:

    \[ \sigma_{21}=-\frac{VQ}{Ib} \]

Since \sigma_{21} is equal to \sigma_{12}, the shear stress on the cross section at the point in question is equal to:

(6)   \begin{equation*} \sigma_{12}=-\frac{VQ}{Ib} \end{equation*}

A similar proof can be found in the Mechanics ebook here.


Calculating the shear stress at the point represented by the black circle by considering the equilibrium of the shaded area in the figure

Figure 5. Calculating the shear stress at the point represented by the black circle by considering the equilibrium of the shaded area in the figure

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:

    \[ \frac{\partial \sigma_{11}}{\partial X_1}+\frac{\partial \sigma_{21}}{\partial X_2}+\frac{\partial \sigma_{31}}{\partial X_3}=0 \]

Notice that \sigma_{31} is assumed to be zero and that we are using the undeformed coordinates X_i since small deformations are assumed. Also, no body forces are assumed to be acting on the beam. By setting \sigma_{31}=0 and substituting for \sigma_{11} and \sigma_{21} using Equations 3 and 6 respectively, we get the following:

    \[ \frac{\partial \sigma_{11}}{\partial X_1}+\frac{\partial \sigma_{21}}{\partial X_2}=\frac{\partial \left(-\frac{MX_2}{I}\right)}{\partial X_1}-\frac{\partial \frac{VQ}{Ib}}{\partial X_2} \]

Assuming I is constant along X_1 and since V and I are functions only of X_1 we get:

    \[ \frac{\partial \sigma_{11}}{\partial X_1}+\frac{\partial \sigma_{21}}{\partial X_2}=-\frac{X_2}{I}\frac{\partial M}{\partial X_1}-\frac{V}{I}\frac{\partial \frac{Q}{b}}{\partial X_2} \]

Therefore:

    \[ \frac{\partial \sigma_{11}}{\partial X_1}+\frac{\partial \sigma_{21}}{\partial X_2}=\frac{V}{I}\left(-X_2 -\frac{\partial \frac{Q}{b}}{\partial X_2}\right) \]

The left hand side is only equal zero if b 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)

  1. Plane sections perpendicular to the neutral axis before deformations remain plane, but not necessarily perpendicular to the neutral axis after deformation (Figure 6).
  2. The deformations are small.
  3. The beam is linear elastic isotropic and Poisson’s ratio effects are ignored.


Figure 6. Timoshenko Beam deformation shape. The cross sections perpendicular to the neutral axis before deformation stay plane after deformation but are not necessarily perpendicular to the neutral axis after deformation.

Figure 6. Timoshenko Beam deformation shape. The cross sections perpendicular to the neutral axis before deformation stay plane after deformation but are not necessarily perpendicular to the neutral axis after deformation.

Similar to the Euler Bernoulli beam, the Timoshenko beam is assumed to lie such that its long axis is aligned with X_1. The deformation is controlled by two functions: y(X_1) which describes the deformation of the neutral axis (Figure 6) and \gamma(X_1) which describes the shear rotation of the cross section. In essence, the total rotation of the cross section is equal to y'(X_1)+\gamma(X_1). The coordinates of an arbitrary point p before deformation are given by:

    \[ X=\left(\begin{array}{c}X_1^p\\X_2^p\\X_3^p\end{array}\right) \]

Ignoring any effect due to Poisson’s ratio, the coordinates of the point p after deformation (Figure 6) are given by:

    \[ x=\left(\begin{array}{c}X_1^p-X_2^p\sin{(y'+\gamma)}\\y+X_2^p\cos{(y'+\gamma)}\\X_3^p\end{array}\right) \]

Setting \psi=y'+\gamma, the small deformations assumption allows the approximations \psi\approx \sin{\psi}\approx \tan{\psi} and \cos{\psi}\approx 1. Therefore, the displacement function after incorporating the small deformations assumption is given by:

    \[ u=x-X=\left(\begin{array}{c}-X_2^p\psi\\y\\0\end{array}\right) \]

We can drop the superscript p since this applies to any arbitrary point. We also note that y'=\tan{\theta}\approx \theta. Therefore, the displacement function has the form:

    \[ u=x-X=\left(\begin{array}{c}-X_2\psi\\y\\0\end{array}\right) \]

we remind the reader that y and \gamma (and hence \psi=y'+\gamma) are functions of X_1 only. Therefore, The gradient of the displacement matrix is given by:

    \[ \nabla u=\left(\begin{matrix}-X_2\frac{\mathrm{d}\psi}{\mathrm{d}X_1} & -\psi & 0 \\ \frac{\mathrm{d}y}{\mathrm{d}X_1} & 0 & 0\\0&0&0\end{matrix}\right) \]

Therefore, the small strain matrix has the form:

    \[ \varepsilon=\frac{1}{2}\left(\nabla u + \nabla u^T\right)=\left(\begin{matrix}-X_2\frac{\mathrm{d}^\psi}{\mathrm{d}X_1} & -\frac{\gamma}{2} & 0 \\-\frac{\gamma}{2} & 0 & 0\\0&0&0\end{matrix}\right) \]

In essence, the deformation assumptions result in all of the strain components being zero except for \varepsilon_{11} and \varepsilon_{12}=\varepsilon_{21}. Furthermore, we will look carefully at the variation of \varepsilon_{11} on a given cross section. For a given cross section of the beam we have X_1 is equal to constant. Therefore, the quantity \frac{\mathrm{d}\psi}{\mathrm{d}X_1} is constant on the cross section and therefore \varepsilon_{11} is linear in the vertical direction (Figure 2). The value of \varepsilon_{11} is equal to zero at the neutral axis X_2=0. The positive convention for \psi and for \psi' lead to a positive strain below the neutral axis with the maximum positive value \varepsilon_a 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 \varepsilon_b at the top fibre of the cross section. On the other hand \varepsilon_{12} 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 \sigma_{11}=E\varepsilon_{11} (Figure 2). Therefore, the normal stress component distribution on the cross section is given by:

(7)   \begin{equation*} \sigma_{11}=-EX_2\frac{\mathrm{d}\psi}{\mathrm{d}X_1}=-EX_2\left(\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}+\frac{\mathrm{d}\gamma}{\mathrm{d}X_1}\right) \end{equation*}

The internal bending moment on the cross section can be obtained using the integral:

(8)   \begin{equation*} M=\int\int\! -\sigma_{11}X_2 \,\mathrm{d}X_2\mathrm{d}X_3=\int\int\! E\frac{\mathrm{d}\psi}{\mathrm{d}X_1}X_2^2 \,\mathrm{d}X_2\mathrm{d}X_3=EI\frac{\mathrm{d}\psi}{\mathrm{d}X_1} \end{equation*}

The shear stress distribution on the cross section is given by:

    \[ \sigma_{12}=kG(2\varepsilon_{12})=-kG\gamma \]

where G is the shear modulus and k is a constant introduced to allow for averaging or smearing the shear stress over the full cross section. k was introduced by Timoshenko in his book. For more information on k you can view this reference. The shear force is the integral of this average \sigma_{12} over the cross section which yields:

    \[ V=\int\int \! -\sigma_{12} \, \mathrm{d}X_2\mathrm{d}X_3=\int\int \! kG\gamma \, \mathrm{d}X_2\mathrm{d}X_3=kA G\gamma \]

In essence, the term kA gives the corresponding “effective” shear area such that the shear stress is equal to \frac{V}{kA}.

Equilibrium Equations:

The equilibrium Equation 4 developed for the Eulber Bernoulli beam also applies to the Timoshenko beam (Figure 4). By substituting for M using Equation 8 in the equilibrium equations, we reach:

(9)   \begin{equation*} EI\frac{\mathrm{d}^3\psi}{\mathrm{d}X_1^3}=q \end{equation*}

We also have the following differential equation in terms of the beam displacement y:

    \[ y'=\psi-\gamma\Rightarrow y'=\frac{\mathrm{d}y}{\mathrm{d}X_1}=\psi - \frac{V}{kAG}\Rightarrow \]

(10)   \begin{equation*} \frac{\mathrm{d}y}{\mathrm{d}X_1}=\psi - \frac{EI}{kAG}\frac{\mathrm{d}^2\psi}{\mathrm{d}X_1^2} \end{equation*}

The solution for y and \psi have the following forms:

    \[\begin{split} EI\psi & =\int\int\int \!q \, \mathrm{d}X_1 \mathrm{d}X_1 \mathrm{d}X_1 + C_1\frac{X_1^2}{2}+C_2X_1+C_3\\ EI y & =\int\int\int\int \!q \, \mathrm{d}X_1 \mathrm{d}X_1 \mathrm{d}X_1\mathrm{d}X_1+C_1\frac{X_1^3}{6}+C_2\frac{X_1^2}{2}+C_3X_1+C_4\\ & -\frac{EI}{kAG}\left(\int\int \!q \, \mathrm{d}X_1 \mathrm{d}X_1+C_1X_1\right) \end{split} \]

The above two equations (9 and 10) can be solved if four boundary conditions for the cross section rotation \psi, the moment M=EI\frac{\mathrm{d}\psi}{\mathrm{d}X_1}, the shear V=EI\frac{\mathrm{d}^2\psi}{\mathrm{d}X_1^2}, and/or the displacement y are available. At least one boundary condition for y has to be given in order to find the integration constant C_4. Note that if kAG>>EI 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 axial-loading 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 u_1=u_1(X_1) 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:

    \[ x=\left(\begin{array}{c}X_1+u_1\\X_2\\X_3\end{array}\right) \]

In this case, the displacement function is:

    \[ u=x-X=\left(\begin{array}{c}u_1\\0\\0\end{array}\right) \]

Since u_1 is a function of only X_1, the only nonzero component of the strain is \varepsilon_{11}:

    \[ \varepsilon=\left(\begin{matrix}\frac{\mathrm{d}u_1}{\mathrm{d}X_1}&0&0\\0&0&0\\0&0&0\end{matrix}\right) \]

On any cross section, X_1 is constant and therefore the values of u_1 and its derivatives are constant. Therefore, \varepsilon_{11} is constant on any cross section perpendicular to X_1.


Deformation assumption of beams under axial loading. Cross section perpendicular to the neutral axis move horizontally in the direction o.

Figure 7. Deformation assumption of beams under axial loading. Cross section perpendicular to the neutral axis move horizontally in the direction of X_1 a distance u_1(X_1).

Equilibrium Equation:

The beam is assumed to be under a distributed normal load of value p of units of force/unit length. Figure 8 shows the applied distributed load along with the equilibrium of a beam slide of length \Delta X_1. The equilibrium equation can be written in terms of the normal force N which is the integral of the stress component \sigma_{11} over the cross sectional area:

    \[ N=\int\int \! \sigma_{11}\,\mathrm{d}X_2\mathrm{d}X_3 \]

Since \sigma_{11} is constant we have:

    \[ N=\sigma_{11}A \]

where A 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 N+\Delta N acting on the right cross section and directed to the right. The second force is N acting on the left cross section and directed to the left. The third is due to the distributed load p. Therefore:

    \[\begin{split} F_1&=N+\Delta N = \sigma_{11}A+\frac{\partial N}{\partial X_1}\Delta X_1=\sigma_{11}A+\frac{\partial (\sigma_{11}A)}{\partial X_1}\Delta X_1\\ F_2& =N=\sigma_{11}A\\ F_3&=p\Delta X_1 \end{split} \]

Equilibrium dictates:

    \[ F_1-F_2+F_3=0 \]

Therefore:

(11)   \begin{equation*} \frac{\partial (\sigma_{11}A)}{\partial X_1}+p=0 \end{equation*}

Assuming a linear elastic beam, the differential equation can be written in terms of the displacement by converting the stresses into strains:

    \[ \sigma_{11}=E\varepsilon_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

Therefore, Equation 11 can be written as follows:

(12)   \begin{equation*} \frac{\mathrm{d} \left(EA \frac{\mathrm{d}u_1}{\mathrm{d}X_1}\right)}{\mathrm{d} X_1}+p=0 \end{equation*}

Or equivalently, if E is constant:

    \[ E\frac{\mathrm{d} A}{\mathrm{d} X_1}\frac{\mathrm{d}u_1}{\mathrm{d}X_1}+EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2} +p=0 \]

In the case when the cross sectional area A is constant, the equilibrium equation take the simple form:

    \[ EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2} +p=0 \]


Figure 8. The beam is assumed to be under a distributed normal load p of units of load/distance. The equilibrium of a slice of the beam is used to write the equation of equilibrium.

Figure 8. The beam is assumed to be under a distributed normal load p of units of load/distance. The equilibrium of a slice of the beam is used to write the equation of equilibrium.

 

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 EI and kAG. 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 kAG (higher shear stiffness) the displacement of both beams are almost the same. Can you comment on what ratio of kAG/EI 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 P and a bending moment M (See figure below) follows the following equation:

    \[ y=\frac{3MX_1^2-3LPX_1^2+PX_1^3}{6EI} \]

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]
Beam2
The deflected shape can be visualized for EI=500 units and L=10 units with a beam height of 1 unit using the following tool with -2\leq P \leq 2 and -15\leq M \leq 15

Timoshenko Cantilever Beam:

The deflection y and cross section rotation \psi of the centerline of a cantilever Timoshenko beam loaded at the right end with a shearing force P and a bending moment M (similar to the Euler Bernoulli beam above) follows the following equations:

    \[ y=\frac{-3EIPX_1+3kAGMX_1^2-3kAGLPX_1^2+kAGPX_1^3}{6EIkAG}\qquad \psi=\frac{2MX-2LPX_1+PX_1^2}{2EI} \]

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 EI=500 units, kAG=0.01EI units, and L=10 units with a beam height of 1 unit using the following tool with -1\leq P \leq 1 and -10\leq M \leq 10

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 mid-span 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 three-dimensional 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.


The direction of the maximum principal stresses in a beam under lateral loading obtained using a commercial finite element analysis (ABAQUS 6.9).

Figure 9. The direction of the maximum principal stresses in a beam under lateral loading obtained using a commercial finite element analysis (ABAQUS 6.9).

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 L=1 units, q=-1 units and EI=1 unit.

Solution:

The DSolve function in Mathematica is used to solve the differential equation of equilibrium with four boundary conditions M=0 and y=0 on both the left and right ends of the beam.
Example1Beam
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 L=1 units, q=-1 units, kAG=1 unit and EI=1 unit.

Solution:

The DSolve function in Mathematica is used to solve the differential equation of equilibrium with four boundary conditions M=0 and y=0 on both the left and right ends of the beam.

Example2Beam
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 E=20 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 \sigma_{11}, \sigma_{12}, \sigma_{vM} 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 X_1=0. 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 q=0:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}X_1^4}=q=0\Rightarrow y=C_1\frac{X_1^3}{6}+C_2\frac{X_1^2}{2}+C_3X_1+C_4 \]

where I=\frac{bt^3}{12} is the moment of inertia and b and t are given. The constants C_i can be solved using the boundary conditions:

    \[\begin{split} @X_1=0 & :y=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\\ @X_1=5 & :V=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=125\\ @X_1=5 & :M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0 \end{split} \]

Therefore, the displacement function y is:

    \[ y=\frac{-125}{6EI}(3LX_1^2-X_1^3) \]

From which the moment and the shear can be obtained as:

    \[ M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=125(X_1-5) kN.m.\qquad V=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=125 kN \]

The stresses along the beam length and depth can be obtained as functions of X_1 and X_2 as follows:

    \[ \sigma_{11}=\frac{-MX_2}{I}=-\frac{125(X_1-5)X_2}{I}\qquad \sigma_{12}=\frac{-VQ}{Ib}=-\frac{125\left(\frac{t}{2}-X_2\right)\left(\frac{X_2}{2}+\frac{t}{4}\right)}{I} \]

the vector of displacement of the beam is given by:

    \[ u=x-X=\left(\begin{array}{c}-X_2\frac{\mathrm{d}y}{\mathrm{d}X_1}\\y\\0\end{array}\right) \]

The required plots are shown below.
Example3Beam

As shown in the plots, the bending stress \sigma_{11} is positive at the top (tension) and negative at the bottom, and increases in magnitude towards the support. The shearing stress component \sigma_{12} 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 \sigma_{11} is zero, and the only nonzero component is \sigma_{12}, thus, the principal stresses at the neutral axis are aligned at 45 degrees. The direction of the principal stress with the highest magnitude at X_2=0.125 m is along the vector {0.9993, –0.0374212,0}. At X_2=-0.125m, the principal stress with magnitude –30.042MPa is aligned with the vector {0.9993,0.0374212,0}.
Example3Beam2

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];
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/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]

Example 4: Bending and Shear Stresses in an Euler Bernoulli Simple Beam

Consider an Euler Bernoulli simple beam with Young’s modulus E=20 GPa. 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 \sigma_{11}, \sigma_{12}, and \sigma_{vM} 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 q=-125:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}X_1^4}=q=-125\Rightarrow EIy=-125\frac{X_1^4}{24}+C_1\frac{X_1^3}{6}+C_2\frac{X_1^2}{2}+C_3X_1+C_4 \]

where I=\frac{bt^3}{12} is the moment of inertia and b and t are given. The constants C_i can be solved using the boundary conditions:

    \[\begin{split} @X_1=0 & :y=0\\ @X_1=0 & :M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0\\ @X_1=8 & :y=0\\ @X_1=8 & :M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0 \end{split} \]

Therefore, the displacement function y is:

    \[ y=\frac{-125}{24EI}(L^3X_1-2LX_1^3+X_1^4) \]

From which the moment and the shear can be obtained as:

    \[ M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=-\frac{125X_1}{2}(X_1-8) kN.m.\qquad V=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=-125(X_1-4) kN \]

The stresses along the beam length and depth can be obtained as functions of X_1 and X_2 as follows:

    \[ \sigma_{11}=\frac{-MX_2}{I}=-\frac{-\frac{125X_1}{2}(X_1-8)X_2}{I}\qquad \sigma_{12}=\frac{-VQ}{Ib}=-\frac{-125(X_1-4)\left(\frac{t}{2}-X_2\right)\left(\frac{X_2}{2}+\frac{t}{4}\right)}{I} \]

the vector of displacement of the beam is given by:

    \[ u=x-X=\left(\begin{array}{c}-X_2\frac{\mathrm{d}y}{\mathrm{d}X_1}\\y\\0\end{array}\right) \]

The required plots are shown below.
BeamExample4
As shown in the contour plots, the contribution of the shear stresses \sigma_{12} in \sigma_{vM} in this beam is very small compared to the bending stress \sigma_{11}. 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/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,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 X_1= L. 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 \sigma_{11} and the vector plot of the displacement along the beam.
Axial1

Solution:

The displacement vector function of the beam has the form:

    \[ u=\left(\begin{array}{c}u_1(X_1)\\0\\0\end{array}\right) \]

The only nonzero strain component is \varepsilon_{11}:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

Ignoring the effect of Poisson’s ratio, the only nonzero stress component is \sigma_{11}:

    \[ \sigma_{11}=E\varepsilon_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

The distributed axial load follows the equation:

    \[ p=125\frac{X_1}{L} \]

The equilibrium equation (Equation 12) for beams under axial loading can be used with constant EA to write:

    \[ EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}=-125\frac{X_1}{L} \]

The boundary conditions are:

    \[ \begin{split} @X_1=0 & :u_1=0\\ @X_1=L & :u_1=0 \end{split} \]

Therefore, the solution to the above ordinary differential equation yields the following forms for u_1, \sigma_{11}, and the normal force N=\sigma_{11}A:

    \[ u_1=\frac{125(X_1L^2-X_1^3)}{6EAL}\qquad \sigma_{11}=\frac{125(L^2-3X_1^2)}{6AL}\qquad N=\frac{125(L^2-3X_1^2)}{6L} \]

The required contour plots and the Mathematica code used are below.
AxialSolution1

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 X_1=L. 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 \sigma_{11} and the vector plot of the displacement.
axial 2

Solution:

The displacement vector function of the beam has the form:

    \[ u=\left(\begin{array}{c}u_1(X_1)\\0\\0\end{array}\right) \]

The only nonzero strain component is \varepsilon_{11}:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

Ignoring the effect of Poisson’s ratio, the only nonzero stress component is \sigma_{11}:

    \[ \sigma_{11}=E\varepsilon_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

The distributed axial load follows the equation:

    \[ p=125\frac{X_1}{L} \]

The equilibrium equation (Equation 12) for beams under axial loading can be used with constant EA to write:

    \[ EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}=-125\frac{X_1}{L} \]

The boundary conditions are:

    \[ \begin{split} @X_1=0 & :u_1=0\\ @X_1=L & :\sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=0 \end{split} \]

Therefore, the solution to the above ordinary differential equation yields the following forms for u_1, \sigma_{11}, and the normal force N=\sigma_{11}A:

    \[ u_1=\frac{125(3X_1L^2-X_1^3)}{6EAL}\qquad \sigma_{11}=\frac{125(3L^2-3X_1^2)}{6AL}\qquad N=\frac{125(3L^2-3X_1^2)}{6L} \]

The required contour plots and the Mathematica code used are below. Compare the plots with the previous example.
AxialSolution2

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 r=2 units at the top varying linearly to r=1 unit at the bottom. The beam is used to carry a concentrated load P and has a weight density \rho. Young’s modulus E=1000 units. Find the vertical displacement field of the beam.
axial 3

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 r and the area A vary with X_1 according to the equations:

    \[ r=2-\frac{X_1}{L}\qquad A=\pi\left(2-\frac{X_1}{L}\right)^2 \]

Similar to the two previous examples, the only nonzero strain component is \varepsilon_{11} and ignoring Poisson’s ratio we get:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}\qquad \sigma_{11}=E\varepsilon_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1} \]

where u_1 is the unknown displacement component along the axis X_1. The equilibrium equation (Equation 12) for beams under axial loading can be used with varying EA to write:

    \[ E\frac{\mathrm{d}A}{\mathrm{d}X_1}\frac{\mathrm{d}u_1}{\mathrm{d}X_1}+EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}+\rho A=0 \]

Substituting for A in the above equation we get:

    \[ E\frac{\pi}{2}\left(\frac{X_1}{4}-2\right)\frac{\mathrm{d}u_1}{\mathrm{d}X_1}+E\pi\left(2-\frac{X_1}{4}\right)^2\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}+\pi\rho\left(2-\frac{X_1}{4}\right)^2=0 \]

The boundary conditions are:

    \[ \begin{split} @X_1=0 & :u_1=0\\ @X_1=L & :\sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=\frac{P}{A}\bigg|_{X_1=L} \end{split} \]

Mathematica can be used to solve the above differential equation to find the solution for u_1 as follows:

    \[ u_1=\frac{-12PX_1-112\pi\rho X_1+24 \pi\rho X_1^2-\pi\rho X_1^3}{6E\pi(X_1-8)} \]

View Mathematica Code
Clear[L,A,Al,ro,Ee,u]
L=4;
A=Pi*(2-x/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 0\leq P \leq 500 and 0\leq \rho \leq 100. What is the combination of \rho and P 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*(2-x/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:

  1. 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.
    ExamplesBeam1
  2. Rewrite the final equations presented for the Euler Bernoulli beam and the Timoshenko beam assuming that the cross-sectional area parameters I and A are not constant.
  3. Conduct a literature search to find the values of the shear coefficient k of the Timoshenko beam for different cross-sectional shapes. (Square, circle, I beam etc.)
  4. A beam is 6m long with a rectangular cross-sectional 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 fixed-ends beam:
    1. Deflection, bending moment, and shearing force distribution along the beam.
    2. \sigma_{11}, \sigma_{12}, \sigma_{vM}, and the hydrostatic stress (Draw the contour plots and comment on the difference).
  5. The shown beam has a rectangular cross-sectional 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 \sigma_{11} and that plane sections perpendicular to the X_1 axis remain plane and perpendicular to the X_1 axis after deformation. Find the following:
    1. A generic form for the position and the displacement functions.
    2. The expression for the infinitesimal strain tensor components.
    3. The differential equation of equilibrium and its boundary conditions.
    4. Solve the differential equation of equilibrium to find the displacement of the beam at every point.
    5. Draw the contour plots of \sigma_{11}.
    6. Draw the vector plot of the displacement vector u=\left(\begin{array}{c}u_1\\u_2\end{array}\right).

    ExampleBeam2

  6. 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:
    1. When the beam is resting on the ground.
    2. When the crane is carrying half the weight of the beam but the bottom of the beam is still on the ground.
    3. When the beam is suspended in the air.

  7. A cantilever beam is subjected to a uniformly distributed load. A rigid block is used to support the free end if the displacement of the free end reaches 0.1 units. Consider the uniformly distributed load as 10 units, Young’s modulus as 100 units, moment of inertia as 1000 units, length of the beam as 10 units.

    1. Find the reaction force of the rigid block.
    2. Find the deflection of the neutral axis of the beam.
    3. Draw the displacement vector field inside the beam.
    4. Draw the contour plot of the stress component \sigma_{11} inside the beam.

Leave a Reply

Your email address will not be published.