Analysis of Deep to Moderately Deep Beams Using Timoshenko Theory and 2D Plane-Stress Elasticity

1. Introduction

A Google Scholar search for the term “Timoshenko beam” currently yields approximately 127,000 results—an increase of 1.6 times compared to the 78,000 entries reported in 2020 by Elishakoff [1].

Although much of the early research on the Timoshenko beam theory focused on vibrations [2,3] and beams on elastic foundations [4], a substantial body of work has since addressed static analysis, including displacement and stress fields, as documented in a recent monograph [5]. According to Conway et al. [6], the analysis of deep beams can be broadly classified into three categories: (i) infinite span with periodic loading, (ii) infinite span with non-periodic loading, and (iii) spans of finite length. In most cases, the Airy stress function is expressed as an infinite series of hyperbolic functions or Fourier integrals. Alternatively, the principle of least work applied to the strain energy functional is sometimes employed [6,7]. Using the Schwarz–Christoffel transformation, Theocaris [8] derived an analytical solution for the stress distribution in a semi-infinite strip subjected to a concentrated axial load, represented by a few terms of hyperbolic functions rather than an infinite series. The same author later extended the approach to cases involving distributed loads [9]. For the semi-infinite strip, Benthem [10] solved the clamped-end case by applying the Laplace transform to eliminate the longitudinal coordinate x from the differential equation governing the Airy stress function.

In the context of mixed boundary-value problems, Horvay and Born [11] proposed solutions in the form of infinite series composed of exponential and trigonometric functions for a semi-infinite elastic strip. Their work addressed two specific cases: (a) a prescribed quadratic shear displacement with zero normal stress, and (b) a prescribed cubic normal displacement with zero shear stress. Expanding on this approach, Johnson and Little [12] employed a series expansion involving ten hyperbolic terms to address all possible combinations of prescribed stresses and/or displacements along the short edge of the strip. This includes the pure displacement case, commonly referred to as the displacement problem.

Shear coefficients—defined as the ratio of the average shear strain over a cross-section to the shear strain at the centroid—have been computed by Cowper [13] for various cross-sectional geometries, including the rectangular section. For rectangular beams, Labuschagne et al. [14] conducted a comparative study of linear beam theories, namely Euler–Bernoulli, Timoshenko, and two-dimensional elasticity theory. For the same cross-sectional type, Levinson [15] proposed a shear coefficient of $5/6$, a value that falls within the commonly reported range of 0.822 to 0.870 found in the literature for rectangular cross sections.

The present paper was primarily motivated by the fact that, within the context of computational mechanics, the exact deflection and associated stresses in an end-loaded cantilever are widely used as a benchmark test. Nevertheless, a question that persists to this day concerns the exact value of the maximum transverse displacement (also referred to as deflection). More specifically, the commonly used closed-form expression proposed by Timoshenko and Goodier [16] assumes a linear distribution of the normal stress ${\sigma }_x$ and a parabolic distribution of the shear stress ${\tau }_{xy}$.

To prevent rigid-body motion, three kinematical constraints were imposed: two displacement components at the midpoint of the fixed end and the axial displacement slope at the same location. As a result, since the slope of the neutral axis is non-zero, this model allows for shear strain at the support. Consequently, the resulting analytical expression for total deflection includes contributions from both bending and shear forces.

However, it is well known that this analytical solution is characterized by warping of the displacement field near the built-in end. Due to this effect, the support region will henceforth be referred to as the fixed section, where the quotation denotes the presence of warping (with displacement components: $U\ne 0$ and $V\ne 0$) rather than complete constraint. In contrast, for the case where the beam’s short edge is fully fixed ($U=V=0$), we adopt the term clamped beam in the present study.

According to the above explanations, the analytical solution given by Timoshenko and Goodier [16] does not exactly describe the example of a horizontally oriented cantilever beam in which the right vertical edge is fully constrained (clamped). For this reason, the boundary conditions necessary to match the exact solution must be followed, otherwise the FEM numerical solution converges to a smaller value [17]. The latter issue will be thoroughly discussed in this paper.

Alternatively, Livesley [18] modifies Timoshenko and Goodier [16] as follows. In addition to the two zero displacement components at the middle of the ‘fixed’ section, he also reduces the horizontal displacement component at the upper end. Then he found that the lower end of the fixed section is also fixed. Note that this issue is also cited as Exercise-5 in [16]. In his textbook, Livesley [18] uses his formula as a reference to compare the quality of one-dimensional finite elements for the deep beam. But the question is how accurate the latter analytical formula is.

Readers may also refer to Roark’s Formulas, which state that the shear displacement, $V_s$, can be determined using standard methods such as the method of unit loads or Castigliano’s theorem. This result, however, requires multiplication by a form factor F, which is 6/5 for a rectangular cross-section [19].

The purpose of this paper is to clarify the mechanics of Timoshenko’s beam, focusing on the following issues:

• To evaluate the accuracy of the calculated deflection obtained using the correcting factor $F=\frac{6}{5}$ (as reported in Ref. [19]).
• To investigate which of the established formulas provide lower or upper bounds for the finite element solution.
• To propose a practical algorithm and accompanying computer code for calculating the maximum deflection of a clamped beam subjected to a tip load.

To enhance reader comprehension, and despite Professor Timoshenko’s contributions to all three formulas investigated herein, we will refer to them distinctly as Roark’s, Timoshenko’s, and Livesley’s formulas. A comprehensive discussion can be found in Section 5.

2. Energy Conservation

