An asymptotic higher-order theory for rectangular beams

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.


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, 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 (crosssection-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][23][24][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.

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 where τ ij = τ ij (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 1 along the symmetry axis. The beam cross section is the rectangle {(x 2 , x 3 ) : −h 2 ≤ x 2 ≤ h 2 , −h 3 ≤ x 3 ≤ h 3 } (figure 1). The boundary conditions on the free faces of the beam may then be written as τ 2i = 0 (i = 1, 2, 3) at x 2 = ±h 2 and τ 3i = 0 (i = 1, 2, 3) at x 3 = ±h 3 . (

2.2)
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. 1 This suggests introducing a natural small parameter ε = h/l and re-scaling the problem in terms of the following nondimensional coordinates: An additional assumption that the beam bends along the axis Ox 3 prompts an appropriate rescaling for the displacement and stress components: u 1 = εlu, u 2 = ε 2 lv, u 3 = lw, τ 11 = εEσ 11 , τ 33 = ε 3 Eσ 33 , τ 1i = ε 2 Eσ 1i , τ 2i = ε 3 Eσ 2i (i = 2, 3). (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 1 A possible way to introduce h could be to define it as h = 1 2 (h 2 + h 3 ) or, alternatively, as h = max(h 2 , h 3 ). In our derivations, we implicitly assume that h 2 and h 3 are of the same asymptotic order; however, the resulting expansions remain valid if one of them is substantially smaller than the other.  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: and (2.6) Boundary conditions (2.2) can also be appropriately re-formulated as where a = h 2 /h and b = h 3 /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: where f (k) = f (k) (ξ , η, ζ , τ ) and

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(ε 2 ) (and smaller), which yields The appropriate solution is given by wherew 0 =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 and ,η + νw (2) ,ζ , ,η + νw (2) ,ζ , ,ζ , (3.8) ,η + v ,ξ , (3.9) ,ζ + w (2) ,ξ , (3.10) ,ζ + w (2) ,η . (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: wherew 2 =w 2 (ξ , τ ). Governing equation (3.6) then gives Functionw 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 (3.14) In the view of the boundary conditions at η = ±a and ζ = ±b, this yields which can be recognized as the classical Euler-Bernoulli equation for a rectangular beam. Indeed, given thatw 0 ∼ w = u 3 /l, we can formulate the dimensional equivalent of equation (3.16) as Keeping in mind that the second moment for rectangular cross section I x 2 = 2h 2 (2h 3 ) 3 /12 and the cross-sectional area A = 4h 2 h 3 , the coefficient in front of the fourth-order term can now be recognized as the usual constant EI x 2 /ρA. To obtain an expression for u (2) , we first express σ 13 | ζ =±b = 0. As a result of these manipulations, we are led to the following boundary value problem: u (2) ,ηη + u (2) ,ζ ζ = 2w 0 ζ , (3.18) u (2) ,η | η=±a = ∓νw 0 aζ , u , We are seeking the solution of problem (3.18) in the form It is easy to verify that if one takes (3.20) then the problem (3.18) is reduced to the boundary value problem The solution of (3.21) is given by within which As a result, we have It can be readily verified that substitution of (3.24) 2 into (3.5) and integration over the cross section R once again gives us equation (3.16) that governsw 0 . Therefore, we have now satisfied all of the O(ε 2 ) 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.

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 0 (ν, q), see (4.17) 2 . First, given a fixed q, the infinite series within (4.17) 2 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 A useful representation of (4.17) 2 for q 1 can be formulated by using the fact that coth x = 1 + 2/(e 2x − 1); the results outlined in appendix A, especially equation (A 4), indicate that 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 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][29][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 where Equation (5.4) predicts no dispersion when k 0 (ν, q) = 0; of course, this actually means that dispersion in this case would be governed by the next O(I 2 x 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 0 (ν, 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.
As we assumed the Poisson ratio ν = 0.3 and since for lower bending mode q = 20, the correction coefficient k 0 (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).

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 2 /h 3 1. The infinite sum on the right-hand side of equation (5.2) is exponentially small in this case; hence, S(∞) = 372ζ (5)/π 5 ≈ 1.2604978 and 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 2 and thickness 2h 3 ; beam bending mode would correspond to the fundamental symmetric mode of the strip. Motion of the plate strip of thickness 2h 3 is governed by the equation The boundary conditions that ensure that strip edges are stress-free can be written in the form and where x 2 = ±h 2 . These conditions are unusual due to the presence of O(h 3 ) 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 That is, in our notation K ≡ S(∞), which immediately suggests the potential relevance of the correction term. Solutions of the form 3 cosh(α 1 kh 2 ) + U 3 cosh(α 2 kh 2 )) exp( ik(x 1 − vt)), (6.6) satisfy governing equation (6.2) as long as 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: within which S i = sinh(α i kh 2 ), C i = cosh(α i kh 2 ) (6.9) and 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: (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 2 /h 3 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: For our new beam bending equation, this corresponds to the limit when q ≡ h 2 /h 3 1; hence, in view of expansion (5.1), we can say that 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 0 (ν, 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 4 ) 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.

Shear correction factor
In derivation of his beam equation, Timoshenko relied upon two key constitutive assumptions: 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) 1 encapsulates the same leading-order relationship between bending moment M and the radius of curvature (∂φ/∂x 1 ) −1 as in the Euler-Bernoulli theory. Assumption (7.1) 2 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 κ: (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) 2 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) 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) 2 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 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 2 k 2 /A. This procedure yields the following approximate dispersion 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 coth(qπ 2n−1 2 ) (2n−1) 5 , (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: 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 κ.  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.
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)]: 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) 2 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: coth(qπ 2n−1 2 ) (2n−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)]: 36q 4 π 6 (2m+1) 2 n 2 (4n 2 +(2m+1) 2 q 2 ) , (7.9) and by Yu & Hodges [14, eqn (68) 1 ]: 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 0 (ν, 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 (  dubious; however, rather than casting the doubt on definitions of κ NPK = κ SL , this simply reflects the fact that (7.1) 2 is not a consistent statement from the point of view of long-wave low-frequency asymptotics. 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 = √ 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 This identity is, of course, approximate, but it becomes more and more accurate as 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.

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 2 ≥ h 3 or, equivalently, that q ≡ h 2 /h 3 ≥ 1, then the lower bending mode of the beam is governed by in which I x 2 = 2h 2 (2h 3 ) 3 /12, A = 4h 2 h 3 and the non-dimensional coefficient k 0 (ν, q) is defined by equations (4.20) and (4.17) 2 . The governing equation for higher bending mode can be obtained by replacing I x 2 and q with I x 3 = (2h 2 ) 3 2h 3 /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: coth(qπ 2n−1 2 ) (2n−1) 5 .
( 8.2) We showed that equation (8.2), although apparently new, is equivalent to the expressions for κ previously derived in [21][22][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][37][38] and references therein). We are hoping to address the sparsity of similar results for higher-order beam theories in our future work. 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: coth(π (2k − 1)q/2) (2k − 1) 5 (A 1) = 1 π ∞ k=1 2 6 (2k − 1) 6  In the view of identity (A 7), the equivalence of expressions (7.5) and (7.6) can now be established by elementary algebraic manipulations.