Flexural vibrations of geometrically nonlinear composite beams with interlayer slip

This paper addresses moderately large vibrations of immovably supported three-layer composite beams. The layers of these structural members are elastically bonded, and as such, subjected to interlayer slip when excited. To capture the moderately large response, in the structural model a nonlinear axial strain-displacement relation is implemented. The Euler–Bernoulli kinematic assumptions are applied layerwise, and a linear interlaminar slip law is utilized. Accuracy and efﬁciency of the resulting nonlinear beam theory are validated by selective comparative plane stress ﬁnite element calculations. The outcomes of application examples demonstrate the grave effect of interlayer slip on the geometrically nonlinear dynamic response characteristic of layered beams.

The dynamic behavior of composite beams with interlayer slip has been studied less frequently. In the pioneering paper of Girhammar and Pan [11], lateral vibrations of layered beams with interlayer slip were analyzed. Later, the same group provided an extension and generalization of this theory for the exact dynamic analysis of such structural members [13]. Adam et al. [3] developed an accelerated modal series solution of vibrations of flexibly bonded layered beams. Piezoelectric and thermo-piezoelectric vibrations of beams with interlayer slip were studied in [16,17]. Challamel [5] treated lateral torsional vibrations of composite beams with partial interaction. The random dynamic response of an uncertain compound bridge subjected to moving loads was analyzed in [4]. Recently, Di Lorenzo et al. [22] presented an efficient procedure for the computation of lateral vibrations of discontinuous layered elastically bonded and non-classically supported beams.
A moderately large response in the presence of immovably supports causes nonlinear axial strains, i.e., the beam response becomes geometric nonlinear. In various studies, the influence of the membrane stress due to stretching of the central fiber on the dynamic response of homogeneous (e.g., [9,31]) and layered beams (e.g., [15,20]) has been considered. However, to the best knowledge of the authors, moderately large vibrations of beams with interlayer slip have not been analyzed yet. To fill this gap, in this contribution a theory for the analysis of the dynamic flexural response of composite beams with elastically bonded layers is presented, whose supports are rigidly held apart. In this theory, nonlinear axial strain-displacement relations are considered, which originate from moderately large vibrations of the member with horizontally restrained supports. In a common assumption, for each layer, the Euler-Bernoulli assumptions are assigned separately. To keep the derivations simple and clear, a three-layer composite beam with symmetrically disposed layers is considered. However, it should be noted that the proposed theory can be extended to beams with an arbitrary number of asymmetrically disposed layers.
The following paper is structured as follows. At first, the governing kinematic equations and their relation to the layerwise cross-sectional resultants of a three-layer composite beam are established. Conservation of momentum, in combination with the cross-sectional resultants expressed in terms of the kinematic variables, yields a set of coupled geometrically nonlinear equations of motion. Corresponding boundary conditions for various horizontally restrained beam ends are formulated. The proposed solution procedure is based on a modal expansion of the lateral deflection into the first few mode shapes of the corresponding linear beam. In an illustrative example, the transient response of a soft-hinged immovably supported member to a half-wave sinusoidally distributed time harmonic excitation is analyzed. In order to verify the proposed beam theory, the derived response is compared to the outcomes of an elaborate plane stress finite element analysis. The second example studies the nonlinear resonance of a beam with partial interaction, to show the effect of different load amplitudes and varying interlaminar stiffness on the amplitude frequency response of various kinematic variables and stress resultants.

