Extension of Boley’s continuum mechanics-based successive approximation method to two-layer rectangular beams

An extension of Boley’s continuum mechanics-based successive approximation method is presented for rectangular beams composed of two isotropic linear elastic layers. The solution is cast into the form of tables, in complete analogy to the tables originally presented by Boley and Tolins for single-layer strips. The first column in these tables corresponds to the classical Bernoulli–Euler theory of beams. The further columns represent comparatively fast converging correction terms of an increasing refinement. Our two-layer formulation automatically satisfies the stress continuity conditions at the interface of the two layers. Enforcing displacement continuity at the interface between the layers, we derive results that do satisfy the equilibrium field equations, the stress continuity conditions at the interface and the stress boundary conditions at the upper and lower edges. When converged, the field constitutive relations and the displacement continuity at the interface between the two layers are also satisfied. We present a compact formulation, which allows writing down the results for more than the three successive steps considered by Boley and Tolins. The elasticity solutions presented subsequently can be used as novel analytic benchmarks for comparison with refined structural mechanics beam theories. Interior solutions for beams with a finite axial extent can be obtained by assigning approximate boundary conditions at the lateral ends. Comparisons to finite element computations for a clamped–clamped beam give strong evidence for the correctness of our analytic results.


Introduction
The continuum-mechanics-based theory of elasticity, see, e.g., Timoshenko and Goodier [28], serves as benchmark for the assessment of the accuracy of various structural mechanics beam theories of increasing complexity. The latter beam theories intend to approximate multi-dimensional elasticity solutions by one-dimensional solutions, stemming from the fundamental formulations of continuum mechanics, see Altenbach and Öchsner [5], by one-dimensional solutions, see, e.g., Krommer and Vetyukov [21]. Particularly, the plane stress formulation of two-dimensional linear elasticity ever since has been used for a comparison in case of plane bending and stretching of beams with rectangular cross sections. This strategy dates back to von Karman [16] and Seewald [26], who, using integral transform techniques, studied an isotropic homogeneous rectangular strip in plane stress according to the elasticity theory and compared their results to the classical Bernoulli-Euler (BE) theory of beams. Von Karman [16] eventually showed that the curvature of the axis of an infinitely extended strip can be expressed as an infinite series in the bending moment and its spatial derivatives. The second (non-vanishing) term in the curvature series, containing the second derivative of the bending moment, Communicated by Andreas Öchsner.
J. Gahleitner (B) · H. Irschik Institute of Technical Mechanics, JKU Linz, Austria E-mail: johannes.gahleitner@jku.at was explicitly communicated in [16]. Von Karman's work [16] was accompanied by a study of Seewald [26], who gave a corresponding series representation for the stresses across the strip. The first terms in the series representations of [16,26] for the curvature, the axial stress and the shear stress turned out to be equal to the formulations provided by the classical BE theory of beams. For smoothly varying external loadings and for slender beams, this explained the success of the BE theory in a relatively wide range of practical applications. The next terms in [16,26] represent corrections of the classical beam theory that are of an increasingly minor refinement, if the loadings, and thus the bending moments, are smooth enough. The work of von Karman [16] and Seewald [26] suggested developing step-by-step procedures for deriving these correction terms for various load cases. A review on such successive approximation techniques was provided by Ghugal and Shimpi [12]. A first step-by-step procedure for stresses was introduced in a brief note by Donnell [8], who studied transverse loadings that are continuously distributed at both the upper and the lower edges of rectangular beams. The corresponding computation of the beam displacements was also sketched in [8]. Independently, Boley and Tolins [7] used a convenient iterative procedure in order to derive, from the continuum mechanics-based two-dimensional elasticity theory, the stresses and deflections due to transverse and shear loadings that vary smoothly along the outer edges of rectangular beams. The iterative procedure used in [7] was derived previously for thermally loaded rectangular beams by Boley [6]. For an introduction into the iterative computation of stresses in rectangular beams, we refer to the book of Altenbach and Naumenko [4] where some instructive applications to cantilever beams are to be found also. In Boley's step-by-step technique, use is made of the properties of the Airy stress function, where the traction (stress) boundary conditions at the upper and lower faces of the beam are exactly satisfied. In [7], results for the first three steps were presented in the form of tables, the highly informative form of which represent a strong motivation for the present work. The tables given by Boley and Tolins [7] demonstrate how displacements and stresses within the rectangular beam do depend on bending moment and axial force, i.e., on integrals or derivatives of the latter, where a single order only appears in each of the steps of the successive computation of the stresses. From the latter, displacements are derived by integration the elastic constitutive relations. Again, only one order of derivative is kept in the columns of the tables presented in [7]. This means that, for the displacements, some results from one respective step of successive approximation were shifted tacitly to the next column in [7]. Since derivatives of bending moment and axial force are related to the known loading of the beam by ordinary differential equations via equilibrium considerations, the higher-order correction terms in the columns of the tables of [7] are known. For an application of Boley's method to the problem of piezoelectric rectangular beams, see Krommer and Irschik [20]. Duva and Simmonds [9] modified the expressions in the step-by-step expansions in [7,8] to account for orthotropic beams weak in shear and gave an explicit presentation of computable error estimates. The step-by-step formulations in [7,8,16,26] refer to beams that are infinitely extended in axial direction. They satisfy exactly the stress boundary conditions at the upper and lower edges of the rectangular beam, as well as the equilibrium relations within the beam, while constitutive relations are generally satisfied only approximately. For polynomial loadings of a finite order, only a limited number of terms in the step-by-step procedures is necessary to satisfy constitutive relations exactly. For beams of finite length, three boundary conditions can be assigned at each of the lateral ends in a more or less straightforward manner, e.g., normal force, shear force and bending moment at a free end or bending moment and two displacements of the axis point at a simply supported end. The resulting solutions do coincide with some analytic results for plane stress problems presented in Timoshenko and Goodier [28]. It was shown by Irschik [15] that the corresponding two-dimensional elasticity results for beams of finite length can be obtained conveniently from enhanced BE beam problems by analogy, since BE problems are also assigned with three boundary conditions at a lateral end. Enhancement suggested in [15] concerns fictitious, but known additional eigenstrain loadings, the latter corresponding to the correction terms in the above step-by-step solutions in [7,8]. However, the solutions with only three boundary conditions at each of the lateral beam ends generally cannot satisfy the boundary conditions that must be required to hold at each of the cross-sectional points on a lateral end, in order to obtain a completely exact two-dimensional elasticity solution. Solutions for beams of a finite length with only three selected boundary conditions at each lateral end therefore often are denoted as interior solutions. The difference between an interior solution and the completely exact two-dimensional elasticity solution should vanish away from the lateral ends according to St. Venant's principle, which, however, must be considered with certain care, see, e.g., Karp and Durban [17] and the literature cited there. For elastic interior solutions for orthotropic beams, see Tullini and Savoia [29], who provided a successive approximation technique involving Fredholm equations of the second kind and stated a sufficient condition for the convergence of their series expansion. An exact interior solution for a beam with span-wise constant transverse loading was presented by Karttunen and von Hertzen [18], see [19] for the anisotropic case. In extension on studies on homogeneous beams, an extensive literature has been published on problems of the deformation of functionally graded and layered or laminated beams. For reviews, see Irschik [14], and Altenbach [3], Ghugal and Shimpi [12], as well as Sayyad and Ghugal [23]. For advanced topics in laminated beams concern, e.g., deformable interfaces, see Schulze [25], Weps [30], Adam and Heuer [2], or for spatially varying interfaces see Li and Liu [22]. Concerning Boley's successive method, an anisotropic beam theory has been given recently by Gahleitner and Schöftner [11], and an application to bimorph piezoelectric beams of finite length composed of two layers with equal material parameters was recently presented by Schöftner and Benjeddou [24]. Nir Enuma [10] presented an extension of the works of von Karman [16], Seewald [26], Donnell [8] and Boley and Tolins [7] by developing a successive procedure for functionally graded orthotropic rectangular beams with transverse loading over the upper surface. A strategy for extending Boley's method to laminated rectangular beams was also sketched in [10], but no explicit analytic formulas were presented.
Although the more recent literature refers to advanced topics in composite beams, such as mentioned above, there is a lack of systematic extensions of Boley's successive method to layered rectangular beams. Successive methods would allow to obtain analytic solutions of increasing accuracy for the sake of comparison. It is the scope of the present contribution to extend the successive presentation given in Boley and Tolins [7] to the case of a composite rectangular beam of constant total height, made of two isotropic homogeneous layers with different material parameters and heights, the two layers being connected rigidly by a common straight interface. This problem has the advantage of including the original single-layer problem treated in [7,8,16,28] as a special case. As simple as the problem may appear, the analytic results for the two-layer beam turn out to be quite complex. Therefore, from space limitations, we subsequently restrict to the corresponding extension of Case I of [7], where results for a single-layer beam under transverse loading applied at the upper edge of the beam were considered. Our intent is to enable formulations for the two-layer beam in a compact manner, in the form of tables, analogous to the ones given in [7]. The method presented below, however, can be easily applied to the other two load cases treated in [7] by analogy. Our subsequent methodology is different from the one in [10], since we satisfy the stress continuity conditions at the interface between the two layers from the onset. We also take the opportunity of clarifying some mathematical aspects of the Boley method of successive approximation, which were not explicitly addressed in the above cited literature. Particularly, we explain the consequences of the above-mentioned shifting process of bringing the results into a form, in which only one order of derivative appears in a single step, i.e., in each of the columns of the tables in [7]. Moreover, as a main result, we present formal mathematical expressions, which allow to compute the results for any number of considered successive steps, enabling to obtain complete analytic results for any order of polynomial loading. This is not directly possible from the results presented originally, e.g., in Boley and Tolins [7], where results for the first three steps were presented, which, however, do not allow extensions by complete induction.
Our paper is organized as follows: Notations and basic relations of our problem are stated in Sect. 2. Boley's successive approximation technique is shortly discussed in Sect. 3. Extended Boley-type tables for the single layer of the composite two-layer beam, which are needed subsequently, are provided in Sect. 4. A strategy for extending Boley's method to two-layer beams is given in Sect. 5. Displacement and stress continuity conditions at the interface are treated, eventually leading to a form that is indeed in complete analogy to the results for a single-layer beam treated in Table I of Boley and Tolins [7], namely in dependence from the bending moment and from the axial force, as well as from their respective integrals and derivatives. A Boleytype table representation for considering three successive steps is presented for the two-layer beam, as well as more general formulas allowing the consideration of any number of steps, see Table 1 and Eqs. (60)-(62) for the Airy stress functions. In Sect. 6, emphasis is given to the computation of displacement fields and the necessary shifting process for obtaining the desired successive form, see Table 1 for three successive steps, and the more general representation in Eqs. (74) and (75) for formal displacement representations that are applicable to any number of steps. Completeness of the successive representation for a given polynomial order of the applied loading is discussed. Also, for the sake of comparison with numerical computations, specific entries to Table 1 are presented for some special ratios between the elastic parameters, as well as between the thicknesses of the layers. In Sect. 7, using these specific entries, comparisons to numerical computations are presented, where a comparatively thick clamped-clamped beam of finite length with a linearly varying applied loading is used. For this case, the three steps considered in Table 1 give a complete analytic set of solutions for the infinitely extended beam. For the beam of finite length, analytic results are obtained by imposing additionally three advantageous boundary conditions for a clamped end that were suggested by Szabo [27]. Since these three boundary conditions are approximate, displacement distributions are obtained at the beam ends, which do not vanish completely there, although they satisfy Szabo's three boundary conditions. Numerical comparison is performed utilizing the computer-code ABAQUS [1], on the basis of a sufficiently fine net of 1000 Elements. A x y q(x) Fig. 1 Composite beam made of two layers under transverse loading q(x) first set of comparisons is obtained by prescribing the non-vanishing end displacements of the analytic solution to the finite element model at the beam ends. In this case, the exterior solution should vanish, and the analytic solution should completely coincide with the numerical solution. It is found that this is indeed the case, for both stresses and displacements. A second comparison is performed by prescribing vanishing displacements at the beam ends to the ABAQUS model. Then, an exterior solution should represent the differences between the ABAQUS solution and the analytic solution, the latter one representing the interior solution. These differences can be expected to fade away from the beam ends, due to the principle of St. Venant, see, e.g., [17] and Ziegler [31]. This also is nicely reflected by the presented results. These successful comparisons give strong evidence for the appropriateness of our Boley-type analytic solutions for the two-layer composite beam.