Let us consider a rectangular beam of height h in the y-direction while its length in the x-direction is L. This cantilever has a narrow rectangular cross section of width b (not shown) bent by a force P applied at the end O (Figure 1).

Following Timoshenko and Goodier [16], we adopt stress distribution according to Bernoulli theory, thus the normal stresses ${\sigma }_x$ and ${\sigma }_y$ follow the linear law:

$${\sigma }_x={\displaystyle \frac{M}{I}}y=-{\displaystyle \frac{P}{I}}xy,\mathrm{and}{\sigma }_y=0,,$$

(1)

while the shear stress follows the well-known parabolic law:

$${\tau }_{xy}=-{\displaystyle \frac{6P}{b{h}^3}}\left({\displaystyle \frac{h^2}{4}}-y^2\right)=-{\displaystyle \frac{P}{2I}}\left({\displaystyle \frac{h^2}{4}}-y^2\right).$$

(2)

The abovementioned second moment of inertia of the cross section is given by:

$$I={\displaystyle \frac{b{h}^3}{12}}.$$

(3)

In general, the sign in the normal stress ${\sigma }_x$ is positive when tension occurs. The lateral force P is taken upwards (i.e., in the positive direction of y-axis), thus the upper fibers ($y>0$) are in compression (${\sigma }_x<0$). Therefore, since the associated y coordinate of these fibers is positive ($y=+h/2$) a minus sign becomes necessary to be set in Eq. (1). Moreover, the usual convention of positive shear stress ${\tau }_{xy}$ dictates a negative sign in Eq. (2) as well. The strain energy of the rectangular sheet (of dimensions $L\times b\times h$) is estimated by the volume integral of the strain energy density:

$$W={\int }_V\left({\sigma }_x{ϵ}_x+{\sigma }_y{ϵ}_y+{\tau }_{xy}{\gamma }_{xy}\right)dV.$$

(4)

To determine the strain ${ϵ}_x$ involved in Eq. (4) we must consider the material constitutive law for the plane-stress state,

$$\left[\begin{array}{c}{\sigma }_x\\ {\sigma }_y\\ {\tau }_{xy}\end{array}\right]={\displaystyle \frac{E}{1-{\nu }^2}}\left[\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& {\displaystyle \frac{(1-\nu )}{2}}.\end{array}\right]$$

(5)

According to Eq. (1) we have ${\sigma }_y=0$, thus eventually Eq. (5) gives:

$${\sigma }_x=E{ϵ}_x,\mathrm{and}{\tau }_{xy}=G{\gamma }_{xy}.$$

(6)

Taking the variables in the intervals $0\le x\le L$, $-\frac{h}{2}\le y\le \frac{h}{2}$, and $0\le z\le b$, and substituting Eqs. (1), (2) and (6) into Eq. (4), the strain energy splits in two volume integrals which can be analytically written as follows:

$$W_b={\displaystyle \frac{1}{2E}}{\int }_V{\sigma }_x^{2}dxdydz={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{P^2}{E}}{\displaystyle \frac{4L^3}{b{h}^3}}\right)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{P^2}{E}}{\displaystyle \frac{L^3}{3I}}\right),$$

(7)

$$W_s={\displaystyle \frac{1}{2G}}{\int }_V{\tau }_{xy}^{2}dxdydz={\displaystyle \frac{1}{2}}{\displaystyle \frac{P^{2}L}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right].$$

(8)

Obviously, the total strain energy $W=W_b+W_s$ due to bending and shear components, must be equal to the external work of the shear force P plus the work $W_{\mathrm{support}}$ done by the reaction stresses at the warped supporting boundary (`fixed’ section), thus we have:

$${\displaystyle \frac{1}{2}}P{V}_0+W_{\mathrm{support}}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{P^2{L}^3}{3EI}}\right)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{P^{2}L}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right].$$

(9)

For the rectangular sheet under study, the obtained factor $F=\frac{6}{5}$ in Eq. (9) is a standard for a cantilever beam. It has been widely written that it comes from the parabolic profile of the shear stress described by Eq. (2).

In more details, the well-known Shear Correction Factor (SCF) $k_s$ is defined as a factor which when applied to the actual cross-sectional area A produces an artificial area $A_s=k_{s}A$ associated to a virtual constant shear stress ${\tau }_{xy}^{\prime }=P/A_s$ which yields the same strain energy as the actual strain distribution. For a rectangle in which the actual strain distribution is parabolic, i.e., ${\tau }_{xy}={\displaystyle \frac{6P}{b{h}^3}}({\displaystyle \frac{h^2}{4}}-y^2)$, it is trivial to validate that $k_s=5/6$ (see, e.g., [5]). Therefore, by replacing the actual cross-sectional area ($A=bh$) with the equivalent area $A_s=k_{s}A$, the inverse $F=1/k_s=6/5$ suddenly appears as a factor in front of the classical shear deflection ${\displaystyle \frac{PL}{AG}}$ if the whole rectangular beam is considered as a sole elastic element. To be more accurate, the latter occurs for every slice of dimensions `thickness × height’ = $\Delta x\times h$, but thanks to the constant value of shear force $Q=P$ to every uniform cross section we can apply the factor F to the average shear $\tau =P⁄A_s=G\gamma =G\left({\displaystyle \frac{V_s}{L}}\right)$ thus solving the latter equation for $V_s$ we have:

$$V_s={\displaystyle \frac{PL}{A_{s}G}}={\displaystyle \frac{PL}{\left(k_{s}A\right)G}}={\displaystyle \frac{1}{k_s}}\left({\displaystyle \frac{PL}{AG}}\right)=F\left({\displaystyle \frac{PL}{AG}}\right).$$

(10)

Equation (10), which was derived from scratch, is also suggested by Roark’s formulas [12].

One may observe that the last term in the right-hand side of Eq. (9), which is the shear work $W_s$ of Eq. (8), refers to the abovementioned work $\frac{1}{2}P{V}_s$. However, due to the still unknown work $W_{\mathrm{support}}$ in the left-hand side, we cannot yet solve Eq. (9) for the total deflection $V_0$. Nevertheless, if we neglect the work $W_{\mathrm{support}}$ done by the reaction stresses at the supporting boundary, solving Eq. (9) for the maximum deflection $V_0$, we receive the well accepted approximation:

$$\mathrm{Roark}\text{’}\mathrm{s}\mathrm{formula}\text{:}\begin{array}{|c|}\hline {\left(V_0\right)}_{\mathrm{approx}}\approx {\displaystyle \frac{P{L}^3}{3EI}}+{\displaystyle \frac{PL}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right]\\ \hline\end{array}.$$

(11)

In other words, the maximum approximate deflection ${\left(V_0\right)}_{\mathrm{approx}}$ in Eq. (11) consists of the classical Bernoulli solution $V_b=P{L}^3/\left(3EI\right)$ plus the shear term:

$$V_s={\displaystyle \frac{PL}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right].$$

(12)

Again, Eq. (12)–which was obtained here from scratch–coincides with that proposed by Roark’s formulas [19].

On the other hand, Timoshenko and Goodier (see, [16], Eq. (r) on page 46 therein]) promote a larger shear value according to the formula:

$$V_s^{\prime }={\displaystyle \frac{P(h^2/4)L}{2IG}}={\displaystyle \frac{PL}{\left(bh\right)G}}\left[{\displaystyle \frac{3}{2}}\right].$$

(13)

Comparing Eq. (12) with Eq. (13), one may observe that the calculated shear deflections are quite different ($V_s<V_s^{\prime }$) unless the shear modulus G in Timoshenko & Goodier’s formula, expressed by Eq. (13), is technically replaced by:

$$G_s={\displaystyle \frac{5}{4}}G=1.25G.$$

(14)

The reason for the discrepancy between Eq. (12) and (13) is that the boundary conditions imposed in [16] lead to relatively large warping, and thus the associated mechanical work, $W_{\mathrm{support}}$, cannot be neglected.

3. Analytical Solutions

To calculate the mechanical work $W_{\mathrm{support}}$ done by the reaction stresses at the supporting `fixed’ boundary, we need to determine the displacement field in the beam and particularly along the `fixed’ section AMB (see, Figure 1).

3.1. Displacement Field

First it is simple to verify that the stresses of Eqs. (1) and (2) satisfy the well-known equilibrium conditions:

$${\displaystyle \frac{\partial {\sigma }_x}{\partial x}}+{\displaystyle \frac{\partial {\tau }_{xy}}{\partial y}}=0,\mathrm{and}{\displaystyle \frac{\partial {\tau }_{xy}}{\partial x}}+{\displaystyle \frac{\partial {\sigma }_y}{\partial y}}=0.$$

(15)

Considering the material constitutive law given by Eq. (5), the stress equilibrium equations are converted into the so-called Navier-Kichhoff equations in terms of only the displacement components:

$${\displaystyle \frac{{\partial }^{2}U}{\partial x^2}}+{\displaystyle \frac{1-\nu }{2}}{\displaystyle \frac{{\partial }^{2}U}{\partial y^2}}+{\displaystyle \frac{1+\nu }{2}}{\displaystyle \frac{{\partial }^{2}V}{\partial x\partial y}}=0,,$$

$${\displaystyle \frac{{\partial }^{2}V}{\partial x^2}}+{\displaystyle \frac{1-\nu }{2}}{\displaystyle \frac{{\partial }^{2}V}{\partial y^2}}+{\displaystyle \frac{1+\nu }{2}}{\displaystyle \frac{{\partial }^{2}U}{\partial x\partial y}}=0,.$$

(16)

Integrating the Navier-Kichhoff equations (16) gives the general solution:

$$\mathbf{u}=\left[\begin{array}{c}U(x,y)\\ V(x,y)\end{array}\right]={\displaystyle \frac{P}{EI}}\left[\begin{array}{c}{\displaystyle \frac{x^{2}y}{2}}-y^3{\displaystyle \frac{(2+\nu )}{6}}-C_{1}y+C_2\\ -{\displaystyle \frac{\nu x{y}^2}{2}}-{\displaystyle \frac{x^3}{6}}+[h^2{\displaystyle \frac{(1+\nu )}{4}}+C_1]x+C_3\end{array}\right].$$

(17)

with $C_1$, $C_2$ and $C_3$ being constants of integration. Note that since the boundary conditions of this problem are of Neumann-type (prescribed tractions), the solution is not unique unless we prevent the rigid-body motion. But even so, the solution for the displacement field (not the stresses) depends on the selected points at which boundary conditions of Dirichlet type are imposed. In an analogous way, prescribed flux in a closed domain governed by Laplace/Poisson equation does not ensure a unique solution unless the potential is prescribed at least to one boundary point. Again, now that the problem is plane stress elasticity, we need the restriction of three degrees of freedom to prevent rigid-body motion of the rectangular beam. This constraint can be realized in various ways, which however lead to different displacement fields thus to different warping patterns. Below we shall study two typical cases:

3.1.1. Case-1: Timoshenko & Goodier [16]

To prevent the abovementioned rigid-body motion, Timoshenko & Goodier [16] impose zero displacement components at the middle point M ($U_M=V_M=0$) of the fixed section as shown in Figure 1, and hence we have $C_2=0$ and $C_3=L^3/6-[h^2(1+\nu )/4+C_1]L$. In addition, they also fix a vertical segment of the cross section at the same point M($x=L,y=0$). Hence, this particular condition ${(\partial U/\partial x)}_{x=L,y=0}=0$, leads to $C_1=L_2/2$, so $C_3=-L^3/3-h^{2}L(1+\nu )/4$. The results are also shown in the first line of Table 1.

Remark: It should become clear that Timoshenko & Goodier [16] eventually applied the condition ${(\partial U/\partial x)}_{x=L,y=0}=0$. Previous to that, the showed that the alternative condition ${(\partial V/\partial x)}_{x=L,y=0}=0$ leads to the well-known value $P{L}^3/\left(3EI\right)$ (Bernoulli solution), usually derived in elementary books on the strength of materials.

Based on the constants shown in the second row of Table 1 (labelled as Timoshenko & Goodier [16]), Eqs. (17) give the following displacements (note that the subscript ‘T’ stands for Timoshenko):

$$\mathrm{Maximum}\mathrm{deflection}\text{:}{V}_T(0,0)={\displaystyle \frac{P{L}^3}{3EI}}\left[1+3(1+\nu ){\displaystyle \frac{h^2}{\left(4L^2\right)}}\right],$$

(18)

which may also be written as:

$$\begin{array}{cc}\hfill V_T(0,0)& ={\displaystyle \frac{P{L}^3}{3EI}}+{\displaystyle \frac{P{(h/2)}^{2}L}{2GI}}\hfill \\ \hfill & ={\displaystyle \frac{P{L}^3}{3EI}}+{\displaystyle \frac{3PL(1+\nu )}{Ebh}},\hfill \end{array}$$

(19)

Displacement component at `fixed’ section in x-direction:

$$U_T(L,y)={\displaystyle \frac{(2+\nu )P}{6EI}}{y}^3.$$

(20)

Displacement component at `fixed’ section in y-direction:

$$V_T(L,y)={\displaystyle \frac{\nu PL}{2EI}}{y}^2.$$

(21)

3.1.2. Case-2: According to Livesley [18]

Following Livesley [18], the boundary conditions are modified as follows. The point M($x=L,y=0$) in the middle of support section AB is fully fixed and also the upper point A of the same section (see, Figure 1) is restrained to move horizontally, thus we have:

$$\begin{array}{cc}\hfill \left(\mathrm{i}\right)\mathrm{Point}\mathrm{M}\text{:}& U=V=0\mathrm{at}(x=L,y=0),\hfill \\ \hfill \left(\mathrm{ii}\right)\mathrm{Point}\mathrm{A}\text{:}& U=0\mathrm{at}(x=L,y=h/2).\hfill \end{array}$$

(22)

These conditions determine the arbitrary constants in Eq. (17) as follows, $C_1={\displaystyle \frac{L^2}{2}}-{\displaystyle \frac{(2+\nu )h^2}{24}},C_2=0,C_3=-{\displaystyle \frac{L^3}{3}}-{\displaystyle \frac{h^{2}L(4+5\nu )}{24}}$, eventually giving:

$$\mathbf{u}=\left({\displaystyle \frac{P}{2EI}}\right)\left[\begin{array}{c}y(x^2-L^2)+{\displaystyle \frac{(2+\nu )y\left(\frac{h^2}{4}-y^2\right)}{3}}\\ {\displaystyle \frac{(L^3-x^3)}{3}}+\left[L^2+{\displaystyle \frac{h^2(4+5\nu )}{12}}\right](x-L)-\nu x{y}^2\end{array}\right].$$

(23)

Note that, because of symmetry, U is also zero at the lower point B $(x=L,y=-h/2)$ shown in Figure 1. Thus, the condition $U=0$ applies to both extreme points A and B at $y=\pm h/2$.

The vertical displacement of the tip O of the cantilever is given by

$$\begin{array}{cc}\hfill V_L(0,0)& =\left({\displaystyle \frac{P{L}^3}{3EI}}\right)\left[1+(2+2.5\nu ){\displaystyle \frac{h^2}{4L^2}}\right]\hfill \\ \hfill & =\left({\displaystyle \frac{P{L}^3}{3EI}}\right)+{\displaystyle \frac{(4+5\nu )PL}{2Ebh}}.\hfill \end{array}$$

(24)

Moreover, regarding the warping along the fixed edge the useful formulas are given by

$$\begin{array}{cc}\hfill x-\mathrm{direction}\text{:}& U_L(L,y)=-{\displaystyle \frac{(2+\nu )P}{6EI}}y(h^2/4-y^2)\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill y-\mathrm{direction}\text{:}& V_L(L,y)={\displaystyle \frac{\nu PL}{2EI}}{y}^2.\hfill \end{array}$$

(26)

3.2. Comparison Between the Two Cases