Basic equations
A single-span beam of length l in principal bending about the out of plane (y-)axis, composed of three elastically bonded layers with constant rectangular cross section along the central beam (x-)axis, is considered, as for instance shown in Fig. 1. The geometry and material parameters of the external layers are the same. That is, with h i denoting the thickness of the ith layer, A i is the ith layer cross-sectional area, E i the ith layer Young's modulus, ρ i the ith layer density, E A i the ith layer extensional stiffness, and E J i the ith layer bending stiffness. The subscript i = 1 refers to quantities of the top layer, i = 2 to the central layer, and i = 3 to the bottom layer. A time-varying distributed lateral load q(x, t) excites the structural member to flexural vibrations. Since the bond between the layers is flexible, the layers displace against each other when the beam is deflected. This relative displacement between the top layer and the central layer, Δu 12 , and the central layer and the bottom layer, Δu 23 , is referred to as interlayer slip. Assuming that the layers are rigid in shear, Euler-Bernoulli theory is applied to each layer separately. This yields the lateral deflection w i and horizontal displacement u i at distance z i from the neutral axis of the ith layer as [16] with (.) ,x indicating partial differentiation with respect to x, and u (0) i the ith axial displacement at z i = 0 (see Fig. 2). Expressing the axial displacements of the outer layers, u (0) 1 and u (0) 3 , in terms of the central axial  displacement u (0) 2 , cross-sectional rotation w ,x , and the interlayer slips Δu 12 and Δu 23 , leads to [16] In these relations, d is the distance between the central axis and the neutral axis of the top/bottom layer, as shown in Fig. 2. In case of rectangular outer layers d = (h 1 + h 2 )/2. Moderately large lateral vibrations in the presence of immovable supports strain the central axis, and thus, the axial strain-displacement relation becomes nonlinear (see e.g., [32]), The longitudinal strain at any fiber of the beam is therefore Assuming a constant slip modulus K , which is the same for both interfaces, the interlaminar shear traction t s12 (t s23 ) is linearly related to the interlayer slip Δu 12 (Δu 23 ), and in combination with the kinematic relations of Eq. (3) it follows that Application of Hooke's law delivers the relation between the layerwise bending moment, M i , and the layerwise axial force, N i , respectively, and the kinematic quantities (see e.g., [32]), The global stress resultants for the entire cross section (compare with Fig. 3) are composed of the layerwise quantities according to No external axial load is applied to the beam, and thus, the global axial force N (also referred to as membrane force) emerges from the nonlinear axial strains due to moderately large deflection w only. The effect of N on the dynamic response is captured through second-order analysis. That is, conservation of momentum in axial (x-) and transverse (z-)direction and conservation of angular momentum about the y-axis is applied to an infinitesimal beam element in its deformed state (shown in Fig. 3), yielding where μ = 2ρ 1 A 1 + ρ 2 A 2 denotes the mass per unit length,ẅ is the lateral acceleration, and T is the transverse cross-sectional force. In Eq. (11) the effect of longitudinal inertia and in Eq. (13) the effect of rotatory inertia have been omitted, thus limiting the analysis to the lower frequency range. Furthermore, in a common assumption of second-order analysis, the longitudinal force S has been replaced by the axial force N [27].
Equation (13) is differentiated with respect to x and subsequently combined with Eq. (12). Considering that according to Eq. (11) the axial force N is only a function of time t (and not of x) this results in In x-direction, omitting the longitudinal inertial, layerwise application of conservation of momentum yields [16] 3 Governing equations

Equations of motion
The equations of motion of this beam problem are obtained by expressing Eqs. (11), (14)- (17) in terms of the governing kinematic variables, i.e., lateral deflection w, interlayer slips Δu 12 and Δu 23 , and axial displacement of the central axis u which results from combining Eq. (10) with constitutive relation Eq. (8) and kinematic relation Eq. (3), into Eq. (11) leads to the first governing equation: Next, in Eq. (9), M i , N 1 , and N 3 are substituted by Eqs. (7) and (8), respectively, and subsequently, variables u 3,x are replaced by the derivative of Eq. (3) with respect to x. This yields the total bending moment as with denoting the bending stiffness of the rigidly bonded beam, and the bending stiffness of the non-composite beam, i.e., K = 0. The second derivative of Eq. (21) with respect to x inserted into Eq. (14) leads to the second equation of motion, Now, Eqs.
One of the five coupled nonlinear partial differential equations of motion, Eqs. (20), (24)- (27) in terms of the four governing kinematic variables w, Δu 12 , Δu 23 , and u (0) 2 are redundant because the sum of the underlying Eqs. (15) to (17) yields Eq. (11). Consequently, Eqs. (20), (25), and (26) are condensed to two equations, which are easier to solve. To this end, Eqs. (25) and (26) are summed up with the outcome Then, Eq. (25) is subtracted from Eq. (26). In the resulting expression, quantity u 2,x x + w ,x w ,x x is eliminated by means of Eq. (20), yielding The parameter is proportional to K , and hence, can be considered as an indicator of the degree of composite action in longitudinal direction.

