A direct asymptotic integration of the full three-dimensional problem of elasticity is employed to derive a consistent governing equation for a beam with the rectangular cross section. The governing equation is consistent in the sense that it has the same long-wave low-frequency behaviour as the exact solution of the original three-dimensional problem. Performance of the new beam equation is illustrated by comparing its predictions against the results of direct finite-element computations. Limiting behaviours for beams with large (and small) aspect ratios, which can be established using classical plate theories, are recovered from the new governing equation to illustrate its consistency and also to illustrate the importance of using plate theories with the correctly refined boundary conditions. The implications for the correct choice of the shear correction factor in Timoshenko's beam theory are also discussed.

1. Introduction

The classical Euler–Bernoulli equation can be interpreted in two ways. On the one hand, it can be understood as an approximate equation that broadly reproduces behaviour of thin beams at sufficiently low frequencies. This interpretation is valid and, in fact, this is how the Euler–Bernoulli equation is most commonly introduced, typically justified by ad hoc strength of materials style arguments. On the other hand, another, more mathematically precise definition is also possible. One can start with a full three-dimensional boundary value problem for a prism composed of an elastic material and consider a long-wave low-frequency asymptotic expansion under the assumptions that frequencies of interest are low and characteristic dimensions of the cross section are negligible in relation to the wavelength of interest. It can be shown that the Euler–Bernoulli equation recovers the leading-order term of this expansion. In this sense, the classical theory is precisely defined and fully mathematically consistent. This also suggests that further terms in the same asymptotic expansion can be considered to, hopefully, construct more precise, yet still asymptotically consistent approximations to the full three-dimensional solution. This is precisely the topic of our discussion.

It may seem surprising that higher-order equations for a beam can still warrant an investigation. After all, more geometrically complex objects, such as plates and shells, have well-developed higher-order theories (e.g. [1–4]). The explanation is that cross-sectional problems for both plates and shells are always one-dimensional and typically possess simple analytical solutions, whereas cross-sectional problems for beams are necessarily two-dimensional, which makes their solution substantially harder. It is also worth noting that deriving new higher-order equations is not the goal in itself; instead, one can use them to construct higher-order asymptotic solutions for much harder boundary value problems, e.g. vibration of elliptic plates [5] or scattering of acoustic waves by an immersed shell [6] or solution of initial-value problems [7].

The most popular refined beam equation was obtained by Timoshenko [8] (see recent review [9]). It is derived using several ad hoc assumptions and, as a result, does not generally predict the correct higher-order asymptotic behaviour of the beam. It features two vibration modes, a lower (fundamental) mode that describes bending, as well as a higher mode that bears only superficial similarity to the true three-dimensional solution, see [10]. Nevertheless, predictions by Timoshenko's governing equation can be made partly consistent because the equation features an imprecisely defined constant κ, usually called the shear correction factor. It turns out that κ can be chosen in such way that the fundamental mode of Timoshenko's theory matches the appropriate expansion of the corresponding exact three-dimensional solution, see [11].

Berdichevskii & Kvashnina [12] were apparently the first who attempted to construct a more consistent equivalent of Timoshenko's equation using the variational asymptotic method. Unfortunately, when specialized to the case of rectangular beams, their equation appears to contain an erroneous coefficient, this will be further discussed in §7. The methodology of Berdichevskii and Kvashnina was generalized in [13–15], focusing mainly on derivations of the approximations of the potential energy.

The goal of this paper is to present the derivation of a truly asymptotically consistent higher-order theory for rectangular beams. Our derivation is based upon the direct asymptotic integration of the exact three-dimensional problem of elasticity. This method has been originally developed by Goldenveizer, see e.g. [16], as well as more recent paper [3]. Such an approach allows one to construct internally consistent higher-order approximations for the stress and strain fields across the cross section. This makes our approach largely similar to [17], in which a more general problem for a beam composed of an anisotropic elastic material has been considered. The main difference is that the coefficients of beam equation in [17] were obtained by solving a cross section problem numerically, whereas we construct a governing equation with fully explicit analytical expressions for all coefficients.

The restriction to the rectangular cross section enables us to obtain a fully analytical description of the simplest non-trivial configuration of the beam, the classical benchmark problem that we believe is still not resolved in the literature conclusively. We validate new beam equation by showing that the dispersion curves by our theory faithfully reproduce the long-wave low-frequency behaviour of the exact dispersion curves computed using the SAFE method for several aspect ratios of the beam. The choice of rectangular cross section also allows us to provide an analytical validation of the results. Specifically, motion of rectangular beams with large or small aspect ratios can be described by the classical plate theories: the out-of-plane bending mode can be reproduced as symmetric fundamental mode of a Kirchhoff plate strip and the in-plane bending mode can be reproduced by studying antisymmetric fundamental mode of a plate strip in the state of plane stress. The Kirchhoff plate strip limit is particularly interesting, because we show that the use of refined boundary conditions for the plate theory by Kolos [18] allows to reproduce the correct asymptotic behaviour of even a moderately thick flat beam.

It is instructive to compare the predictions by our equation to the response of Timoshenko's beam equation. We do this by matching the relevant dispersion curves; this allows us to propose a new expression for the shear correction factor κ. In his papers, Timoshenko used constant (cross-section-independent) values of κ, see [8,19]; the most commonly cited reference [20] provides values of κ for a number of canonical cross sections that also depend on the Poisson ratio ν. However, the coefficients given in [20] are known to be asymptotically inconsistent, at least for circular and rectangular cross sections, see [11]. A particular expression of κ, first obtained in [21], was independently confirmed in [22–25]. The expression for value of κ for rectangular beams obtained in this paper confirms the consistency of the coefficient proposed in [21]. In addition, we demonstrate that a number of other definitions of shear correction factor given in the recent literature do not lead to a consistent form of Timoshenko's equation.

2. Governing equations

Consider an infinite rectangular beam composed of a homogeneous isotropic elastic material. Motions of isotropic media are governed by the equations of three-dimensional elasticity

τ_{i j, j} = ρ {\ddot{u}}_{i}, τ_{i j} = \frac{E}{1 + ν} (\frac{ν ϵ_{k k} δ_{i j}}{1 - 2 ν} + ϵ_{i j}), ϵ_{i j} = \frac{1}{2} (u_{i, j} + u_{j, i}),

2.1

where

τ_{i j} = τ_{i j} (x, t)

is the stress tensor,

u_{i} = u_{i} (x, t)

the displacement vector, ρ the mass density per unit volume, E the Young modulus and ν the Poisson ratio. Comma subscripts in (2.1) indicate differentiation with respect to implied spatial coordinates; henceforth we will also use comma subscripts to denote differentiation with respect to explicitly shown non-dimensional spatial coordinates and time. Also assumed specifically in (2.1) are the summation over repeated suffices and overdots denoting differentiation with respect to time t.

We place the origin of the Cartesian coordinate system at the symmetry axis of the beam and direct axis Ox₁ along the symmetry axis. The beam cross section is the rectangle {(x₂,x₃):−h₂≤x₂≤h₂,−h₃≤x₃≤h₃} (figure 1). The boundary conditions on the free faces of the beam may then be written as

\begin{aligned} τ_{2 i} & = 0 (i = 1, 2, 3) at x_{2} = \pm h_{2} \\ and τ_{3 i} & = 0 (i = 1, 2, 3) at x_{3} = \pm h_{3} . \end{aligned}

2.2

Figure 1. The dimensional and non-dimensional coordinate systems for the considered problem.

Our primary interest is with the propagation of long waves, i.e. waves whose wavelength l is much greater than a characteristic dimension of the cross section h.¹ This suggests introducing a natural small parameter ε=h/l and re-scaling the problem in terms of the following non-dimensional coordinates:

x_{1} = l ξ, x_{2} = ε l η, x_{3} = ε l ζ, t = \frac{ε^{- 1} l τ}{\sqrt{E / ρ}} .

2.3

An additional assumption that the beam bends along the axis Ox₃ prompts an appropriate re-scaling for the displacement and stress components:

\begin{aligned} u_{1} & = ε l u, u_{2} = ε^{2} l v, u_{3} = l w, \\ τ_{11} & = ε E σ_{11}, τ_{33} = ε^{3} E σ_{33}, \\ τ_{1 i} & = ε^{2} E σ_{1 i}, τ_{2 i} = ε^{3} E σ_{2 i} (i = 2, 3) . \end{aligned}

2.4

The introduction of scalings (2.4) can be partly motivated by referring to the classical Euler–Bernoulli theory and also by the asymptotic analysis of the exact dispersion relations for similar problems (e.g. for the bending of a circular beam). Here, for the sake of brevity, we introduce the ansatz constructively and prove its consistency a posteriori, by observing that the proposed form of the asymptotic series satisfies appropriate equations and boundary conditions at all orders.

In terms of the newly introduced non-dimensional quantities, we can re-write governing equations (2.1) in the following form:

\begin{aligned} σ_{11, ξ} + σ_{12, η} + σ_{13, ζ} & = ε^{2} u_{, τ τ}, \\ σ_{21, ξ} + σ_{22, η} + σ_{23, ζ} & = ε^{2} v_{, τ τ}, \\ σ_{31, ξ} + σ_{32, η} + σ_{33, ζ} & = w_{, τ τ}, \end{aligned}

2.5

and

\begin{aligned} ε^{2} (1 + ν) (1 - 2 ν) σ_{11} & = ε^{2} (1 - ν) u_{, ξ} + ε^{2} ν v_{, η} + ν w_{, ζ}, \\ ε^{4} (1 + ν) (1 - 2 ν) σ_{22} & = ε^{2} ν u_{, ξ} + ε^{2} (1 - ν) v_{, η} + ν w_{, ζ}, \\ ε^{4} (1 + ν) (1 - 2 ν) σ_{33} & = ε^{2} ν u_{, ξ} + ε^{2} ν v_{, η} + (1 - ν) w_{, ζ}, \\ 2 ε^{2} (1 + ν) σ_{12} & = u_{, η} + ε^{2} v_{, ξ}, \\ 2 ε^{2} (1 + ν) σ_{13} & = u_{, ζ} + w_{, ξ}, \\ and 2 ε^{4} (1 + ν) σ_{23} & = ε^{2} v_{, ζ} + w_{, η} . \end{aligned}

2.6

Boundary conditions (2.2) can also be appropriately re-formulated as

\begin{aligned} σ_{2 i} & = 0 (i = 1, 2, 3) at η = \pm a \\ and σ_{3 i} & = 0 (i = 1, 2, 3) at ζ = \pm b, \end{aligned}

2.7

where a=h₂/h and b=h₃/h.

The solution f={u,v,w,σ_ij} of the boundary value problem (2.5)–(2.7) is now sought in form of the following asymptotic ansatz:

f (ξ, η, ζ, τ) = f^{(0)} + ε^{2} f^{(2)} + ε^{4} f^{(4)} + \dots,

2.8

where f^(k)=f^(k)(ξ,η,ζ,τ) and

f^{(k)} (ξ, η, ζ, τ) = \sum_{m} {\tilde{f}}_{m}^{(k)} (ξ, τ) F_{m}^{(k)} (η, ζ) .

2.9

3. The leading-order theory

The leading-order displacements are governed by the leading-order terms of governing equations (2.6). The relevant equations can be obtained by substituting expansion (2.8) into equations (2.6) and neglecting all terms that are of O(ε²) (and smaller), which yields

w_{, ζ}^{(0)} = 0, u_{, η}^{(0)} = 0, u_{, ζ}^{(0)} + w_{, ξ}^{(0)} = 0, w_{, η}^{(0)} = 0.

3.1

The appropriate solution is given by