So far, we have seen that none of the above two cases, `T’ and `L’, ensures full fixation of the initially supposed `fixed’ section but eventually found to be warped.

The warped profiles in the x-direction are compared in Figure 2 (points A,M,B are the same as those illustrated in Figure 1, but rotated in clockwise direction by ${90}^{\circ }$). There, one may observe that in Case 1 (Timoshenko & Goodier: Eqs. (18) to (21)) the warping vanishes only at the middle point M and is opposite to the existing tension and compression associated to parts MB and MA, respectively (clockwise reaction moment). In contrast, in Case 2 (Livesley: Eqs. (24) to ()) the warping is more reasonable (vanishes at three points, A, M and B) but is consistent to the dominating tension and compression (anti-clockwise reaction moment). In both cases, Eqs. (20) and (25) suggest that, although the expressions are different, the variation of the longitudinal $U\left(y\right)$-displacement with respect to the offset y is cubic ($U\propto y^3$) on the clamp.

Regarding the warping profiles in the y-direction, Eq. (21) coincides with Eq. (), and thus in both cases we receive identical vertical (transverse) displacements ($V\left(y\right)$) of quadratic shape ($V\propto y^2$) on the clamp.

4. Work Associated with the Warping

Below we provide analytical formulas in both cases, ‘T’ and ‘L’.

4.1. Timoshenko & Goodier Formula [16] (`T’)

The warped section at $(x=L)$ is described by the displacement components $U_T$ and $V_T$, which are described by Eqs. (20) and (21). The work done to create this warping is made by the stresses ${\sigma }_x$ and ${\tau }_{xy}$, which are associated with the displacements $U_T(L,y)$ and $V_T(L,y)$, respectively. The stresses are given by Eq. (1) and Eq. (2).

Due to the symmetric profile of stress and displacement with respect to the middle point M, in both directions the work along the entire `fixed’ section AMB is double of the work along the half-side MA. Therefore, considering that the subscript `b’ stands for bending whereas `s’ for shear, we have:

$$\begin{array}{cc}\hfill {\left(W_{AMB}\right)}_{T,b}& =2W_{MA}\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}{\sigma }_x{U}_T(L,y)bdy\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}\left[\left(-{\displaystyle \frac{PL}{I}}y\right)\left({\displaystyle \frac{(2+\nu )P}{6EI}}{y}^3\right)\right]bdy\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{3P^2(2+\nu )L}{10Ebh}},\hfill \end{array}$$

(27)

and

$$\begin{array}{cc}\hfill {\left(W_{AMB}\right)}_{T,s}& =2{\left(W_{MA}\right)}_s\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}{\tau }_{xy}{V}_T(L,y)bdy\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}\left[-{\displaystyle \frac{P}{2I}}(h^2/4-y^2)\right]\left[{\displaystyle \frac{\nu PL}{2EI}}{y}^2\right]bdy\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{3P^2\nu L}{10Ebh}}.\hfill \end{array}$$

(28)

The algebraic sum of the two above work components gives the total work for the linear and parabolic stresses (given by Eqs. (1) and (2), respectively) to create the warping, thus we have:

$$\begin{array}{cc}\hfill {\left(W_{\mathrm{support}}\right)}_T& ={\left(W_{AMB}\right)}_{T,b}+{\left(W_{AMB}\right)}_{T,s}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{3P^{2}L(1+\nu )}{5Ebh}}.\hfill \end{array}$$

(29)

The consideration of the above work ${\left(W_{\mathrm{support}}\right)}_T$ leads to a change of the work done by the total external shear force P (of which an equivalent may be taken at the middle of the left edge) thus:

$${\displaystyle \frac{1}{2}}P{(\Delta V_o)}_T={\left(W_{\mathrm{support}}\right)}_T.$$

(30)

It is worthy to mention that if the above work ${\left(W_{\mathrm{support}}\right)}_T$ is substituted in the approximate Eq. (9), the produced result $V_T(0,0)$ becomes

$$\begin{array}{cc}\hfill V_T(0,0)& ={\left(V_o\right)}_{\mathrm{approx}}-{(\Delta V_o)}_T\hfill \\ \hfill & =\left({\displaystyle \frac{P}{E}}{\displaystyle \frac{L^3}{3I}}\right)+{\displaystyle \frac{PL}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right]+{\displaystyle \frac{3PL(1+\nu )}{5Ebh}}\hfill \\ \hfill & =\left({\displaystyle \frac{P}{E}}{\displaystyle \frac{L^3}{3I}}\right)+{\displaystyle \frac{3PL(1+\nu )}{Ebh}}.\hfill \end{array}$$

(31)

Remark: Comparing Eq. (31) with Eq. (19), we see that after the above correction of Roark’s formula [19] considering the warping work, we receive exactly the solution proposed by Timoshenko & Goodier [16].

4.2. Livesley’s Formula [18] (`L’)

In this subsection we consider the displacement field according to the kinematical assumptions made by Livesley [18], in which the warping along the edge AMB is given according to Eqs. (25) and (). In a similar way, due to the symmetric profile of stress ${\sigma }_x$ and displacement $U(L,y)$ with respect to the middle M, the first work is double of the work along the half-side MA in the x-direction, and due to Eqs. (1) and (2) as well as Eqs. (25) and (), we have (`b’=bending, `s’=shear):

$$\begin{array}{cc}\hfill {\left(W_{ABC}\right)}_{L,b}& =2{\left(W_{AB}\right)}_{L,b}\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}{\sigma }_{x}U(L,y)bdy\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}\left[\left(-{\displaystyle \frac{PL}{I}}y\right)\left(-{\displaystyle \frac{(2+\nu )P}{6EI}}y(h^2/4-y^2)\right)\right]bdy\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}·{\displaystyle \frac{P^2(2+\nu )Lb{h}^5}{720E{I}^2}}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}·{\displaystyle \frac{P^2(2+\nu )L}{5Ebh}}.\hfill \end{array}$$