Boundary conditions
The equations of motion are solved in combination with the actual boundary conditions. Because this mixed initial-boundary value problem comprises one fourth-order differential equation (Eq. (24)) and three secondorder differential equations (Eqs. (20), (28), (29)), for its solution five boundary conditions are to be specified at each end. The ends of the considered structural members are immovably supported, i.e., at both ends the horizontal displacement of the central axis and the lateral deflection are zero, and are either hinged supported without shear restraints, hard hinged supported, or clamped. Free-end boundary conditions are not considered because in moderately large beam vibrations the axial force N is primarily a result of nonlinear stretching of the central fiber due to longitudinally fully restrained ends. In Eqs. (31) and (32), the subscript b indicates quantities of the boundaries at x = 0 and x = l.
Hinged support without shear restraints When at a hinged end no shear restraints are applied (i.e., the slip between the layers is not restrained), the rotations of the layerwise cross section are not restrained, and consequently, the layerwise bending moments are zero, (M i ) b = 0, i = 1, 2, 3. Eq. (7) facilitates to express this dynamic boundary condition in terms of the kinematic quantity w as At any hinged end, the overall bending moment is zero, M b = 0. That is, according to Eq. (21) and considering Eq. (33), Since the interlayer slip is not constrained at the boundaries, at the beam end the horizontal support reaction, which corresponds to the axial force N , is fully transferred into the central layer, i.e., (N 2 ) b = N . Inserting Eqs. (8) and (18) into this relation yields the fifth boundary condition, which couples all kinematic variables of the problem.

Hard hinged support
At a hard hinged support, an end plate prevents the relative displacement of the layers at the interface, i.e., Thus, the shear tractions at the interfaces (t s12 ) b and (t s23 ) b are also zero. The boundary condition M b = 0 can be expressed alternatively as compared with Eq. (21), Rigidly clamped end A rigidly clamped end yields the slope of the lateral deflection zero, and slip at both interfaces is constrained, It should be noted here that the coupled differential equations Eqs. (24) and (28) can be condensed to a single sixth-order equation of motion in terms of deflection w. In this equation (Eq. (55)), the corresponding boundary conditions and stress resultants in terms of w are derived in Appendix A.

Solution procedure
The coupled set of differential equations of motion (Eqs. (20), (24), (28), (29)) is solved by expanding the lateral deflection w(x, t) into the first N M mode shapes of the corresponding linearized beam problem (i.e., N = 0) (see e.g., [1]), The mode shapes φ n (x) are obtained easily from the homogeneous differential equation of motion of sixth order in terms of lateral deflection w of the corresponding linear beam problem (i.e., Eq. (55) with N = 0 and q = 0), as described comprehensively in [13]. The modal expansion of w is inserted into the ordinary differential equations (20) and (28), which are subsequently solved in combination with Eq. (29) and the corresponding boundary conditions for Δu 12 , Δu 23 , and u (0) 2 as a function of Y n , n = 1, . . . , N M. With these variables available, the axial force N is evaluated according to Eq. (18), yielding a nonlinear series in terms of Y n . For this analysis, any value of 0 ≤ x ≤ l can be employed because N is constant along the span l. The required derivatives of these in such a manner obtained series for w, Δu 12 , Δu 23 , and N are inserted into the partial differential equation Eq. (24). According to the rule of Galerkin (see e.g., [32]), this equation is successively multiplied by φ m (x), m = 1, . . . , N M and integrated over the span l. The orthogonality conditions of the mode shapes φ n (x) simplify considerably the resulting N M nonlinear ordinary differential equations in time t for unknown coordinates Y m , m = 1, . . . , N M, which are coupled through the effect of the axial force N . Eventually, these equations are solved by means of standard procedures of numerical analysis.
Subsequently, this procedure of analysis is exemplarily outlined for a hinged beam without shear restraints. Since the corresponding mode shapes are simply sine waves [3], the solution of the ordinary differential equations Eqs. (20), (28), (29) in combination with the corresponding boundary conditions Eqs.
whereas Δu 12 and Δu 23 contain a linear and a quadratic term of unknown coordinates Y n : The first derivative of the latter equations and of Eq. (40) is inserted into Eq. (18). Evaluation of this expression at, for instance, x = 0 yields after some algebra the axial force N in terms of Y 2 r (r = 1, . . . , N M), Multiplying the mode expanded fourth-order partial differential equation (24) by the mth mode shape, φ m , integrating subsequently over the span l, and considering the orthogonality relations of the mode shapes yields a coupled set of N M nonlinear ordinary differential equations for the modal coordinates Y n , In this equation, ω n denotes the nth natural circular frequency of the corresponding linearized simply supported beam with interlayer slip (see e.g., [3]), with the parameter α according to Eq. (56). P n is the nth modal load, Viscous damping has been added modally to Eq. (47) via the modal damping coefficient ζ n [3].