Problem statement
Consider an initially straight rectangular beam under the action of a distributed transverse loading that is applied at the upper edge of the beam. The beam is composed of two rectangular, homogeneous, isotropic and linear elastic layers, which are firmly bonded 1 together at their common interface, see Fig. 1.
Young's modulus, Poisson's ratio and thickness of a layer are denoted by E, ν and c, respectively. These entities are assumed to be different in the two layers. Index l refers to the lower layer, while index u corresponds to the upper one. Numbers α, β, γ are introduced in order to denote the ratios between the values of E, ν and c in the two layers: In the following, the thickness of the upper layer is abbreviated with c u = c, and we write shortly E u = E, ν u = ν. The two layers have an equal width b, see Fig. 1, which is used as unit reference length, b = 1 in the following. It is more over assumed that the thickness of each layer is sufficiently larger than b, i.e., that the rectangular cross section is sufficiently narrow, such that a state of plane stress can be taken into account. A global Cartesian (x, y)-coordinate system is attached to the interface of the layers, where the axial coordinate is denoted by x, and y is the transverse coordinate. The transverse loading is denoted as q(x); it is applied at y = c, see Fig. 1, and it is assumed to be sufficiently smooth. Suppressing the index that indicates a specific layer, the plane stress field equations in each of the layers read, see, e.g., Timoshenko and Goodier [28]: Normal stresses are denoted as σ x and σ y , and τ xy stands for the shear stress; axial and transverse displacements are denoted by u and v, respectively. The stress boundary conditions at the upper and lower edges of the beam are, see Fig. 1:

Boley's successive approximation technique
Our strategy for solving the above stated problem is to apply and to extend a methodology originally developed by Bruno A. Boley for single-layer rectangular beams, Boley and Tolins [7], see Boley [6] also. Boley's methodology consists in expanding Airy's stress function ϕ in the form of a step-by-step series solution: This formulation ensures that ϕ satisfies the bi-harmonic differential equation, as this should be, see Timoshenko and Goodier [28]. Note that Eq. (7) are to be understood as ordinary differential equations in the coordinate y.
The nonzero stress boundary conditions at the upper and lower edges of the single layer are satisfied by means of the first term, ϕ 1 . The stress boundary conditions for ϕ i , i > 1, thus become trivial at the upper and lower edges, which means that there is Having solved these one-dimensional boundary value problems step by step, the stresses follow as, see [28]: Airy stress functions subsequently will be shortly denoted as potential functions or potentials. From the stresses, displacements are eventually computed by integrating the constitutive relations, Eq. (3), and by adjusting the arbitrary functions which arise so as to satisfy the third of Eq. (3). In [7], stresses and displacements were computed step-wise for three loading cases, and the results were presented in the form of three tables, each showing the respective results of the first three steps of the procedure. The loading cases considered in [7] were a transverse loading at the upper edge of the single-layer beam (Table I of [7]), and a symmetric as well as an anti-metric shear loading applied at the upper and lower edges (Tables II and III of [7]). Some short explanations seem to be in order. The tables presented in [7] are formulated in terms of bending moments and normal forces and their integrals and derivatives. As already mentioned, for the displacements a single order of derivative only was used also in each of the columns of the tables in [7], the rest of the results of the respective steps being shifted tacitly to the next columns. The first columns in the tables by Boley and Tolins [7] correspond to ϕ 1 in Eq. (7) and coincide with the results of the elementary Bernoulli-Euler beam theory, see, e.g., Ziegler [31] for the latter theory. The second column, corresponding to ϕ 2 in Eq. (7), dates back to fundamental contributions by Karman [16] and Seewald [26], who obtained their results via an integral transform technique. The third columns in the tables of [7], corresponding to the third step in Eq. (7), were presented for the first time in [7]. Nowadays, results of a higher-order multi-step procedure can be easily obtained by means of symbolic computation. It is the scope of our present paper to extend the Boley-type integrate technique to two-layer beams.