w^{(0)} = {\tilde{w}}_{0}, u^{(0)} = - {\tilde{w}}_{0}^{'} ζ,

3.2

where

{\tilde{w}}_{0} = {\tilde{w}}_{0} (ξ, τ)

and a dash denotes differentiation with respect to ξ.

The governing equations for the next-order displacements and remaining leading-order quantities have the form

\begin{aligned} σ_{11, ξ}^{(0)} + σ_{12, η}^{(0)} + σ_{13, ζ}^{(0)} & = 0, \end{aligned}

3.3

\begin{aligned} σ_{21, ξ}^{(0)} + σ_{22, η}^{(0)} + σ_{23, ζ}^{(0)} & = 0, \end{aligned}

3.4

\begin{aligned} σ_{31, ξ}^{(0)} + σ_{32, η}^{(0)} + σ_{33, ζ}^{(0)} & = {\tilde{w}}_{, τ τ}, \end{aligned}

3.5

and

\begin{aligned} (1 + ν) (1 - 2 ν) σ_{11}^{(0)} & = (1 - ν) u_{, ξ}^{(0)} + ν v_{, η}^{(0)} + ν w_{, ζ}^{(2)}, \end{aligned}

3.6

\begin{aligned} 0 & = ν u_{, ξ}^{(0)} + (1 - ν) v_{, η}^{(0)} + ν w_{, ζ}^{(2)}, \end{aligned}

3.7

\begin{aligned} 0 & = ν u_{, ξ}^{(0)} + ν v_{, η}^{(0)} + (1 - ν) w_{, ζ}^{(2)}, \end{aligned}

3.8

\begin{aligned} 2 (1 + ν) σ_{12}^{(0)} & = u_{, η}^{(2)} + v_{, ξ}^{(0)}, \end{aligned}

3.9

\begin{aligned} 2 (1 + ν) σ_{13}^{(0)} & = u_{, ζ}^{(2)} + w_{, ξ}^{(2)}, \end{aligned}

3.10

\begin{aligned} 0 & = v_{, ζ}^{(0)} + w_{, η}^{(2)} . \end{aligned}

3.11

Using the leading-order solution (3.2) and taking into account the (anti-) symmetry of the problem allow one to integrate governing equations (3.7), (3.8) and (3.11) with the following result:

v^{(0)} = ν {\tilde{w}}_{0}^{″} η ζ, w^{(2)} = \frac{1}{2} ν {\tilde{w}}_{0}^{″} (ζ^{2} - η^{2}) + {\tilde{w}}_{2},

3.12

where

{\tilde{w}}_{2} = {\tilde{w}}_{2} (ξ, τ)

. Governing equation (3.6) then gives

σ_{11}^{(0)} = - {\tilde{w}}_{0}^{″} ζ .

3.13

Function ${\tilde{w}}_{0}$ underlying the leading-order solution still has not been specified. It can be determined by multiplying equation (3.3) by ζ and integrating over the cross section represented by rectangular region R, where

R = {(η, ζ) | - a \leq η \leq a, - b \leq ζ \leq b} .

3.14

In the view of the boundary conditions at η=±a and ζ=±b, this yields

\iint_{R} \frac{\partial σ_{13}^{(0)}}{\partial ξ} d η d ζ = - \frac{4}{3} a b^{3} {\tilde{w}}_{0}^{‴} .

3.15

Similar integration of (3.5) over R gives then the equation for

{\tilde{w}}_{0}

\frac{1}{3} b^{2} \frac{\partial^{4} {\tilde{w}}_{0}}{\partial ξ^{4}} + \frac{\partial^{2} {\tilde{w}}_{0}}{\partial τ^{2}} = 0,

3.16

which can be recognized as the classical Euler–Bernoulli equation for a rectangular beam. Indeed, given that

{\tilde{w}}_{0} \sim w = u_{3} / l

, we can formulate the dimensional equivalent of equation (3.16) as

\frac{E h_{3}^{2}}{3 ρ} \frac{\partial^{4} u_{3}}{\partial x_{1}^{4}} + \frac{\partial^{2} u_{3}}{\partial t^{2}} = 0 .

3.17

Keeping in mind that the second moment for rectangular cross section I_x₂=2h₂(2h₃)³/12 and the cross-sectional area A=4h₂h₃, the coefficient in front of the fourth-order term can now be recognized as the usual constant EI_x₂/ρA.

To obtain an expression for u⁽²⁾, we first express $σ_{12}^{(0)}$ and $σ_{13}^{(0)}$ from equations (3.9) and (3.10), and then substitute the resulting expressions into (3.3) and boundary conditions $σ_{12}^{(0)} |_{η = \pm a} = 0$ and $σ_{13}^{(0)} |_{ζ = \pm b} = 0$ . As a result of these manipulations, we are led to the following boundary value problem:

\begin{matrix} u_{, η η}^{(2)} + u_{, ζ ζ}^{(2)} = 2 {\tilde{w}}_{0}^{‴} ζ, \\ u_{, η}^{(2)} |_{η = \pm a} = \mp ν {\tilde{w}}_{0}^{‴} a ζ, u_{, ζ}^{(2)} |_{ζ = \pm b} = \frac{1}{2} ν {\tilde{w}}_{0}^{‴} (η^{2} - b^{2}) - {\tilde{w}}_{2}^{'} . \end{matrix}

3.18

We are seeking the solution of problem (3.18) in the form

u^{(2)} = {\tilde{w}}_{0}^{‴} U^{(2)} (η, ζ) - {\tilde{w}}_{2}^{'} ζ .

3.19

It is easy to verify that if one takes

U^{(2)} = \frac{1}{6} (2 - ν) ζ^{3} - b^{2} ζ + \frac{1}{2} ν ζ η^{2} + Ψ (η, ζ),

3.20

then the problem (3.18) is reduced to the boundary value problem

\begin{matrix} Ψ_{, η η} + Ψ_{, ζ ζ} = 0, \\ Ψ_{, η} |_{η = \pm a} = \mp 2 ν a ζ, Ψ_{, ζ} |_{ζ = \pm b} = 0 . \end{matrix}

3.21

The solution of (3.21) is given by

Ψ = \frac{32 ν a b^{2}}{π^{3}} \sum_{n = 1}^{\infty} C_{n} \sin (\frac{(2 n - 1) π}{2 b} ζ) \cosh (\frac{(2 n - 1) π}{2 b} η),

3.22

within which

C_{n} = \frac{{(- 1)}^{n}}{{(2 n - 1)}^{3}} csch (\frac{(2 n - 1) π a}{2 b}) .

3.23

As a result, we have

σ_{12}^{(0)} = \frac{{\tilde{w}}_{0}^{‴}}{2 (1 + ν)} [Ψ_{, η} + 2 ν η ζ] and σ_{13}^{(0)} = \frac{{\tilde{w}}_{0}^{‴}}{2 (1 + ν)} [Ψ_{, ζ} + ζ^{2} - b^{2}] .

3.24

It can be readily verified that substitution of (3.24)₂ into (3.5) and integration over the cross section R once again gives us equation (3.16) that governs

{\tilde{w}}_{0}

. Therefore, we have now satisfied all of the O(ε²) equations (3.3)–(3.11). To obtain a non-trivial correction term(s) for the governing equation (3.16), we now have to proceed to the next asymptotic order.

4. The higher-order theory

The consideration of O(ε⁴) terms after substituting asymptotic ansatz (2.8) into equations (2.5) and (2.6) results in the following system of governing equations:

\begin{aligned} σ_{11, ξ}^{(2)} + σ_{12, η}^{(2)} + σ_{13, ζ}^{(2)} = u_{, τ τ}^{(0)}, \end{aligned}

4.1

\begin{aligned} σ_{31, ξ}^{(2)} + σ_{32, η}^{(2)} + σ_{33, ζ}^{(2)} = w_{, τ τ}^{(2)} \end{aligned}

4.2

and

\begin{aligned} (1 + ν) (1 - 2 ν) σ_{11}^{(2)} & = (1 - ν) u_{, ξ}^{(2)} + ν v_{, η}^{(2)} + ν w_{, ζ}^{(4)}, \end{aligned}

4.3

\begin{aligned} (1 + ν) (1 - 2 ν) σ_{22}^{(0)} & = ν u_{, ξ}^{(2)} + (1 - ν) v_{, η}^{(2)} + ν w_{, ζ}^{(4)}, \end{aligned}

4.4

\begin{aligned} (1 + ν) (1 - 2 ν) σ_{33}^{(0)} & = ν u_{, ξ}^{(2)} + ν v_{, η}^{(2)} + (1 - ν) w_{, ζ}^{(4)} \end{aligned}

4.5

(note that we omitted here the equations that will not be required for deriving the refined governing equation for

\tilde{w}

We begin solving system (4.1)–(4.5) by integrating equation (4.2) over the cross section R. If conditions on free faces are taken into account, we obtain

\iint_{R} σ_{13, ξ}^{(2)} d η d ζ = \iint_{R} w_{, τ τ}^{(2)} d η d ζ .

4.6

As the next step, we multiply equation (4.1) by ζ and integrate over the cross section to determine that

\iint_{R} σ_{13}^{(2)} d η d ζ = \iint_{R} ζ (σ_{11, ξ}^{(2)} - u_{, τ τ}^{(0)}) d η d ζ .

4.7

The stress component

σ_{11}^{(2)}

on the right-hand side can be expressed using equations (4.3)–(4.5) in the form

σ_{11}^{(2)} = ν (σ_{22}^{(0)} + σ_{33}^{(0)}) + u_{, ξ}^{(2)},

4.8

and substituted back into the right-hand side of (4.7). In addition, by multiplying equations (3.4), (3.5) and (3.5) by ηζ,

\frac{1}{2} η^{2}

and

\frac{1}{2} ζ^{2}

, respectively, and integrating the obtained equations over the cross-section, we obtain the following identities:

\begin{aligned} \iint_{R} ζ σ_{22}^{(0)} d η d ζ & = \iint_{R} η ζ σ_{12, ξ}^{(0)} d η d ζ - \iint_{R} η σ_{23}^{(0)} d η d ζ, \end{aligned}

4.9

\begin{aligned} \iint_{R} η σ_{23}^{(0)} d η d ζ & = - \iint_{R} \frac{1}{2} η^{2} ({\tilde{w}}_{0, τ τ} - σ_{13, ξ}^{(0)}) d η d ζ \end{aligned}

4.10

\begin{aligned} and \iint_{R} ζ σ_{33}^{(0)} d η d ζ & = - \iint_{R} \frac{1}{2} ζ^{2} ({\tilde{w}}_{0, τ τ} - σ_{13, ξ}^{(0)}) d η d ζ . \end{aligned}

4.11

By combining these results together, one arrives at the equality

\begin{aligned} \iint_{R} ζ σ_{11}^{(2)} d η d ζ = \iint_{R} {ν [η ζ σ_{12, ξ}^{(0)} + \frac{1}{2} (η^{2} - ζ^{2}) ({\tilde{w}}_{0, τ τ} - σ_{13, ξ}^{(0)})] + ζ u_{, ξ}^{(2)}} d η d ζ . \end{aligned}

4.12

In view of the specific structure of expressions for

σ_{i j}^{(k)}

, see (2.8) and (2.9), we differentiate equation (4.12) with respect to ξ and substitute the result into the right-hand side of equation (4.7). Then we differentiate equation (4.7) with respect to ξ and substitute the result into the left-hand side of (4.6). The resulting expression for (4.6) is given by

\iint_{R} {\frac{\partial^{6} {\tilde{w}}_{0}}{\partial ξ^{6}} Φ (η, ζ) + \frac{\partial^{4} {\tilde{w}}_{0}}{\partial ξ^{2} \partial τ^{2}} [ν η^{2} + (1 - ν) ζ^{2}] - \frac{\partial^{4} {\tilde{w}}_{2}}{\partial ξ^{4}} ζ^{2} - \frac{\partial^{2} {\tilde{w}}_{2}}{\partial τ^{2}}} d η d ζ = 0,

4.13

in which

\begin{aligned} Φ (η, ζ) & = ζ Ψ + \frac{ν}{4 (1 + ν)} {2 η ζ \frac{\partial Ψ}{\partial η} + (ζ^{2} - η^{2}) \frac{\partial Ψ}{\partial ζ} \\ + (6 ν + 1) η^{2} ζ^{2} + \frac{(4 + 5 ν - 2 ν^{2})}{3 ν} ζ^{4} - \frac{(4 + 5 ν)}{ν} b^{2} ζ^{2} + b^{2} η^{2}} . \end{aligned}

4.14

The integral within equation (4.13) can be computed explicitly, yielding the governing equation for

{\tilde{w}}_{2}

in the form

\frac{b^{2}}{3} \frac{\partial^{4} {\tilde{w}}_{2}}{\partial ξ^{4}} + \frac{\partial^{2} {\tilde{w}}_{2}}{\partial τ^{2}} + k_{22} \frac{\partial^{4} {\tilde{w}}_{0}}{\partial ξ^{2} \partial τ^{2}} + k_{60} \frac{\partial^{6} {\tilde{w}}_{0}}{\partial ξ^{6}} = 0,

4.15

with the coefficients

\begin{aligned} k_{22} & = - \frac{b^{2}}{3} [1 - ν + ν q^{2}] \\ and k_{60} & = \frac{b^{4}}{45 (1 + ν)} [12 + 27 ν + 3 ν^{2} - 5 q^{2} ν (1 + 3 ν) + 15 ν^{2} q S (q)], \end{aligned}

4.16

defined using

q = \frac{a}{b} \equiv \frac{h_{2}}{h_{3}} and S (q) = \frac{384}{π^{5}} \sum_{n = 1}^{\infty} \frac{\coth (q π (2 n - 1) / 2)}{{(2 n - 1)}^{5}} .

4.17

Keeping in mind the form of the original asymptotic expansion (2.8) and, specifically, the fact that $\tilde{w} = {\tilde{w}}_{0} + ε^{2} {\tilde{w}}_{2} + O (ε^{4})$ , we can combine equations (3.16) and (4.15) into

\frac{b^{2}}{3} \frac{\partial^{4} \tilde{w}}{\partial ξ^{4}} + \frac{\partial^{2} \tilde{w}}{\partial τ^{2}} + ε^{2} [k_{22} \frac{\partial^{4} \tilde{w}}{\partial ξ^{2} \partial τ^{2}} + k_{60} \frac{\partial^{6} \tilde{w}}{\partial ξ^{6}}] = 0,

4.18

which can be re-written in an asymptotically equivalent form

\frac{b^{2}}{3} \frac{\partial^{4} \tilde{w}}{\partial ξ^{4}} + \frac{\partial^{2} \tilde{w}}{\partial τ^{2}} - ε^{2} k_{0} (ν, q) \frac{b^{2}}{3} \frac{\partial^{4} \tilde{w}}{\partial ξ^{2} \partial τ^{2}} = 0,

4.19

where

k_{0} (ν, q) = \frac{1}{5 (1 + ν)} [17 + 27 ν - 2 ν^{2} - 10 ν^{2} q^{2} + 15 ν^{2} q S (q)] .

4.20

Equation (4.19) can also be re-scaled back to the dimensional form, yielding

\frac{E h_{3}^{2}}{3 ρ} \frac{\partial^{4} u_{3}}{\partial x_{1}^{4}} + \frac{\partial^{2} u_{3}}{\partial t^{2}} - \frac{1}{3} h_{3}^{2} k_{0} (ν, q) \frac{\partial^{4} u_{3}}{\partial x_{1}^{2} \partial t^{2}} = 0 .

4.21

This equation, together with the expression for a non-dimensional coefficient k₀(ν,q), constitute the main result of this paper. We will now turn to the assessment of the qualitative and quantitative behaviour of the newly derived beam equation.

5. Numerical examples

Before we discuss the numerical performance of the beam equation (4.21), it is worth making few remarks about practical computation of the function S(q) used to define k₀(ν,q), see (4.17)₂. First, given a fixed q, the infinite series within (4.17)₂ is absolutely convergent and converges very rapidly. The first two terms of the series already provide relative errors below 0.1%; the first three terms result in relative errors below 0.01%. An asymptotic expansion for q≪1 has the form

S (q) = \frac{4}{5 q} + \frac{2}{3} q - \frac{2}{15} q^{3} + \dots .

5.1

A useful representation of (4.17)₂ for q≫1 can be formulated by using the fact that

\coth x = 1 + 2 / (e^{2 x} - 1)

; the results outlined in appendix A, especially equation (A.4), indicate that

S (q) = \frac{372}{π^{5}} ζ (5) + \frac{768}{π^{5}} \sum_{n = 1}^{\infty} \frac{1}{{(2 n - 1)}^{5}} \frac{1}{e^{π q (2 n - 1)} - 1},

5.2

where ζ(x) is the Riemann zeta function and ζ(5)≈1.0369278.

We can now turn our attention to assessing how well the new equation reproduces the dynamics of rectangular beams. This can be done by comparing the dispersion curves of the beam modelled by the full three-dimensional theory with the dispersion curves predicted by the one-dimensional model. The dispersion is characterized by looking at plane wave solutions of the form

u_{i} (x, t) = U_{i} (x_{2}, x_{3}) \exp (i k (x_{1} - v t)), i = 1, 2, 3,

5.3

where k is the wavenumber and v the phase velocity. A closed-form analytical expression for the dispersion relation describing v=v(k) is not available, so the dispersion is usually studied using numerical methods (e.g. [26,27]). We computed numerical solutions of the relevant dispersion relation using the semi-analytical finite-element (SAFE) method, which is a technique that involves the analytic separation of the axial variable followed by the finite-element solution of the suitably parametrized cross-sectional problem (e.g. [28–30] and references therein). All finite-element computations in this paper were performed using the COMSOL Multiphysics 5.1.

The first numerical example that we considered concerns a rectangular beam with cross section 2×1. The dispersion curves for the beam are shown in figure 2a. Our particular interest is in the long-wave behaviour of two lowest (bending) modes. Both of the relevant curves are shown in the shaded area of the graph and, also, presented magnified in figure 2b. Finite-element solutions (the solid curves) are shown together with predictions of the Euler–Bernoulli beam theory (dotted lines), as well as predictions of the new asymptotic equation (4.21) (dash-dotted curves). More specifically, the relevant dispersion relations can be obtained by inserting (5.3) into equation (4.21) and re-arranging the result as

\frac{v^{2}}{E / ρ} = \frac{I_{x_{2}} k^{2} / A}{1 + k_{0} (ν, q) I_{x_{2}} k^{2} / A} = \frac{I_{x_{2}} k^{2}}{A} (1 - k_{0} (ν, q) \frac{I_{x_{2}} k^{2}}{A} + O (\frac{I_{x_{2}}^{2} k^{4}}{A^{2}})),

5.4

where I_x₂=2h₂(2h₃)³/12 and A=4h₂h₃, just like we mentioned in §3. The dispersion relation for the Euler–Bernoulli theory is obtained when the O(I_x₂k²/A) term is omitted on the right-hand side; the dispersion relation for the new beam equation must include this term. Equation (5.4) approximates a single (lower) dispersion curve that corresponds to bending in the plane Ox₁x₃. The bending in the plane Ox₁x₂ can be described by replacing I_x₂ and q in equation (5.4) by I_x₃=(2h₂)³2h₃/12 and 1/q, respectively, see the upper dispersion curve in figure 2b. The improved accuracy of the new equation compared to the accuracy of the classical Euler–Bernoulli equation is evident; the classical theory overestimates the phase velocity, just as expected.

Figure 2. Non-dimensional phase velocity $v / \sqrt{E / ρ}$ as the function of non-dimensional wave number kh₂ for a rectangular beam 2×1 with ν=0.3. For the lower bending mode q=2, for the upper bending mode q=1/2. (a) The dispersion curves obtained using the SAFE method. The shaded area, also shown zoomed as (b), presents a comparison between two bending modes (solid lines) and their approximations using the Euler–Bernoulli theory and asymptotic equation (4.21) (dotted and dash-dotted lines, respectively).

Equation (5.4) predicts no dispersion when k₀(ν,q)=0; of course, this actually means that dispersion in this case would be governed by the next $O (I_{x_{2}}^{2} k^{4} / A^{2})$ term. For example, in a rectangular beam with the cross section 6.2916×1 made of a steel with ν=0.3, the correction term in equation (4.21) vanishes, so one expects to see a substantial improvement in the performance of the Euler–Bernoulli theory for the lower bending mode. For even larger aspect ratios, k₀(ν,q) becomes negative; it does not quite mean the change of the type of dispersion, that e.g. happens for shallow water waves, but it does mean that velocities predicted by the Euler–Bernoulli theory may actually underestimate the true phase velocity of the lower bending mode at low frequencies. This is contrary to the common intuition suggesting that the Euler–Bernoulli theory overestimates the phase velocity in a beam.

In order to illustrate this situation, we also computed dispersion curves for the somewhat extreme case of a beam with cross section 20×1 (figure 3a). Similar to the previous example, we focused on the long-wave low-frequency modes, i.e. modes within the shaded area of the graph that is also presented magnified in figure 3b. The improved approximation accuracy for the upper bending mode, for which q=1/20, is easy to assess from the figure. Phase velocity of the lower bending mode, for which q=20, is so low that all three curves—the numerical dispersion curve, as well as its Euler–Bernoulli approximation and prediction by equation (4.21)—overlap and appear indistinguishable. This is why a separate figure 4 presents numerically computed relative approximation errors for both (a) lower and (b) upper bending modes. The relative error of the Euler–Bernoulli theory for lower mode barely exceeds 1%, which explains our observations (figure 4a). The relative error of the Euler–Bernoulli theory for the higher bending mode is significantly larger (figure 4b). In both cases, the relative errors of the new bending equation (4.21) are one to two orders of magnitude smaller than relative errors of the classical theory.

Figure 4. Relative approximation errors for long-wave low-frequency phase velocity of the bending modes in a rectangular beam 20×1 with ν=0.3. (a) The lower frequency (plate strip) bending mode. (b) The higher frequency in-plane (plane stress) bending mode.

As we assumed the Poisson ratio ν=0.3 and since for lower bending mode q=20, the correction coefficient k₀(0.3,20)≈−46.315<0. Therefore, just as we indicated earlier, for lower bending mode in flat plate-like beams the Euler–Bernoulli theory tends to underestimate the phase velocity (this is related to, but is not the same thing as the negative values for Timishenko's shear correction factor, which will be discussed in §7). Of course, proper modelling of bending waves in such a structure is likely to also require properly accounting for twisting motions of a beam (figure 3b). Therefore, in situations featuring low-frequency propagation/interaction of several bending modes of flat beams, one would probably benefit more from using a plate theory instead of a single mode beam equation such as (4.21).

6. The limiting case of thin flat beams

Let us designate thin, flat flexural elements of rectangular cross section and small (or large) width to thickness aspect ratio as simply thin flat beams. Such beams, e.g. the beam with aspect ratio 20×1 that we studied in §5, are special in that they can be meaningfully modelled using thin plate approximations. It can be gathered from figure 3b that long-wave response of a thin flat beam is dominated by three modes: two bending modes (in-plane and out-of-plane), as well as a twisting mode. Typical beam deformations that correspond to these modes are schematically illustrated in figure 5. The twisting mode shown in figure 5c is distinct from the bending modes in that it is not a low-frequency mode; its long-wave limit corresponds to a relatively low, yet non-zero frequency. This frequency would be a natural limit for the applicability of beam equation (4.21) for thin flat beams. We are not going to otherwise study the twisting mode and leave its analysis to a separate publication.

The lower bending mode (figure 5a) corresponds to the asymptotic limit when q≡h₂/h₃≫1. The infinite sum on the right-hand side of equation (5.2) is exponentially small in this case; hence, $S (\infty) = 372 ζ (5) / π^{5} \approx 1.2604978$ and

k_{0} (ν, q) = - \frac{ν^{2}}{1 + ν} (2 q^{2} - 3 q S (\infty)) + O (1), q ≫ 1 .

6.1

This expansion shows that for all sufficiently flat beams the Euler–Bernoulli theory will underestimate the true phase velocity of plane waves, just as we observed previously for the rectangular beam 20×1.

The correctness of expansion (6.1) can be confirmed using a classical Kirchhoff plate theory. A thin flat beam can be thought of as a plate strip of width 2h₂ and thickness 2h₃; beam bending mode would correspond to the fundamental symmetric mode of the strip. Motion of the plate strip of thickness 2h₃ is governed by the equation

\frac{2 E h_{3}^{3}}{3 (1 - ν^{2})} (\frac{\partial^{4} u_{3}}{\partial x_{1}^{4}} + 2 \frac{\partial^{4} u_{3}}{\partial x_{1}^{2} \partial x_{2}^{2}} + \frac{\partial^{4} u_{3}}{\partial x_{2}^{4}}) + 2 ρ h_{3} \frac{\partial^{2} u_{3}}{\partial t^{2}} = 0 .

6.2

The boundary conditions that ensure that strip edges are stress-free can be written in the form

\frac{\partial^{2} u_{3}}{\partial x_{2}^{2}} + ν \frac{\partial^{2} u_{3}}{\partial x_{1}^{2}} + K h_{3} (1 - ν) \frac{\partial^{3} u_{3}}{\partial x_{1}^{2} \partial x_{2}} = 0

6.3

and

\frac{\partial^{3} u_{3}}{\partial x_{2}^{3}} - (2 - ν) \frac{\partial^{3} u_{3}}{\partial x_{1}^{2} \partial x_{2}} = 0,

6.4

where x₂=±h₂. These conditions are unusual due to the presence of O(h₃) term in (6.3), which represents a first non-trivial correction to the classical Kirchhoff boundary conditions by Kolos [18]. Without this correction term, the classical leading-order theory of plate bending would only be able to capture the leading-order behaviour of expansion (6.1); with the correction term added we expect to recover (6.1) fully. Indeed, Kolos [18, eqn (2.13)] showed that the constant K is defined as

K = \frac{384}{π^{5}} \sum_{n = 1}^{\infty} \frac{1}{{(2 n - 1)}^{5}} .

6.5

That is, in our notation

K \equiv S (\infty)

, which immediately suggests the potential relevance of the correction term.

Solutions of the form

u_{3} = (U_{3}^{(1)} \cosh (α_{1} k h_{2}) + U_{3}^{(2)} \cosh (α_{2} k h_{2})) \exp (i k (x_{1} - v t)),

6.6

satisfy governing equation (6.2) as long as

α_{1, 2}^{2} = 1 \pm \frac{\sqrt{3 (1 - ν^{2})}}{k h_{3}} \frac{v}{\sqrt{E / ρ}} .

6.7

Substitution into boundary conditions (6.3) and (6.4) shows that non-trivial solutions of form (6.6) only exist when the following secular equation is satisfied:

g (α_{1}, α_{2}, ν) C_{1} S_{2} - g (α_{2}, α_{1}, ν) S_{1} C_{2} + k h_{3} K (1 - ν) α_{1} α_{2} (α_{1}^{2} - α_{2}^{2}) S_{1} S_{2} = 0,

6.8

within which

S_{i} = \sinh (α_{i} k h_{2}), C_{i} = \cosh (α_{i} k h_{2})

6.9

and

g (α_{i}, α_{j}, ν) = α_{j} (α_{i}^{2} - ν) (α_{j}^{2} + ν - 2), i \neq j, i, j = 1, 2 .

6.10

If K was equal to zero, equation (6.8) would reduce to the dispersion relation first obtained and analysed by Konenkov [31]. As is, this secular equation represents a generalization of Konenkov's result to the case of non-vanishing ratios h₃/h₂.

Although (6.8) cannot, in general, be solved analytically, formal asymptotic expansions are possible and expansion in the long-wave low-frequency limit has the following form:

\frac{v^{2}}{E / ρ} = \frac{I_{x_{2}} k^{2}}{A} (1 + \frac{ν^{2} (2 q^{2} - 3 K q)}{1 + ν} \frac{I_{x_{2}} k^{2}}{A} + O (\frac{I_{x_{2}}^{2} k^{4}}{A^{2}})),

6.11

which in view of equations (5.4) and (6.5) is clearly equivalent to (6.1). This serves to confirm the validity of (5.4), at least in the asymptotic sense for q≫1, and, simultaneously, this serves as a non-trivial confirmation of the correctness of the boundary layer adjustment term by Kolos [18] for the free boundary conditions in the Kirchhoff theory.

The higher bending mode, sketched in figure 5b, is associated with the values of parameter q≡h₂/h₃≪1. A mode like this can also be extrapolated by using a plate theory; however, the bending deformation in this case would be confined to a plane and, therefore, would require an application of the plane stress theory. Timoshenko [19] and Stephen [11] considered this limiting behaviour and obtained, in different notations, the following expansion for the phase velocity:

\frac{v^{2}}{E / ρ} = \frac{I_{x_{2}} k^{2}}{A} (1 - (2 ν + \frac{17}{5}) \frac{I_{x_{2}} k^{2}}{A} + O (\frac{I_{x_{2}}^{2} k^{4}}{A^{2}})) .

6.12

For our new beam bending equation, this corresponds to the limit when q≡h₂/h₃≪1; hence, in view of expansion (5.1), we can say that

k_{0} (ν, q) = 2 ν + \frac{17}{5} - \frac{2 ν^{2}}{5 (1 + ν)} q^{4} + \dots, q ≪ 1 .

6.13

A quick glance at dispersion relation (5.4) now confirms that Stephen's expansion (6.12) gives the correct leading-order behaviour for the dispersion coefficient k₀(ν,q) when q≪1.

It would be interesting to explore whether a refined plane stress theory would be sufficient to also recover the O(q⁴) correction term in (6.13). The relevant higher-order equations are available in [4] and the appropriately refined boundary conditions are described in [18]. However, this discussion would take us too far outside of the scope of this paper.

7. Shear correction factor

In derivation of his beam equation, Timoshenko relied upon two key constitutive assumptions:

M = - E I_{x_{2}} \frac{\partial ϕ}{\partial x_{1}} and Q = κ μ A (\frac{\partial w}{\partial x_{1}} - ϕ),

7.1

where M is the bending moment, ϕ the angle characterizing rotation of a cross-sectional element, Q the shearing force, A the cross-sectional area and μ=E/2(1+ν) the shear modulus, see [8]. Assumption (7.1)₁ encapsulates the same leading-order relationship between bending moment M and the radius of curvature (∂ϕ/∂x₁)⁻¹ as in the Euler–Bernoulli theory. Assumption (7.1)₂ that relates shearing force at the cross section to transverse shear strain at the centroidal axis is unique to Timoshenko's theory. The proportionality constant κ—‘the shear correction factor’—had to be formally introduced to bring into accord the relationship that was known to Timoshenko to be approximate at best. The associated beam equation, unsurprisingly, also depends on κ:

E I_{x_{2}} \frac{\partial^{4} w}{\partial x_{1}^{4}} + ρ A \frac{\partial^{2} w}{\partial t^{2}} - ρ I_{x_{2}} (1 + \frac{E}{κ μ}) \frac{\partial^{4} w}{\partial x_{1}^{2} \partial t^{2}} + \frac{ρ^{2} I_{x_{2}}}{κ μ} \frac{\partial^{4} w}{\partial t^{4}} = 0,

7.2

(see [9]). It is important to stress that, even though equation (7.2) contains higher-order terms compared to the Euler–Bernoulli equation, the associated higher-order corrections are not necessarily consistent with the appropriately truncated expressions from the three-dimensional theory. This happens because the underlying assumptions (7.1) are not accurate to the same order of truncation, see [9] and also discussions of similar issues in refined theories of plates and shells in [3,32]. More specifically, even if the constitutive relation for shearing force (7.1)₂ is made asymptotically consistent, the recovery of the appropriately refined displacement and stress fields would also require a higher-order generalization of the constitutive relation for bending moment (7.1)₁ as well as appropriate refinements to the transverse inertia. This can be deduced from asymptotic formulae in § 4, see also [3,32] and references therein for more details. As a result, even with the presence of ‘tuning’ parameter κ, one does not expect to be able to recover correct asymptotic behaviour of every aspect of the resulting theory. In fact, later in this section, we will demonstrate how the choice of κ that increases only the numerical accuracy of equation (7.1)₂ is resulting in inconsistent governing equation (7.2).

A sizeable body of the literature sprung up from attempts to either formulate a consistent method to derive κ or to justify a particular expression for it. Timoshenko himself used several cross-section-independent constant values [8,19]; a surprisingly influential paper by Cowper [20] used a series of approximations to arrive at a general expression for κ that was then specialized to a number of canonical cross sections; for rectangular beams, it gives

κ_{C} = \frac{10 (1 + ν)}{12 + 11 ν} .

7.3

Stephen [11] used several long-wave low-frequency expansions (including (6.12)) to show that Cowper's expressions cannot be consistent, at least for circular and rectangular cross sections.

The same idea can also be applied to our beam equation (4.21) to derive an asymptotically consistent expression for the Timoshenko shear correction factor. The appropriate dispersion relation can be obtained by substituting plane wave solution (5.3) into (7.2) and expanding the result for small values of I_x₂k²/A. This procedure yields the following approximate dispersion relation:

\frac{v^{2}}{E / ρ} = \frac{I_{x_{2}} k^{2}}{A} (1 - (1 + \frac{E}{κ μ}) \frac{I_{x_{2}} k^{2}}{A} + O (\frac{I_{x_{2}}^{2} k^{4}}{A^{2}})) .

7.4

A direct comparison of expansions (7.4) and (5.4) indicates that Timoshenko's equation will be asymptotically equivalent to beam equation (4.21) as long as

\begin{aligned} κ_{NPK} & = \frac{2 (1 + ν)}{k_{0} (ν, q) - 1} = \frac{5 {(1 + ν)}^{2}}{6 + 11 ν - ν^{2} - 5 ν^{2} q^{2} + \frac{2880 ν^{2} q}{π^{5}} \sum_{n = 1}^{\infty} \frac{\coth (q π \frac{2 n - 1}{2})}{{(2 n - 1)}^{5}}}, \end{aligned}

7.5

which appears to be a new expression for the Timoshenko shear correction factor. This expression can be compared with other expressions for κ in the literature. Stephen & Levinson [21, eqn (54b)] published an ad hoc derivation that led to the following expression:

κ_{SL} = \frac{5 {(1 + ν)}^{2}}{6 + 11 ν + ν^{2} (5 - q^{4} + \frac{90 q^{5}}{π^{5}} \sum_{n = 1}^{\infty} \frac{\tanh (π n / q)}{n^{5}})},

7.6

which, although appears different to (7.5), is equivalent to (7.5), see the demonstration in appendix A. The same expression for κ was confirmed by Stephen [22, eqn (31)] and, in a somewhat different form, by Hutchinson [23, eqn (47)], see the discussion in [24]. In addition, a direct perturbation in powers of thickness variable was developed by Chan et al. [25], which also confirmed the correctness of (7.6). Based on these comments, we conclude that κ_NPK=κ_SL and that equations (7.5) and (7.6) are equivalent and asymptotically correct forms for the shear correction factor for a rectangular Timoshenko beam.

The infinite series used in equations (7.5) and (7.6) will have to be truncated for numerical computations. If these series are truncated after the same number of terms, equation (7.5) provides a more precise answer for $q ≳ 1$ and equation (7.6) provides a more precise answer for $q ≲ 1$ . However, the precision of (7.6) with truncated series quickly deteriorates for larger values of q, whereas the precision of (7.5) with truncated series stays approximately uniform for all values of q (e.g. the first 10 terms in the series are typically sufficient to obtain five significant digits of κ). Thus, equation (7.5) can be recommended as more universal expression for computing κ.

Berdichevskii & Kvashnina made a similar comparison of their asymptotic beam theory [12, eqn (3.2)] with the Timoshenko equation. This led them to derive a general expression for κ valid for sufficiently symmetric cross sections, see [12, eqn (8.6)]. Their resulting coefficient for circular cross sections [12, eqn (8.7)] matches the coefficient derived from three-dimensional theory by Stephen [11] and since then confirmed by a number of other researchers. At the same time, their coefficient for rectangular cross sections [12, eqn (8.9)]:

\begin{aligned} κ_{BK} = \frac{1}{\frac{191 ν^{2} + 594 ν + 324}{270 {(1 + ν)}^{2}} + \frac{ν^{2}}{{(1 + ν)}^{2}} (\frac{5}{54 q^{4}} + \frac{6 q}{π^{5}} \sum_{n = 1}^{\infty} \frac{\tanh (π n / q)}{n^{5}})}, \end{aligned}

7.7

despite certain similarities, is not equivalent to (7.5) or (7.6) and, therefore, appears to be in error.

Substantial research effort has been spent trying to evaluate the value of κ that would render approximate formula (7.1)₂ numerically exact for rectangular beams. The earliest relevant result in the literature was obtained by Eliseev [33, eqn (33)]; it can be re-written in our notation as:

κ_{E} = \frac{{(1 + ν)}^{2}}{\frac{6}{5} (1 + 2 ν + 2 ν^{2}) - \frac{576 ν^{2} q}{π^{5}} \sum_{n = 1}^{\infty} \frac{\coth (q π \frac{2 n - 1}{2})}{{(2 n - 1)}^{5}} + ν^{2} q^{2}} .

7.8

Fully equivalent, although very different-looking, shear correction factors were obtained by Renton [34], just above eqn (27):

κ_{R} = \frac{1}{\frac{6}{5} + {(\frac{ν}{1 + ν})}^{2} \sum_{m = 0}^{\infty} \sum_{n = 1}^{\infty} \frac{36 q^{4}}{π^{6} {(2 m + 1)}^{2} n^{2} (4 n^{2} + {(2 m + 1)}^{2} q^{2})}},

7.9

and by Yu & Hodges [14, eqn (68)₁]:

κ_{YH} = \frac{1}{\frac{6}{5} + q^{4} {(\frac{ν}{1 + ν})}^{2} (\frac{1}{5} - \frac{18 q}{π^{5}} \sum_{n = 1}^{\infty} \frac{\tanh (π n / q)}{n^{5}})} .

7.10

The equivalence of expressions (7.8)–(7.10) can be easily demonstrated using the identities described in appendix A, so κ_E=κ_R=κ_YH.

Despite some superficial similarities between the expressions for κ_NPK=κ_SL and expressions for κ_E=κ_R=κ_YH, it is important to stress that they are very different functions of the beam aspect ratio q and the Poisson ratio ν. Figure 6 illustrates this point for the Poisson ratio ν=0.3. It is immediately obvious that expressions κ_E=κ_R=κ_YH do not result in the correct limiting behaviour for q≪1, just like κ_C did not. The definition of κ that leads to the asymptotically consistent governing equation is actually discontinuous when k₀(ν,q)=1, see (7.5); when ν=0.3 this happens for q≈5.5671. The asymptotically correct value of κ for even flatter beams is negative, which suggests anti-spring relationship between shearing force and transverse shear strain. The physical meaning of equation (7.1)₂ for negative or discontinuous values of κ is dubious; however, rather than casting the doubt on definitions of κ_NPK=κ_SL, this simply reflects the fact that (7.1)₂ is not a consistent statement from the point of view of long-wave low-frequency asymptotics.

Figure 6. The shear correction factor κ for rectangular beams when the Poisson ratio ν=0.3. Parameter q is the aspect ratio of the cross section rectangle.

In fact, it is possible to estimate the ‘true’ value of κ from the dispersion curves obtained using the finite-element method. Indeed, we expect that Timoshenko's dispersion relation (7.4) must provide a fair approximation of the exact dispersion curves for sufficiently small wavenumbers. Therefore, the difference between the phase velocity of Euler–Bernoulli beam $v_{EB} = \sqrt{E / ρ} I_{x_{2}} k^{2} / A$ and the phase velocity resulting from the full three-dimensional theory v_FEM must be equal, to the leading order, to the higher-order correction term in equation (7.4). A simple algebraic manipulation can then be used to show that

κ_{FEM} \approx \frac{2 (1 + ν)}{(2 \frac{v_{EB} - v_{FEM}}{\sqrt{E / ρ} {(I_{x_{2}} k^{2} / A)}^{3}} - 1)} .

7.11

This identity is, of course, approximate, but it becomes more and more accurate as I_x₂k²/A→0. Thus, the long-wave limit of (7.11) can be used to estimate the values of κ that correspond to the considered dispersion curve. We applied the described procedure to all four bending mode dispersion curves discussed in §5 and obtained estimates for the shear correction factor. These estimates are denoted as κ_FEM in figure 6 and they can also be seen tabulated in table 1. Also included in table 1 are the values of the shear correction factor that are obtained from the expressions discussed in this section. From these results, it is abundantly clear that only the definitions of κ_NPK=κ_SL are truly consistent with the results of finite-element simulations.

Table 1.Specific values of the shear correction factor for ν=0.3.
Collapse
q	κ_C	κ_E=κ_R=κ_YH	κ_NPK=κ_SL	κ_FEM
1/20	0.8497	0.8333	0.8667	0.87
1/2	0.8497	0.8329	0.8671	0.87
2	0.8497	0.7844	0.9267	0.93
20	0.8497	0.04866	−0.05495	−0.055

8. Summary

An asymptotic beam equation constructed in this paper can be used to describe the long-wave low-frequency behaviour of both bending modes of a rectangular beam without making any assumptions about aspect ratio of the beam. If we assume that h₂≥h₃ or, equivalently, that q≡h₂/h₃≥1, then the lower bending mode of the beam is governed by

\frac{E I_{x_{2}}}{ρ A} \frac{\partial^{4} w}{\partial x_{1}^{4}} + \frac{\partial^{2} w}{\partial t^{2}} - k_{0} (ν, q) \frac{I_{x_{2}}}{A} \frac{\partial^{4} w}{\partial x_{1}^{2} \partial t^{2}} = 0 .

8.1

in which I_x₂=2h₂(2h₃)³/12, A=4h₂h₃ and the non-dimensional coefficient k₀(ν,q) is defined by equations (4.20) and (4.17)₂. The governing equation for higher bending mode can be obtained by replacing I_x₂ and q with I_x₃=(2h₂)³2h₃/12 and 1/q, respectively.

We used a variety of tests to demonstrate the asymptotic consistency of the newly derived equation. Specifically, we correctly recovered the asymptotic behaviour in the limit of thin flat beams modelled by a strip of Kirchhoff plate. Very remarkably, when using the classical plate theory with the Kirchhoff boundary conditions refined by Kolos [18], the first two terms of the corresponding dispersion relation were shown to fully match the dispersion relation for the new beam theory. This appears to be the first non-trivial confirmation of the correctness of Kolos’ result. We also showed that the new equation correctly reproduces the limiting case of narrow beams, which can be obtained by analysing a strip problem for a thin plate in the state of plane stress.

A comparison with the recent literature on the Timoshenko shear correction factor κ was also made. A new asymptotically consistent expression for κ was obtained in the following form:

κ = \frac{5 {(1 + ν)}^{2}}{6 + 11 ν - ν^{2} - 5 ν^{2} q^{2} + \frac{2880 ν^{2} q}{π^{5}} \sum_{n = 1}^{\infty} \frac{\coth (q π \frac{2 n - 1}{2})}{{(2 n - 1)}^{5}}} .

8.2

We showed that equation (8.2), although apparently new, is equivalent to the expressions for κ previously derived in [21–23,25], see also discussion [24]. Given the substantial differences between the methods used to obtain these expressions, our study may serve as a confirmation of the asymptotic consistency of the previously reported results.

The relative simplicity of rectangular cross section allowed us to obtain a number of explicit analytical results for rectangular beams and cross-verify predictions that resulted from considering distinct asymptotic limiting processes. The case of more general beam cross sections will necessarily be less explicit and will become the subject of a separate treatment. However, the main benefit of the asymptotically consistent beam theory, such as the one derived in this paper, is the possibility of extending our results to finite beams. Such an extension must necessarily involve an appropriate refinement of the end boundary conditions [35]. There exists a substantial body of research on refinements of edge boundary conditions for plate theories (e.g. [36–38] and references therein). We are hoping to address the sparsity of similar results for higher-order beam theories in our future work.

Data accessibility

This article has no additional data.

Authors' contributions

All authors participated in deriving the results reported in this article. All authors contributed to writing the manuscript. Each of the authors gave their final approval for publication.

Competing interests

We have no competing interests.

Funding

No funding has been received for this article.

Appendix A. Demonstration of equivalence of (7.5) and (7.6)

The shear correction factors (7.5) and (7.6) feature infinite sums that must be expressed in terms of each other if we are to compare them. The necessary result can be established by using the following series for hyperbolic functions [39], 1.421.2 and 1.421.4:

\coth π x = \frac{1}{π x} + \frac{2 x}{π} \sum_{k = 1}^{\infty} \frac{1}{x^{2} + k^{2}}

A 1

and

\tanh \frac{π x}{2} = \frac{4 x}{π} \sum_{k = 1}^{\infty} \frac{1}{{(2 k - 1)}^{2} + x^{2}} .

A 2

We will also need several identities related to the Riemann zeta function. First, we note that for s>1

ζ (s) = \sum_{k = 1}^{\infty} \frac{1}{k^{s}} = \sum_{k = 1}^{\infty} [\frac{1}{{(2 k - 1)}^{s}} + \frac{1}{2^{s} k^{s}}] = \sum_{k = 1}^{\infty} \frac{1}{{(2 k - 1)}^{s}} + \frac{ζ (s)}{2^{s}},

A 3

see [39, Sect. 9.52, 9.54]. Consequently,

\sum_{k = 1}^{\infty} 2^{s} / {(2 k - 1)}^{s} = (2^{s} - 1) ζ (s)

, and since ζ(6)=π⁶/945 we can conclude that

\sum_{k = 1}^{\infty} \frac{2^{6}}{{(2 k - 1)}^{6}} = \frac{π^{6}}{15} .

A 4

Similarly, it is easy to see that for s>1 and r>1

\sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} \frac{1}{n^{s}} \frac{2^{r}}{{(2 k - 1)}^{r}} = (2^{r} - 1) ζ (s) ζ (r);

A 5

hence, because ζ(2)=π²/6 and ζ(4)=π⁴/90,

\sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} \frac{1}{n^{2}} \frac{2^{4}}{{(2 k - 1)}^{4}} = \frac{π^{6}}{36} and \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} \frac{1}{n^{4}} \frac{2^{2}}{{(2 k - 1)}^{2}} = \frac{π^{6}}{180} .

A 6

Identities (A.4) and (A.6) can also be computed directly using a computer algebra system such as Waterloo Maple. With all these equations in mind, one can construct the following sequence of transformations:

\begin{aligned} 2^{5} q \sum_{k = 1}^{\infty} \frac{\coth (π (2 k - 1) q / 2)}{{(2 k - 1)}^{5}} & \overset{(A 1)}{=} \frac{1}{π} \sum_{k = 1}^{\infty} \frac{2^{6}}{{(2 k - 1)}^{6}} [1 + 2 {(2 k - 1)}^{2} \sum_{n = 1}^{\infty} \frac{1}{{(2 k - 1)}^{2} + {(2 n / q)}^{2}}] \\ \overset{(A 4)}{=} \frac{π^{5}}{15} + \frac{2^{7}}{π} \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} \frac{1}{{(2 k - 1)}^{4}} \frac{1}{{(2 k - 1)}^{2} + {(2 n / q)}^{2}} \\ = \frac{π^{5}}{15} + \frac{2^{7}}{π} \sum_{n = 1}^{\infty} \sum_{k = 1}^{\infty} [\frac{{(2 n / q)}^{- 2}}{{(2 k - 1)}^{4}} - \frac{{(2 n / q)}^{- 4}}{{(2 k - 1)}^{2}} + \frac{{(2 n / q)}^{- 4}}{{(2 k - 1)}^{2} + {(2 n / q)}^{2}}] \\ \overset{(A 6)}{=} \frac{π^{5}}{15} + \frac{π^{5} q^{2}}{18} - \frac{π^{5} q^{4}}{90} + \frac{2^{3} q^{4}}{π} \sum_{n = 1}^{\infty} \frac{1}{n^{4}} \sum_{k = 1}^{\infty} \frac{1}{{(2 k - 1)}^{2} + {(2 n / q)}^{2}} \\ \overset{(A 2)}{=} \frac{π^{5}}{90} (6 + 5 q^{2} - q^{4}) + q^{5} \sum_{n = 1}^{\infty} \frac{\tanh (π n / q)}{n^{5}} . \end{aligned}

A 7

In the view of identity (A.7), the equivalence of expressions (7.5) and (7.6) can now be established by elementary algebraic manipulations.

Footnotes

1 A possible way to introduce h could be to define it as $h = \frac{1}{2} (h_{2} + h_{3})$ or, alternatively, as $h = max (h_{2}, h_{3})$ . In our derivations, we implicitly assume that h₂ and h₃ are of the same asymptotic order; however, the resulting expansions remain valid if one of them is substantially smaller than the other.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

References

1
Berdichevskii VL. 1972Equations describing transverse vibration of thin elastic plates. Mech. Solids 7, 138–141. (transl. from Izv. AN SSSR. MTT vol. 1972(6), pp. 152–155). Google Scholar
2
Berdichevskii VL. 1979Variational-asymptotic method of constructing a theory of shells. J. Appl. Math. Mech. 43, 711–736. (transl. from PMM vol. 43(4), pp. 664–687). (doi:10.1016/0021-8928(79)90157-6) Crossref, Google Scholar
3
Goldenveizer AL, Kaplunov JD, Nolde EV. 1993On Timoshenko-Reissner type theories of plates and shells. Int. J. Solids Struct. 30, 675–694. (doi:10.1016/0020-7683(93)90029-7) Crossref, ISI, Google Scholar
4
Kaplunov JD, Kossovich LYu, Nolde EV. 1998Dynamics of thin walled elastic bodies. New York, NY: Academic Press. Google Scholar
5
Kaplunov JD, Nolde EV, Shorr BF. 2005A perturbation approach for evaluating natural frequencies of moderately thick elliptic plates. J. Sound Vib. 281, 905–919. (doi:10.1016/j.jsv.2004.02.046) Crossref, ISI, Google Scholar
6
Belov AV, Kaplunov JD, Nolde EV. 1998A refined asymptotic model of fluid-structure interaction in scattering by elastic shells. Flow Turbul. Combust. 61, 255–267. (doi:10.1023/A:1026437810379) Crossref, ISI, Google Scholar
7
Kaplunov J, Nolde E, Rogerson GA. 2006An asymptotic analysis of initial-value problems for thin elastic plates. Proc. R. Soc. A 462, 2541–2561. (doi:10.1098/rspa.2006.1687) Link, Google Scholar
8
Timoshenko SP. 1921On the correction for shear of the differential equation for transverse vibrations of prismatic bars. Philos. Mag. Ser. VI 41, 744–746. (doi:10.1080/14786442108636264) Crossref, Google Scholar
9
Elishakoff I, Kaplunov J, Nolde E. 2015Celebrating the centenary of Timoshenko's study of effects of shear deformation and rotary inertia. Appl. Mech. Rev. 67, 060802 (11 pages). (doi:10.1115/1.4031965) Crossref, ISI, Google Scholar
10
Stephen NG. 2006The second spectrum of Timoshenko beam theory—further assessment. J. Sound Vib. 292, 372–389. (doi:10.1016/j.jsv.2005.08.003) Crossref, ISI, Google Scholar
11
Stephen NG. 1981Considerations on second order beam theories. Int. J. Solids Struct. 17, 325–333. (doi:10.1016/0020-7683(81)90066-4) Crossref, ISI, Google Scholar
12
Berdichevskii VL, Kvashnina SS. 1976On equations describing the transverse vibrations of elastic bars. J. Appl. Math. Mech. 40, 104–119. (transl. from PMM vol. 40(1), pp. 120–135). (doi:10.1016/0021-8928(76)90117-9) Crossref, Google Scholar
13
Berdichevskii VL, Starosel'skii LA. 1983On the theory of curvilinear Timoshenko-type rods. J. Appl. Math. Mech. 47, 809–817. (transl. from PMM vol. 47(6), pp. 1015–1024). (doi:10.1016/0021-8928(83)90121-1) Crossref, Google Scholar
14
Yu W, Hodges DH. 2004Elasticity solutions versus asymptotic sectional analysis of homogeneous, isotropic, prismatic beams. ASME J. Appl. Mech. 71, 15–23. (doi:10.1115/1.1640367) Crossref, ISI, Google Scholar
15
Hodges DH. 2006Nonlinear composite beam theory. Progress in Astronautics and Aeronautics. New York, NY: AIAA. Google Scholar
16
Goldenveizer AL. 1966The principles of reducing three-dimensional problems of elasticity to two-dimensional problems of the theory of plates and shells. In Applied Mechanics: Proc. of the Eleventh Int. Congress of Applied Mechanics, Munich (Germany) 1964 (eds H Görtler, P Sorger), pp. 306–311. Berlin, Germany: Springer. Google Scholar
17
Sayir M. 1980Flexural vibrations of strongly anisotropic beams. Archive of Appl. Mech. 49, 309–323. (doi:10.1007/BF02426910) Google Scholar
18
Kolos AV. 1965Methods of refining the classical theory of bending and extension of plates. J. Appl. Math. Mech. 29, 914–925. (transl. from PMM vol. 29(4), pp. 771–781). (doi:10.1016/0021-8928(65)90104-8) Crossref, Google Scholar
19
Timoshenko SP. 1922On the transverse vibrations of bars of uniform cross-section. Philos. Mag. Ser. VI 43, 125–131. (doi:10.1080/14786442208633855) Crossref, Google Scholar
20
Cowper GR. 1966The shear coefficient in Timoshenko's beam theory. ASME J. Appl. Mech. 33, 335–340. (doi:10.1115/1.3625046) Crossref, ISI, Google Scholar
21
Stephen NG, Levinson M. 1979A second order beam theory. J. Sound Vib. 67, 293–305. (doi:10.1016/0022-460X(79)90537-6) Crossref, ISI, Google Scholar
22
Stephen NG. 1980Timoshenko's shear coefficient from a beam subjected to gravity loading. ASME J. Appl. Mech. 47, 121–127. (doi:10.1115/1.3153589) Crossref, ISI, Google Scholar
23
Hutchinson JR. 2001Shear coefficients for Timoshenko beam theory. ASME J. Appl. Mech. 68, 87–92. (doi:10.1115/1.1349417) Crossref, ISI, Google Scholar
24
Stephen NG, Hutchinson JR. 2001Discussion: ‘Shear coefficients for Timoshenko beam theory’ (Hutchinson, J.R., 2001, ASME J. Appl. Mech., 68, pp. 87–92). ASME J. Appl. Mech. 68, 959–961. (doi:10.1115/1.1412454) Crossref, ISI, Google Scholar
25
Chan KT, Lai KF, Stephen NG, Young K. 2011A new method to determine the shear coefficient of Timoshenko beam theory. J. Sound Vib. 330, 3488–3497. (doi:10.1016/j.jsv.2011.02.012) Crossref, ISI, Google Scholar
26
Tanaka K, Iwahashi Y. 1977Dispersion relation of elastic waves in bars of rectangular cross section. Bull. Jpn. Soc. Mech. Eng. 20, 922–929. (doi:10.1299/jsme1958.20.922) Crossref, Google Scholar
27
Krushynska AA, Meleshko VV. 2011Normal waves in elastic bars of rectangular cross section. J. Acoust. Soc. Am. 129, 1324–1335. (doi:10.1121/1.3531800) Crossref, PubMed, ISI, Google Scholar
28
Hayashi T, Song WJ, Rose JL. 2003Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example. Ultrasonics 41, 175–183. (doi:10.1016/S0041-624X(03)00097-0) Crossref, PubMed, ISI, Google Scholar
29
Damljanović V, Weaver RL. 2004Propagating and evanescent elastic waves in cylindrical waveguides of arbitrary cross section. J. Acoust. Soc. Am. 115, 1572–1581. (doi:10.1121/1.1687424) Crossref, ISI, Google Scholar
30
Predoi M, Castaings M, Hosten B, Bacon C. 2007Wave propagation along transversely periodic structures. J. Acoust. Soc. Am. 121, 1935–1944. (doi:10.1121/1.2534256) Crossref, PubMed, ISI, Google Scholar
31
Konenkov YuK. 1960Plate waves and flexural oscillations of a plate. Soviet Physics. Acoustics 6, 52–59. (transl. from Akusticheskii Zhurnal vol. 6(1), pp. 57–64). Google Scholar
32
Goldenveizer AL, Kaplunov JD, Nolde EV. 1990Asymptotic analysis and improvements of Timoshenko-Reissner-type theories of plates and shells. Mech. Solids 25, 126–139. (transl. from Izv. AN SSSR. MTT vol. 1990(6), pp. 124–138). Google Scholar
33
Eliseev VV. 1989Constitutive equations for elastic prismatic bars. Mech. Solids 24, 66–71. (transl. from Izv. AN SSSR. MTT vol. 1989(1), pp. 70–75). Google Scholar
34
Renton JD. 1991Generalized beam theory applied to shear stiffness. Int. J. Solids Struct. 27, 1955–1967. (doi:10.1016/0020-7683(91)90188-L) Crossref, ISI, Google Scholar
35
Nicholson JW, Simmonds JG. 1977Timoshenko beam theory is not always more accurate than elementary beam theory. J. Appl. Mech. 44, 337–338. (doi:10.1115/1.3424048) Crossref, ISI, Google Scholar
36
Gregory RD, Wan FYM. 1985On plate theories and Saint-Venant's principle. Int. J. Solids Struct. 21, 1005–1024. (doi:10.1016/0020-7683(85)90052-6) Crossref, ISI, Google Scholar
37
Gol’denveizer AL. 1994Algorithms of the asymptotic construction of linear two-dimensional thin shell theory and the St Venant principle. J. Appl. Math. Mech. 58, 1039–1050. (transl. from PMM vol. 58(6), pp. 96–108). (doi:10.1016/0021-8928(94)90120-1) Crossref, Google Scholar
38
Babenkova E, Kaplunov J. 2004Low-frequency decay conditions for a semi-infinite elastic strip. Proc. R. Soc. Lond. A 460, 2153–2169. (doi:10.1098/rspa.2003.1275) Link, ISI, Google Scholar
39
Gradshteyn IS, Ryzhik IM. 2007Table of integrals, series, and products(ed. A. Jeffrey, D. Zwillinger), 7th edn. New York, NY: Elsevier. Google Scholar