Illustrative examples
In the following, the nonlinear dynamic response of a three-layer beam whose ends are resting on hinged supports without shear restraints, subjected to a time harmonic half-wave sine load distribution with excitation frequency ν, is analyzed. The load q and the fundamental mode shape φ 1 are affine, and thus, all modal loads except the fundamental one are zero, For internal resonance, higher modes may contribute significantly to the total response due to mode interaction caused by the cubic terms in the modal equations. In the present study, internal resonance is not studied, and since P n = 0 ∀ n = 2, . . . , ∞, the contribution of the higher modes to the nonlinear forced dynamic response is negligible. Consequently, the coupled modal equations of motion reduce to a single equation in terms of the fundamental modal coordinate only: In particular, a three-layer beam of length l = 1.0 m, layer thickness h 1 = h 3 = 0.01 m, h 2 = 0.0102, and width b = 0.1 m is considered. Young's modulus of the external layers, E 1 = E 3 = 7.0 · 10 10 N/m 2 , is seven times larger than of the core, E 2 = 1.0 · 10 10 N/m 2 . To the interfaces, a slip modulus of K = 1.0 · 10 9 N/m 2 is assigned. This setup of the cross section corresponds to a lateral composition parameter α (Eq. (56)) times l of αl = 13.3, indicating a moderate lateral layer interaction [12]. The mass density of the top and the bottom layer ρ 1 = ρ 3 is = 2700 kg/m 3  where v is the velocity at a point with current spatial coordinates x [7]. For the numerical integration of D, the approach presented in [19] is used.
An eigenfrequency analysis yields the first five natural frequencies of the corresponding geometric linear FE model as follows, Comparing these outcomes with the corresponding results of the beam theory, Eq. (53), shows that the first natural frequency of both approaches differs only by about 0.2%. This small difference increases to about 0.8% for the fifth frequency.
The subsequent figures show the undamped nonlinear transient response of the beam. At first, in Fig. 4 the lateral deflection at midpan, w(x = 0.5l), normalized with respect to the nonlinear static response at midspan, w S (x = 0.5l), is shown as a function of the ratio time t over fundamental period T 1 = 2π/ω 1 . The bold full black line represents the nonlinear outcome of the proposed beam theory. As observed, the midspan deflection increases until the maximum is obtained at t/T 1 = 4.763, then decreases to almost zero, then increases again, and so on. This kind of beat effect is a result of the geometric nonlinearity of the problem that yields the fundamental frequency amplitude dependent. In this figure, also the normalized midspan deflection of the  Fig. 5, also verifies for this example problem the accuracy of the proposed beam theory. Figure 6 shows the upper and lower interlayer slip, Δu 12 and Δu 23 , at the left beam end (i.e., x = 0) with respect to t/T 1 . Δu 12 and Δu 23 are also normalized with respect to nonlinear static deflection at midspan, w S (0.5l), to keep the difference in the order of magnitude of the considered displacement response quantities apparent. One interesting observation is that the nonlinearity due to moderately large vibrations causes the upper (black line) and lower (blue line) interlayer slip to become different. In contrast, in the linear beam both interlayer slips (red graph) are identical. The distribution of the normalized interlayer slips of the nonlinear beam over x at t/T 1 = 4.763 reveals that close to the beam ends Δu 12 and Δu 23 become different due to the geometric nonlinear membrane force N . Close to the beam center, where the beam deformation is small, both quantities are the same. The FE solution shown in Figs. 6 and 7 confirms again the accuracy of the proposed theory.
The longitudinal displacement of the central fiber u (0) 2 at x = 0.08l, again normalized with respect to w S (0.5l), is displayed in Fig. 8. The peak value of u (0) 2 appears at about x = 0.08l, compared with the red graph in Fig. 7, which shows the distribution of u  In contrast to the previously discussed kinematic response variables, however, the local FE peak values of u 2 (0.08l) are slightly larger than those of the beam response, as seen in Fig. 8. It should be noted that u 2 is the result of stretching of the central fiber, and thus, for the corresponding linear beam zero.
Subsequently, the resulting internal forces are displayed and discussed. Figure 9 shows the time history of the membrane force N , and the layerwise axial forces N 1 , N 2 , and N 3 , respectively, at midspan. The quantities are normalized to the nonlinear axial force N 3 at x = 0.5l. While membrane force N and central layer force N 2 are always in tension (or zero), the axial forces of the faces, N 1 and N 3 , alternate between tension and pressure. However, in contrast to the linear beam, the peaks in tension are larger than in compression, due  In Fig. 11, the overall bending moment M, and the layerwise bending moments M 1 , M 2 , and M 3 at midspan are depicted as a function of time ratio t/T 1 . The static nonlinear overall bending moment M S at x = 0.5l is used to normalize these dynamic bending moments. In addition, the overall bending moment in the center of span l of the corresponding linear beam is also shown. The moments of the external layers, M 1 and M 3 , which are proportional to the beam curvature w ,x x , are identical because these layers have the same bending stiffness. It is readily observed that the dynamic amplification of the overall bending moment at midspan is of the same magnitude as the one of the midspan deflection, compared with Fig. 4. To complete the insight into the nonlinear response behavior of the considered structural member, in Fig. 12 for a certain time instant the normalized overall and layerwise bending moments are plotted as a function of x.
From the results of this example, it can be concluded that the proposed theory for immovably supported composite beams with interlayer slip predicts very accurately the nonlinear response of those members, validated through a comparative FE analysis.