Extended Boley-type tables for the single layer of the composite two-layer beam
In a first step, the composite beam presented in Fig. 1 is cut along the common interface of the two layers into two single-layer beams, in y = 0. The resulting free-body diagrams are shown in Figs. 2 and 3, which result in five loading cases, three for the upper layer and two for the lower one. In Fig. 2, the interfacial stresses are denoted by σ = σ y (x, 0) and τ = τ xy (x, 0), respectively. Note that the interfacial stress continuity conditions at y = 0 are satisfied automatically by this convention. Results for the Airy stress functions and the displacements for the five loading cases are gathered in abbreviated Boley-type tables, which are given in "Appendix A." Captions for these tables are denoted as follows: Table 2 belongs to  Fig. 2a, and Table 3 is due to the loading case in Fig. 2b. Two Tables, 4 and 5, are needed for the shear loading case in Fig. 2c. These follow by representing the τ -loading case as the sum of a symmetric and an antimetric loading, likewise to Case II and Case III in [7]. Correspondingly, the lower layer, see Fig. 2d, e, is represented in Tables 6, 7 and 8, respectively. The potentials and displacements presented in the tables of "Appendix A" are to be understood for the single layer of the composite beam, see Fig. 3 for positive directions of bending moments M, shear forces Q and normal forces P. Note that bending moments in "Appendix A" are to be understood with reference to the respective axis of the layer.
Bending moments, shear and axial forces are assigned with additional indices in order to distinguish between the load cases indicated in Fig. 2; e.g., M τ u and P τ u denote bending moments and normal forces due to the loading case τ in the upper layer, see Fig. 2c and Tables 4 and 5, respectively. The tables of "Appendix A" have been independently derived using the step-wise original method of Boley and Tolins [7], but they are formulated already in the global interfacial coordinate system shown in Fig. 1 Fig. 3 Stress resultants in the two layers of a differential beam element introduced there. The following convention for abbreviating x-derivatives is used: Particularly, we have for the upper layer Analogously, there is for the lower layer The following layer-wise differential equilibrium relations do hold, see Fig. 3 The tables of "Appendix A" can be used for single-layer homogeneous beams, in addition to the ones that were presented before by Boley and Tolins [7]. As in the latter reference, displacements in these tables are to be understood modulo rigid-body displacements.