(32)

The second work is also double of the work along the half-side AB in the y-direction, thus

$$\begin{array}{cc}\hfill {\left(W_{AMB}\right)}_{L,s}& =2{\left(W_{MA}\right)}_{L,s}\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}{\tau }_{xy}V(L,y)bdy\hfill \\ \hfill & =2·{\displaystyle \frac{1}{2}}{\int }_0^{h/2}\left[-{\displaystyle \frac{P}{2I}}(h^2/4-y^2)\right]\left[{\displaystyle \frac{\nu PL}{2EI}}{y}^2\right]bdy\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{P^2\nu Lb{h}^5}{480E{I}^2}}\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}·{\displaystyle \frac{3P^2\nu L}{10Ebh}}.\hfill \end{array}$$

(33)

As previously happened, the algebraic sum of the two above work components gives the total work for the linear and parabolic stresses (given by Eqs. (1) and (2), respectively) to create the warping, thus we have:

$$\begin{array}{cc}\hfill {\left(W_{\mathrm{support}}\right)}_L& ={\left(W_{AMB}\right)}_{L,b}+{\left(W_{AMB}\right)}_{L,s}\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}·{\displaystyle \frac{P^{2}L(4-\nu )}{10Ebh}}.\hfill \end{array}$$

(34)

The consideration of the above work ${\left(W_{\mathrm{support}}\right)}_L$ leads to a change of the work done by the total external shear force P (of which an equivalent sum may be taken at the middle of the left edge) thus:

$${\displaystyle \frac{1}{2}}P{(\Delta V_o)}_L={\left(W_{\mathrm{support}}\right)}_L.$$

(35)

In a similar manner with subSection 4.1, if ${\left(W_{\mathrm{support}}\right)}_L$ is substituted in the approximate Eq. (9), now the produced result $V_L(0,0)$ becomes

$$\begin{array}{cc}\hfill V_L(0,0)& ={\left(V_o\right)}_{\mathrm{approx}}-{(\Delta V_o)}_L\hfill \\ \hfill & =\left({\displaystyle \frac{P{L}^3}{3EI}}\right)+{\displaystyle \frac{(4+5\nu )PL}{2Ebh}}.\hfill \end{array}$$

(36)

As previously, comparing Eq. (36) with Eq. (), we see that after the above correction we receive exactly the solution proposed by Livesley [18].

Remark: Comparing Eq. (29) with Eq. (34) one may observe that the two works done at the warped support are different each other. Obviously, their difference $\Delta W$ determines the difference $\Delta V$ between the associated deflections at the free end, through the equality $\frac{1}{2}P\Delta V=\Delta W$. This fact comes from the application of Eq. (9) twice, the first time for the `T’-model and the second time for the `L’-model and then their subtraction by parts.

5. Deflection Bounds

From the above discussion, we have determined three possible formulas, as follows:

• Roark’s formula, given by Eq. (11).
• Timoshenko’s formula, given by Eq. (31).
• Livesley’s formula, given by Eq. (36).

The classification of the beam depends on the aspect ratio (length:height), which is defined as follows:

$$\lambda ={\displaystyle \frac{L}{h}}.$$

(37)

Introducing the common coefficient:

$$a={\displaystyle \frac{PL}{Ebh}}={\displaystyle \frac{\lambda P}{Eb}},$$

(38)

the above estimations are proportional to $\lambda$ and can be written as follows:

$$\begin{array}{cc}\hfill \mathrm{Roark}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(11\left)\right)\text{:}{(V_{\mathrm{s}},\mathrm{Roark})}_{\max }& ={\displaystyle \frac{12}{5}}(1+\nu )a\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill \mathrm{Timoshenko}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(31\left)\right)\text{:}{(V_{\mathrm{s}},\mathrm{Timoshenko})}_{\max }& =3(1+\nu )a\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \mathrm{Livesley}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(36\left)\right)\text{:}{(V_{\mathrm{s}},\mathrm{Livesley})}_{\max }& ={\displaystyle \frac{4+5\nu }{2}}a.\hfill \end{array}$$

(41)

Taking the first derivative of the above set with respect to $\lambda$, we obtain:

$$\begin{array}{cc}\hfill \mathrm{Roark}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(11\left)\right)\text{:}{\displaystyle \frac{d{(V_{\mathrm{s}},\mathrm{Roark})}_{\max }}{d\lambda }}& ={\displaystyle \frac{12}{5}}(1+\nu )\left({\displaystyle \frac{P}{Eb}}\right)\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \mathrm{Timoshenko}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(31\left)\right)\text{:}{\displaystyle \frac{d{(V_{\mathrm{s}},\mathrm{Timoshenko})}_{\max }}{d\lambda }}& =3(1+\nu )\left({\displaystyle \frac{P}{Eb}}\right)\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \mathrm{Livesley}\text{’}\mathrm{s}\mathrm{formula}(\mathrm{Eq}.(36\left)\right)\text{:}{\displaystyle \frac{d{(V_{\mathrm{s}},\mathrm{Livesley})}_{\max }}{d\lambda }}& ={\displaystyle \frac{4+5\nu }{2}}\left({\displaystyle \frac{P}{Eb}}\right).\hfill \end{array}$$

(44)