Example 2: Nonlinear resonance
After having validated the proposed theory for a particular member, subsequently linear and nonlinear frequency response functions of the same but 5% damped composite beam are derived by sweeping the excitation frequency ν in the vicinity of the fundamental frequency ω 1 of the corresponding linear structure. At time t = 0, the member is subjected to the harmonic load according to Eq. (50), and a time history analysis is conducted. After decay of the transient vibrations, the maximum of the steady-state response is recorded. and 10q (re f ) 0 /3, the resonance curves exhibit in a certain frequency range multivalued amplitudes, and the entire solution splits into two stable and one unstable branch. However, in this and the subsequent figures only the stable response branches are depicted because the response has been found by time history analyses. Two different frequency sweeps are conducted to determine the two stable solutions. In the first sweep, starting at ν/ω (re f ) 1 = 0.1, the excitation frequency is stepwise increased, and the last response of the current step is used as initial condition for the subsequent analysis with increased excitation frequency. In the second sweep, starting at ν/ω (re f ) 1 = 2.5 the excitation frequency is stepwise reduced. Both sweeps are continued until that frequency where the well-known jump phenomenon occurs, i.e., the tangent of the amplitude functions becomes vertical. For the load case 2q A similar response behavior is observed for the normalized steady-state longitudinal displacement max |u (0) 2 p (0.08l)|/w SL (0.5l), as seen in Fig. 15. This quantity is zero for the linear member and grows with increasing nonlinearity related to an increase in the load amplitude. Figure 16 represents the resonance functions of the overall axial force N , and the layerwise axial forces N 1 , N 2 and N 1 at midspan only for reference load q (re f ) 0 . All steady-state internal forces are divided by the static axial force in the bottommost layer, at x = l/2 of the corresponding linear beam subjected to q (re f ) 0 , referred to as reference axial force N (re f ) 3SL (0.5l). Thus, in this representation the difference in magnitude of the various axial forces is maintained. As observed, the maximum of the layer quantity N 3 is about 2.69 times larger than of the resultant axial force N , and 38.8 times larger than N 2 . To quantify the increase in the non-dimensional steady-state overall axial force with increasing load, in Fig. 17 the ratio max |N p | over N decreases slightly for the two larger load amplitudes, as revealed by Fig. 18. The amplification of this quantity in the linear beam is 10.0, and thus, smaller than in the three considered nonlinear responses. The resonance curves of the overall bending moment at midspan shown in Fig. 19 are similar to the ones of the midspan deflection, compared with Fig. 13.
Subsequently, the linear and nonlinear frequency response functions of the considered flexibly bonded beam (αl = 13.3) subjected to reference load q The deflection frequency response functions depicted in Fig. 20 represent the dynamic (de-)amplification with respect to the midspan deflection of the flexibly bonded linear beam, w (re f ) SL (0.5l). This figure reveals the large effect of the interlayer stiffness K on the global stiffness, and consequently, on the nonlinear dynamic structural response. As observed, the nonlinear quasistatic midspan deflection of the beam without bonded layers is 6.3 times larger than the static deflection of the flexibly bonded beam, and the nonlinear peak amplitude ratio is 21.0 compared to 9.15 of the reference structure. Because of the large flexibility of this member, its response characteristics becomes more nonlinear compared to the reference beam. That is, in a certain frequency range the response is multivariate with two stable and one unstable branch, and the effect of subharmonics at about one third of the fundamental frequency is more pronounced. By contrast, the fully bonded beam exhibits a nonlinear quasistatic midspan deflection of 0.63 times w (re f ) SL (0.5l), and a peak deflection amplitude amplification of 5.54 is predicted. In the flexibly and the rigidly bonded beam, the nonlinear peak deflection amplification is slightly smaller than in the linear counterpart, however, in the unbonded beam the difference between the linear and nonlinear peak amplification is large.
According to Fig. 21, the nonlinear steady-state interlayer slip amplitude max |Δu 23 p | at x = 0 of the unbonded structure is 6.66 times larger than the one of the flexibly bonded reference beam. It is also seen that in the unbonded nonlinear beam the peak amplification max |Δu 23 p |(0)/w (re f ) SL (0.5l) is much smaller than in the linear member, whereas in the flexibly bonded beam this response behavior is the other way around. Figure 22 shows frequency response functions for the overall axial force N , normalized with respect to the static axial force of the bottom layer at midspan of the linear flexibly bonded beam, N (re f ) 3SL (0.5l). The displayed results show that the maximum nonlinear overall axial force N of all beam configurations is of similar order of magnitude.
Since in the unbonded beam, no shear tractions can be transferred from the central layer to the external layers across the interfaces, the axial forces in the external layers, N 1 and N 3 , are zero. Thus, Fig. 23 shows  22 Frequency response function of the normalized overall axial force. Variation of the interlayer stiffness the midspan frequency response function of N 3 for the flexibly and the rigidly bonded beam. One interesting observation is that for both members the quasistatic internal force is almost the same. The maximum nonlinear response amplification of N 3 is 11.1 for the flexibly bonded member, compared to 12.2 for the beam without interlayer slip. The nonlinear peak of the frequency response function of N 3 (0.5l) exceeds the linear one, which is in contrast to the response behavior of the peak deflection.
In the last figure, Fig. 24, the nonlinear and linear steady-state overall bending moment amplitude at x = l/2 divided by the static bending moment of the linear flexibly bonded beam is shown. While the nonlinear peak bending moments of the flexibly and the rigidly bonded structure (of about 9.4) are about the Nonlinear and linear response same, the maximum amplification of the unbonded beam is more than three times smaller compared to the other members. The reason is the large flexibility of the unbonded beam, by which means the membrane force N = N 2 becomes dominant in the transfer of the external load to the supports. Note that the resonance peak of the overall moment M is identical for the three linear beams because the system is statically determined, and thus, M is unaffected by the structural stiffness.