A strategy for extending Boley's method to two-layer beams
In the Boley-type tables of "Appendix A," the expressions for potential functions and displacements in the two layers are expressed as functions of the respective internal forces and moments and their axial integrals and derivatives. Continuity of the interfacial stress components σ and τ , see Fig. 2, is automatically taken into account. What remains is to assure continuity of the displacements of the two layers at the common interface at y = 0, The index 0 refers to the interface of the two firmly bonded layers. In order to reduce the unknowns in the tables of "Appendix A," namely M qu , M σ u , M τ u , P τ u and their counterparts in the lower layer, M σ l , M τl , P τl , it turns out that it is sufficient to study the difference at y = 0 of the first derivatives of the axial displacements, u 0 , and of the second derivatives of the transverse displacements, v (2) 0 , which also must vanish: For instance, u 0u1 follows as the sum of the first derivatives of the axial displacements in the first columns of Tables 2, 3, 4 and 5 . This also holds for the next columns in the latter tables.
In a first step, we replace the unknowns M qu , M σ u , M τ u , M σ l , M τl and P τl , by four quantities, which we choose as M, M u0 , P and P τ u . This can be performed by using the following relations of static equivalence: The bending moment M u0 denotes the moment of the axial stress σ x in the upper layer with respect to the interface y = 0. M l0 is defined analogously. The total bending moment with respect to y = 0, and the resulting normal force P becomes: Note that we suppress the subscript 0 for M in the following. The subsequent differential relations are also needed: Equation (21) We eventually obtain the following expressions: Recall that in the original tables of Boley and Tolins [7], all entities were expressed by the bending moment M and the normal force P of the single-layer beam, as well as by their integrals and differentials. In order to approach Boley-type tables for the two-layer composite beam, it is straightforward to further express M u0 and P τ u by a series in the bending moment M, its even-order derivatives and the normal force P. Hence, we write Recall that P is taken as constant, see Eq.
In order to find expressions for the unknown coefficients a and b in Eqs. (33)-(38), we re-substitute Eqs. (33)-(38) into Eq. (18). Since the solution must hold for arbitrary functions M, we require that Eq. (18) must be satisfied for P and for each order of derivative of M separately. This yields the following set of equations gathered in Eqs. (39)-(46), which allows to compute the coefficients a and b in a step-by-step manner: One obtains from the first four relations, see Eqs. (39)-(42): This results for a 0 and b 0 can be substituted into the next two relations, Eqs. (43) and (44), in order to compute a 2 and b 2 , etc. The expressions for a 2 , b 2 , a 4 and b 4 are already lengthy, but can be easily obtained via symbolic computation. Now that M u0 and P τ u have been obtained according to Eqs. (31) and (32), the potentials in the upper and lower layers of the composite beam can be computed. We exemplary show this for the first part ϕ u1 of the upper layer's potential. The corresponding transversal stress σ yu1 is computed from the potentials of the first columns of Tables 2, 3, 4 and 5 by using Eq. (9): Substituting derivatives of Eqs. (19), (22), (23) and (24) into Eq. (51), σ yu1 is expressed via M (2) , M u0 and P (2) τ u only: Following an argument by Seewald [26], the functional dependency on the y-coordinate in the potentials and in σ y , which is the second axial derivative of the former, see Eq. (9), can be taken as equal. Hence, we find: Analogously, in order to compute ϕ u2 and ϕ u3 , the relations in Eqs. (19)-(24) as well as their derivatives with respect to x are substituted into the potentials of the second and third column of Tables 2, 3, 4 and 5 . The same procedure has to be applied to the lower layer to get ϕ l . Substituting M u0 and P τ u from Eqs. (31) and (32) and resorting with respect to M, P, M (2) and M (4) , we get the counterparts for the two-layer beam of Case I of the original single-layer table in Boley and Tolins [7]. In order to demonstrate the procedure, we note the result for the first term ϕ u1 and ϕ l1 in more detail: 3 (55) For a compact notation, a dimensionless parameter λ has been introduced: It is noted that Eλc 3 can be interpreted as an effective flexural rigidity, reducing to E2c 3 /3 for the homogeneous beam of thickness 2c. Equations (54) and (55) now transform to: For instance, the following abbreviations have been introduced in Eq. (57): Repeating the procedure for the next successive steps and applying the method of complete induction, the potentials ϕ u and ϕ l can be represented in the following compact series form, which allows to consider any number n of steps in the successive procedure: with i as the number of iterations and f , g as polynomials in y: The coefficients for f and g can be taken from Tables 9 and 10 in "Appendix B," where k i j ,k 1 j correspond to the upper layer and g i j ,ḡ 1 j to the lower layer. Now, using Eqs. (60) and (61) and restricting to the first three steps lead to a Boley-type table for the two-layer composite beam, see Table 1. Together with the formal series representation in Eqs. (60) and (61), Table 1 represents the main results of our paper.  (20). For the sake of comparison to numerical computations, coefficients corresponding to a two-layer composite beam with α = 4, β = 1, γ = 3 are presented in Tables 11 and 12 in "Appendix C." We now add a brief note on criteria to truncate the number of iterations. For polynomial loadings of order n ≥ 0, the number of iterations is necessary for an exact solution which follows from Table 1. Furthermore, the magnitude of the correction terms for i > 1 depends on the thickness-to-length ratio (aspect ratio) μ. This is exemplarily shown subsequently for the non-dimensional potential for a homogeneous beam (Case I of [7]) with reference bending moment M, non-dimensional transverse coordinate η = y/c, aspect ratio μ = 2cβ/L < 1, beam length L and parameter β of Eq. (1). Since the methodology is comparatively fast converging, a three-step procedure often is accurate enough, dependent on the aspect ratio, which is exemplarily shown for the upper fiber normal stress σ x (x, c) of a homogeneous beam in Fig. 4. The results are normalized to the elementary (first iteration) normal stress. The exact solution is a polynomial of 11th order which approximates a sine half-wave and leads to i = 7 necessary iterations, see Eq. (63). For non-polynomial loads, such as sine and cosine, the half-wave number must be taken into account in order to achieve convergence of the iteratively obtained results. Exemplarily, this is shown for a simply supported, homogeneous beam loaded by a distributed load q(x) = q 0 sin(mπ x/L) with the integer odd half-wave number m ≥ 1, see Fig. 5. For this special case of loading and supporting conditions, a closed-form solution can be determined with the Ansatz φ(x, y) = M(x) f (y) where the bending moment is computed from M (2) = −q(x) satisfying boundary conditions at the lateral ends. Substituting the Ansatz into the bipotential equation leads to an ordinary differential equation of fourth order for f (y). Four constants of f (y) are to be adjusted so that the stress boundary conditions at the horizontal surfaces are satisfied, see Fig. 5. To demonstrate the convergence behavior of the introduced iterative procedure, the iteratively determined entities are compared to those obtained from the closed-form solution. Exemplarily, this is shown for the upper fiber normal stress σ x (c, L/2) evaluated at x = L/2 where the closed-form solution reads including the parameters μ and m. Convergence studies are presented for two half-wave numbers (m = 1, m = 5) and for two aspect ratios (μ = 1/5, μ = 1/4), see Fig. 6. In Fig. 6a, the iterative and closed-form results agree already at the third and fourth iteration i for a sine half-wave with m = 1. For a half-wave number m = 5, the iteratively obtained results agree for i = 7 and i = 8, see Fig. 6b. This gives evidence for the statement by Boley and Tolins in their Paper [7] concerning the limited usage of the method for non-smoothly varying forces.