A first observation on Eqs. (39)-to-(41) is that for the same Poisson’s ratio $\nu$, for any two formulas out of the three there is a corresponding linear relationship between the relevant shear deflection $V_s$, which is independent of the aspect ratio $\lambda$. For example, omitting the subscript `max’, we have:

$${\displaystyle \frac{V_{\mathrm{s}},\mathrm{Roark}}{V_{\mathrm{s}},\mathrm{Timoshenko}}}={\displaystyle \frac{4}{5}},$$

(45)

and

$${\displaystyle \frac{V_{\mathrm{s}},\mathrm{Roark}}{V_{\mathrm{s}},\mathrm{Livesley}}}={\displaystyle \frac{24(1+\nu )}{5(4+5\nu )}}.$$

(46)

Considering that Poisson’s ratio varies within the interval $0.15\le \nu \le 0.35$ (from concrete to steel), the value predicted by Roark’s formula (Eq. (11) and Eq. (39)) falls between those obtained using Timoshenko’s and Livesley’s boundary conditions (Eqs.(40) and (41), respectively). Therefore, it is imperative to determine which of them provides the closest approximation to the converged finite element solution, that will be regarded as the most accurate reference.

The above observations are validated using a finite element (FEM) model of plane-stress elasticity (4-node bilinear rectangular elements), with $n_x=600$ uniform subdivisions in the longitudinal x-direction. The number of elements in the transverse y-direction, $n_y$, is chosen such that all elements are square. The beam length is kept at a standard value (e.g., $L=3$), and its height $h=L/\lambda$ is continuously regulated accordingly, while all other quantities equal to unit ($b=1,E=1,P=1$).

The arithmetic value of length L does not play any role at all, because –for prescribed values of load P, material properties ($E,\nu$), and thickness b– the maximum deflection depends only on the aspect ratio $\lambda$, according the following formulas:

$$V_{\mathrm{Bernoulli}}={\displaystyle \frac{P{L}^3}{3EI}}={\displaystyle \frac{4P}{Eb}}{\lambda }^3,$$

(47)

where the shear component, second term of Roark’s formula (Eq. (11)), becomes:

$$V_{\mathrm{s}}={\displaystyle \frac{PL}{\left(hb\right)G}}\left[{\displaystyle \frac{6}{5}}\right]={\displaystyle \frac{12P(1+\nu )}{5Eb}}\lambda .$$

(48)

For example, if Poisson’s ratio is taken as $\nu =0.30$, one may observe in Table 2 that for a broad spectrum of aspect ratios $\lambda =L/h$, the FEM-values, which are taken as a reference, lie between Livesley’s model (i.e., Eq. (36)) and Roark’s model (Eq. (11)). In contrast, Timoshenko’s formula (i.e., Eq. (31)) overestimates the maximum deflection. By subtracting the Bernoulli component, $V_{\mathrm{Bernoulli}}$, the same results (in terms of only the shear displacement $V_{\mathrm{s}}$) are shown in Figure 3.

Overall, for a broader range of Poisson’s ratios ($\nu =0.15,0.20,0.25,0.30$), the shear deflection $V_s$, coming from the FEM solution, is illustrated in Figure 4.

One may observe:

• With adequate accuracy, the FEM solution can be represented by a straight line in the graph $V_s$ versus $\lambda$. The slope of this line depends on Poisson’s ratio $\nu$, as also shown in Figure 5. One may observe that as the Poisson’s ratio increases the shear deflection increases as well, and thus the corresponding line moves upwards.
• For the smallest studied aspect ratio $\lambda =1$, and any value of Poisson’s ratio $\nu$, the FEM-solution almost belongs to the Roark’s curve, defined by Eq. (11).
• As Poisson’s ratio $\nu$ increases, the FEM-curve tends to Livesley’s curve, that is defined by Eq. (36).
• For Poisson’s ratio $\nu =0.30$ and the largest studied aspect ratio $\lambda =5$, the FEM-solution almost belongs to the Livesley’s curve, defined by Eq. (36).

6. Curve Fitting

In the previous section, we shaw that for each of the four Poisson’s ratios $\nu \in [0.15,0.30]$ we have determined one polygonal line which consists of nine points. Each of these nine points, is associated with one discrete aspect ratio, which varies in the interval $\lambda \in [1.0,5.0]$, with a constant step $\Delta \lambda =0.5$. All the four curves (one for each Poisson’s ratio) can be fitted in several manners, but in this paper we present only two of them, as follows:

• Linear interpolation. The regression formulas and the R-squared values are shown in Table 3.
• Cubic interpolation (true cubic polynomial interpolation), with $R^2=1.0000000$ is shown in Table 4.

7. Interpolation Algorithms

Since the FEM solution was found to be bounded by Livesley’s and Roark’s curves, and because it monotonically rotates in the clockwise direction, for a given arbitrary pair $(\nu ,\lambda$) we can approximate the FEM solution implementng one of the following algorithms (Section 7.1-to-Section 7.3).

7.1. Algorithm-1: Using only the End-Points and Linear Interpolations

This algorithm is as follows:

• We perform linear interpolation between the four curves at $\lambda =1$, and thus determine the value ${\left(V_s\right)}_{\mathrm{left}}$ (denoted by `Left’ on the blue line of Figure 6).
• We perform linear interpolation between the four curves at $\lambda =5$, and thus determine the value ${\left(V_s\right)}_{\mathrm{right}}$ (denoted by `Right’ on the red line of Figure 6).
• Based on the above two values, ${\left(V_s\right)}_{\mathrm{left}}$ and ${\left(V_s\right)}_{\mathrm{right}}$, we linearly interpolate in terms of $\lambda$, and thus determine the desired value $V_s$.

In general, the linear interpolation between ${\left(V_s\right)}_{\mathrm{left}}$ and ${\left(V_s\right)}_{\mathrm{right}}$ leads to a small error. For example, if we pretend that the value $V_s$ at ($\nu =0.20$ and $\lambda =3.0$) is unknown, the above procedure gives ${\left(V_s\right)}_{\mathrm{interp}}=8.0958$, while the FEM solution is ${\left(V_s\right)}_{\mathrm{FEM}}=8.2245$. Therefore, a more accurate procedure is under demand (see, Section 7.2-to-Section 7.3).

7.2. Algorithm-2: Linear Interpolation Based on Local Data

An alternative procedure is to work as follows:

• For the given $\lambda$-value, we perform linear interpolation and determine four values $V_s$ on the corresponding polygonal curves of the FEM-results associated with four Poisson’s ratios $\nu =0.15,0.20,0.25,0.30$.
• For the given $\nu$-value, which is generally different than the above standardized ones, we perform a linear interpolation among the above four values $V_s$.

A MATLAB^® code, which implements Algorithm-2, is given in Appendix A. The reader may execute this code and will find the anticipated value $V_s=8.2245$.

7.3. Algorithm-3: Higher Degree Interpolation

This is probably the most accurate alogithm of this Section. A similar procedure with Algorithm-2 is applied, with the difference that –instead of the abovementioned linear interpolation (polygonal line by FEM results)– we perform cubic interpolation as follows:

• For the given $\lambda$-value, we apply a cubic interpolation –according to Table 4– on each of the four cubic lines associated with four Poisson’s ratios $\nu =0.15,0.20,0.25,0.30$. Therefore, four discrete value $V_s$ are determined.
• Having found the above four values on the interpolated curves, we then perform linear (or better cubic interpolation) between these four values.

The application of Algorithm-3 leads to the anticipated value $V_s=8.2245$.

8. Discussion

The mechanical behaviour of homogeneous beams mainly depends on the aspect ratio $\lambda =L/h$, where L is the length and h the height. The standard definitions typically break down as follows:

• $\lambda \le 2$: Deep beam — Requires 2D or 3D theory (not simple beam theory). Shear deformation is dominant; plane sections do not remain plane. 2D elasticity needed to capture full stress state.
• $2\lesssim \lambda \lesssim 5\text{–}10$: Timoshenko beam — First-order shear deformation theory. Both bending and shear are significant. Timoshenko or 2D plane stress gives good accuracy.
• $\lambda \gtrsim 10$: Euler–Bernoulli beam — Classical beam theory applies. Shear deformation is negligible; plane sections remain plane.

This study focused on the case of a clamped beam subjected to a transverse tip load, with an aspect ratio ranging from $1\le L/h\le 5$. Over this entire range, the plane-stress model is applicable. We found that the calculated FEM curve (maximum deflection versus aspect ratio) exhibits near-linear behavior and is bounded by the values proposed by Livesley [18] and Roark’s formula [19], which serve as lower and upper limits, respectively. Conversely, a formula cited in the textbook by Timoshenko and Goodier [16] yields higher, thus conservative, values.

For discrete Poisson’s ratios of $\nu =$ 0.15, 0.20, 0.25, and 0.30, we performed linear and cubic regression analyses on the finite element method (FEM) results. These regressions establish a relationship between the maximum shear deflection ($V_s$) and the continuous variable aspect ratio ( $\lambda$). By incorporating the conventional Bernoulli bending term ( ${\displaystyle \frac{P{L}^3}{3EI}}$), the total deflection can be readily determined.

Based on the four closed-form expressions, we propose the following approach for an arbitrary Poisson’s ratio $\overline{\nu }$: For a given aspect ratio $\lambda =L/h$, we utilize readily available cubic fitting formulas to determine the maximum shear deflection $V_s$ for discrete Poisson’s ratio values of $\nu =$0.15,0.20,25,30, and 0.30. Subsequently, we merely interpolate for the specified Poisson’s ratio $\overline{\nu }$.

The case of very small aspect ratios ($\lambda <1$) was not addressed in this paper. In this range, the plane-stress model must be replaced by a plane-strain model, or a 3D-modeling approach must be employed.

The inability of the three analytical solutions (proposed by Timoshenko, Roark, and Livesley, respectively) to accurately represent reality, as better depicted by the numerical FEM solution, stems from the boundary layer developed at the clamped ends. This boundary layer, characterized by rapid stress variation and traction discontinuity, cannot be precisely captured by the smooth polynomial expressions inherent in these analytical formulas.

While the Introduction addresses both displacement and stress analysis, this paper focuses solely on the former. Ongoing investigations into curve-fitting the normal stress along the clamp indicate that polynomial fitting is unsuitable. Instead, a high number of harmonics effectively captures its variation. Further research on stress analysis and the associated shear-lag effect [20,21,22,23,24] will be presented in a future publication.

9. Conclusions

This study demonstrates that the Finite Element Method (FEM) solution for a clamped beam subjected to a tip load is effectively bounded by three established closed-form analytical solutions: one providing a lower bound and two serving as upper bounds. However, a single correction factor cannot reconcile these analytical expressions with the FEM results due to the differing slopes of their respective linear regressions. To address this, we proposed fitted analytical expressions for four discrete Poisson’s ratios, enabling estimation of maximum deflection as a function of aspect ratio within the range of 1 to 5. Additionally, we provided an algorithm and its source code for the general case of arbitrary Poisson’s ratios ($\nu \in [0.15,0.30]$) and aspect ratios ($\lambda \in [1,5]$), offering a practical tool for comprehensive deflection analysis.