Summary and conclusions
In this paper, the equations of motion and corresponding boundary conditions for vibrating flexibly bonded three-layer composite beams on immovable supports have been derived. The proposed theory is based on a layerwise application of kinematic Euler-Bernoulli assumptions, a linear elastic relation between the interlayer slip and the interlayer shear traction, and a nonlinear strain-displacement relation for the central fiber. The latter relation accounts for stretching of the central fiber during moderately large lateral beam vibrations, which develop because the supports are rigidly held apart.
In a first illustrative example, it was shown that the proposed beam solution of the considered simply supported member without shear restraints and the outcomes of a comparative plane stress finite element analysis are in excellent agreement. As such validated, in the second illustrative example the proposed beam theory was used to derive frequency response functions for various kinematic response and internal force quantities. The results of this study demonstrate the significant impact of partial layer interaction on the nonlinear dynamic response of composite beams. Consequently, the benefit of the proposed theory for efficiently predicting moderately large vibrations of composite beams with interlayer slip is confirmed.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A.1 Equation of motion
In Eq. (24), Δu 12,x x x + Δu 23,x x x is substituted by the expression obtained by differentiating Eq. (28) with respect to x and solved for Δu 12,x x x + Δu 23,x x x . Then, the resulting equation is differentiated two times with respect to x, and Δu 12,x x x + Δu 23,x x x is replaced by the expression originating from rearranging Eq. (24). This analysis yields the following sixth-order partial differential equation: where the parameter [16] is proportional to the slip modulus K , and thus, defines the degree of lateral composite action. The practical range of this parameter is discussed in [12].