Computation of the displacement fields
We demonstrate the displacement-computation procedure for the upper layer in more detail, since it performs analogously for the lower layer. The displacements are computed according to Eq. (3), where the normal stresses σ x and σ y from Table 1 are substituted. By integration, one finds for the axial displacement: and for the transversal displacement: where roman superscripts denote derivatives with respect to y. In the course of the integration procedure, two unknown functions, G(y) and F(x), are to be introduced. The latter functions have to be determined so that the third of Eq. (3) is satisfied. Derivatives of Eq. (66) with respect to y, and of Eq. (67) with respect to x, as well as the shear stress τ xy from Table 1  Substituting the coefficients according to "Appendix B" shows that the terms containing Mdx, cPx, M and M (3) are functions of x only. However, the terms containing M (5) and M (7) are y dependent also. To satisfy the requirement that Eq. (68) must contain functions of x and y only, but no mixed terms, we assume that M (5) = const for the time being. One then finds the following ordinary differential equation to compute G(y): In the framework of this reasoning, we obtain for F(x): The displacements are determined by integration of Eqs. (69) and (70) and re-substituting the results into Eqs. (66) and (67). Modulo rigid-body displacements, this leads to The so-determined displacements do satisfy the constitutive relations in Eq. (3). With respect to the displacement continuity conditions, see Eq. (17), we find from the coefficients of Tables 9 and 10 in "Appendix B" that Thus, Eq. (17) are satisfied with M (4) = −q (2) = 0, or for the homogeneous case with α = β = γ = 1. First, some short remarks concerning the homogeneous case seem to be instructive. In order that the interface displacement continuity conditions are satisfied, M (4) and M (5) need not to vanish here. In order to include corresponding non-vanishing terms into a Boley-type table, an additional column n = 4 must be opened, to which the terms containing M (4) and M (5) have to be shifted. But additional terms containing M (4) and M (5) , as well as derivatives of M of higher order, come into the play from integrating the stress contributions in the newly opened columns, and these terms must be included in the table also. So, it is to be noted that the displacement terms in the ith column of a Boley-type table, say, do not only stem from the stresses in this ith column, but also include contributions from column i −1. We have designed our present notation by introducing functions f i in order to make this fact more visible and tractable than this is possible from the original tables in Boley and Tolins [7]. The above remarks on the shifting procedure of course apply to the composite two-layer beam also, c.f. the columns in our Table 1. However, for this case of an inhomogeneous beam, the situation is complicated in so far, as it has turned out above to be necessary to assume that M (4) = −q (2) = 0, in order that the displacement continuity conditions at the interface y = 0 are satisfied in the course of a three-step procedure, n = 3, see Eq. (73). Now, in order that complete solutions can be obtained for polynomial loadings q of a higher order than one, the situation again turns out to be resolvable by opening additional columns and repeating the above procedure. Analogous to our above discussion for n = 3, the vanishing of derivatives of the loading of a properly increased order has to be required, in order to satisfy the respective extended forms of Eq. (73). Following the above reasoning, the displacements presented in Table 1 have been eventually obtained, restricting to the first three columns, n = 3; likewise this was done in Boley and Tolins [7]. Note that our Table 1 reduces to Case I of [7] for the homogeneous case, α = β = γ = 1. However, as already shortly noted above, our present notation allows to formulate the displacements of the composite beam in the following series form, without restricting the number n of columns, with n > 1. By complete induction, the results in Table 1 for the displacements can thus be written as: Mdx 2 + cP f (I I I )  The following convenient definition for a generalized cross-sectional rotation, suggested by Szabo [27], is used as trivial boundary condition in order to obtain the analytic solutions according to the above formulations: The FEA is performed by setting up two different types of boundary conditions, which are denoted as "elastic" and "rigid" in the following, and are abbreviated by the superscripts "E" and "R," respectively. For the elastic case "E," the set of boundary conditions reads at x = ±L/2: where u and v mean the analytic solutions derived by our above theory. Hence, since in our case of a linearly distributed loading, see Fig. 7, the three-step formulation in Table 1 gives a complete solution, the FEA for case "E" should coincide with our analytic solution.
For the rigid case "R," the following trivial constraints are prescribed at x = ±L/2 for the FEA: This corresponds to the model of a completely undeformable clamped end. In order to obtain the analytic solution according to Table 1, taking into account the reduced clamping boundary conditions in Eq. (76), the corresponding bending moment M and normal force P are needed. They are obtained by integration of the beam equilibrium relations, Eq. (21), with the three redundant static integration constants. We now take into account the displacements in Table 1 and substitute the coefficients of Tables 11  and 12, which refer to our present example with α = 4, β = 1, γ = 3. Adding rigid-body displacement fields with three further kinematic unknown constants, satisfaction of the six boundary conditions in Eq. (76) leads to a set of six linear relations for the three static redundant constants and the three kinematic constants. Bending moment and normal force eventually follow in analytic form to:  For determining the analytic stress solutions, we substitute bending moment M and normal force P into Table 1 by using the coefficients of Tables 11 and 12 in "Appendix C." Corresponding non-dimensionalized stress distributions are shown in Fig. 8 where the stresses σ x , σ y and τ xy are evaluated at three axial coordinates: The stresses are normalized by introducing the non-dimensional transverse coordinate η = y/c in the following form: where the terms in the denominators of Eq. (81) denote the results from our theory evaluated at the respective x-coordinate, where σ * x is evaluated at the lower fiber, at η = −1, σ * y is evaluated at the upper fiber, at η = 1, and τ * xy is taken at the common interface, in η = 0, respectively. It turns out that the stress distributions derived from our theory are in full agreement with the FEA results for case "E," as this should be, since in this case the exterior solution does vanish, because boundary conditions do coincide, see Eq. (76). For rigid constraints, case "R," the FEA solution deviates from our analytical solution, since boundary conditions are different, see Eq. (79), and case "R" thus takes into account end effects, which are not covered in this work. It is seen from Fig. 8, however, that, in accordance with Saint-Venant's principle, the non-vanishing exterior solution, which represents the difference with respect to our analytic solution, decays nicely here within a short distance away from the clamped ends, also in the present case of a comparatively thick composite beam. At the clamped ends, the exterior solution obviously shows stress concentrations and a jump-type behavior, which cannot be well reflected by the FEA. That is why we have performed our comparison for a cross section in the vicinity of the clamped end only, at x = −49L/100, see Fig. 8a.

Displacements
Cross-sectional deformations in axial and transverse direction are presented in Fig. 9. The deformations u and v are normalized to   Table 1 only are shown as red lines where u * and v * denote the displacements, derived from our theory evaluated at the upper fiber at η = 1 and at the respective x-coordinates. Analogous to the stress fields presented in Fig. 8, the analytical results are in excellent agreement to the FEA results for case "E," i.e., when elastic constraints are prescribed on the FEA model. Near the clamped end, the analytical solution and the FEA-rigid solution "R" again deviate from the analytic solution due to end effects. Away from the clamped end, see Fig. 9b, c, the transverse displacements do coincide, and the shape of the axial displacements is similar to the FEA-rigid result "R," but shifts are to be observed. The deflection cure v(ξ, 0) of the axis and the axial cross-sectional mid-span displacement u(ξ, 0) are presented in Fig. 10 by introducing the non-dimensional axial coordinate ξ = x/L. The entities are normalized to with u * * = u(1/4, 0) = 0.00086 mm and v * * = v(0, 0) = 0.033 mm as the analytical results for the displacements at ξ = 0 and η = 0 exemplarily. Likewise to the latter results, the FEA results for case "E" are in agreement with the analytic solution. This also holds for the deflection curve of the "R" solution. Additionally, the results obtained from the first column in Table 1 are shown in Fig. 10. A deviation to our analytic solution of 37% can be observed for the deflection at the mid-span ξ = 0.