A.2 Stress resultants
Equations (15) and (16)  3,x obtained from rearranging of Eq. (8) yields The difference of the axial forces in the upper and the lower layer follows from rewriting of Eq. (9) and substituting the constitutive equation (7), Inserting Eq. (58) and its second derivative with respect to x into Eq. (57) leads to Eventually, the overall bending moment M results from combination of Eq. (59) and Eq. (14) as The difference of the axial forces in the outer layers N 1 − N 3 , which is subsequently needed for defining the boundary conditions, is found by substituting in Eq. (58) M i by Eq. (7) and M by Eq. (60), Adding Eq. (15) and Eq. (16) reveals that the sum of the shear tractions at the interfaces t s12 + t s23 is equal to −(N 1,x − N 3,x ). Thus, Eq. (61) is differentiated with respect to x and multiplied by (−1), leading to

A.3 Boundary conditions
Hinged support without shear restraints For the sixth-order boundary problem (55), at each beam end three boundary conditions must be prescribed. For a simply supported end, w b = 0, (w ,x x ) b = 0, and M b = 0, compared with Sect. 3.2. Equation (60) is used to express M b = 0 in terms of the variable w, i.e., (E J 0 w ,x x x x − q) b = 0. In summary, the boundary conditions read Hard hinged support The boundary conditions for a hard hinged support are w b = 0, M b = 0, and (Δu 12 ) b = (Δu 23 ) b = 0. The latter conditions can be alternatively expressed as (t s12 + t s23 ) b = 0. Utilizing Eqs. (60) and (62), respectively, the boundary conditions of a hard hinged support in terms of w read as Rigidly clamped end At a rigidly clamped end w b = 0, (w ,x ) b = 0, and (Δu 12 ) b = (Δu 23 ) b = 0, or alternatively (t s12 + t s23 ) b = 0. Expressed in terms of w these conditions become A.4 Discussion Equation (55) in combination with the pertinent boundary conditions can only be solved for w if the axial force N is the result of an external axial force applied at the boundary of the beam, such as in dynamic buckling analysis. In the current problem, however, the overall axial force N develops due to membrane strains as a result of immovable supports. Thus, according to Eq. (18), additionally to w the interlayer slips Δu 12 and Δu 23 , and the axial displacement u (0) 2 need to be solved simultaneously. This is only possible in an efficient manner with the coupled equations present in Sect. 3. However, the formulation of this boundary problem according to Eq. (55) reveals parameter α, which is an important parameter to assess at the outset the vulnerability of the member to interlayer slip [12].