Conclusion
In the present contribution, analytic plane stress linear elasticity results for two-layer composite rectangular beams of infinite axial extension have been obtained, in extension and in analogy to Boley's successive approximation procedure, see Boley and Tolins [7] for the single-layer case. To the best knowledge of the present authors, our results for the two-layer beam have not been published elsewhere in the open literature so far. Although results have turned out to be quite complex, a compact notation has allowed us presenting them in a form analogous to the one stated by Boley and Tolins [7] for the single-layer beam, see Table 1, which contains Case I of [7] as a special case. Our two-layer formulation automatically satisfies the stress continuity conditions at the interface of the two layers. In order to enable this, we above have started from successive approximation results for appropriate load cases that are applied to the single layer separately and that are put together properly. Eventually enforcing displacement continuity at the interface between the two layers, we have succeeded to derive results for the composite beam that exactly do satisfy the equilibrium field equations, the stress continuity conditions at the interface and the stress boundary conditions at the upper and lower edges of the beam. Moreover, our compact notation allows considering more than the three successive steps that were considered in [7], see Eqs. (60)-(62) and Eqs. (74) and (75) for stress functions and displacements, respectively. In extension of previous successive approximation formulations from the literature, it is thus possible to obtain complete solutions for beams for imposed loadings of an arbitrary polynomial order, meaning that the solutions do satisfy exactly the stress boundary conditions at the upper and lower edges of the composite beam, the stress transition conditions at the interface of the two layers, as well as the equilibrium relations within the layers. When the number of steps needed for the respective polynomial order of the imposed load is taken into account, the displacement continuity relations at the interface as well as the linear constitutive relations within the two layers are also satisfied exactly by our results. However, Boley's successive approximation procedure is known to converge rather fast. Boley and Tolins [7] restricted their presentation on single-layer beams to three steps. We have presented corresponding three-step results for the two-layer beam in more detail in the Boley-type Table 1. Entries for specific ratios between the elastic parameters and between the thicknesses of the layers have also been given, for the sake of subsequent comparison to finite element computations. Interior solutions for beams of finite length have been derived in a straightforward manner. Particularly, we have presented corresponding analytic results for a comparatively thick clamped-clamped beam under the action of a linearly varying applied load, for which Table 1 is sufficient. As can be seen from Figs. 8, 9 and 10, the results of our proposed method agree well with the results of the finite element solutions. As far as the present authors can see, there are nevertheless two main advantages of the proposed analytical formulation: First, from a theory point of view, the method allows a comparison of the present semi-exact results for displacements and stresses with the assumptions of numerous refined beam theories that have been developed in literature, the basic Bernoulli-Euler theory being contained directly in the first column of Table 1, as mentioned above. For analogies between the Bernoulli-Euler theory and some refined beam theories, see, e.g., Irschik [13]. Table 1 allows to relate the appropriateness of assumptions of refined theories to the polynomial order of the loading. Moreover, from the point of engineering practice, formulations of the present type allow to perform parametric symbolic computations in the design process of composite beams. It is hoped that our results will serve as further benchmark solutions for performing comparisons to higher-order theories for composite rectangular beams, and that our above discussion on some mathematical aspects of Boley's successive approximation procedure, such as on the shifting process, will stimulate extensions of this method to more advanced topics, such as non-isotropic and non-homogeneous layers, or non-rigid interfaces. Table 7 Load case τ l M , c l = βc Table 8 Load case τ l P , c l = βc (2) τl Appendix B Table 9 Upper layer coefficients k i j andk 1 j k 1 j k 2 j k 3 jk1 j 0 (a 0 − 1)λ a 2 λ a 4 λ -