Forced response of double beams with initial geometric imperfections and connected by a nonlinear elastic layer

Abstract

This study investigates the forced vibrations of double simply supported Euler-Bernoulli beams with initial deflections connected through a distributed linear and nonlinear elastic layer. It is assumed that the upper beam is subjected to a concentrated transverse harmonic force located at its midspan. The integro-partial differential equations of motion are introduced. The Galerkin approach is utilized, and the mode shapes of uniform simply supported beams are used to establish the nonlinear governing equations of motion that incorporate quadratic and cubic nonlinearities. The harmonic balance method (HB) in conjunction with the pseudo arc-length continuation scheme are applied to obtain the amplitude-frequency curves. To assess the accuracy and validity of the proposed method and solution, some findings obtained by the suggested approach have been compared with those found through numerical integration, where a strong concordance was demonstrated. The influences of several factors such as the linear and nonlinear stiffness parameters, and the initial deflections of the beams on the steady state amplitudes have been examined. It is observed that these factors have considerable influences on the dynamic responses of the double beams. For the sake of generality and convenience, the results are displayed in dimensionless forms.

Keywords

harmonic force nonlinear vibrations harmonic balance method arc-length continuation frequency response initial deflection

Introduction

The rapid advancements in mechanical, marine, civil, and aeronautical engineering have led to the emergence of new and distinctive structures, including double-beam systems. Double-beam systems are employed to model various applications and structures, including composite and sandwich beams, double-walled carbon nanotubes, vibration absorbers, double-beam cranes, pipelines, and railway tracks situated on elastic or viscoelastic foundations.¹ Double-beam systems typically comprise two parallel beams that are constantly linked by elastic linear or nonlinear springs. These systems serve as effective mechanical analogs for a wide range of real engineering and nanoscale structures. Such configurations are found in sandwich and laminated panels in aerospace, marine, and automotive industries, where two stiff face sheets are bonded to a soft or viscoelastic core that exhibits geometric and material nonlinearities under dynamic loading conditions.^2–4 Moreover, similar formulations are used to describe adhesive-bonded joints and double-cantilever beam tests, where a nonlinear interlayer models adhesive softening, hyperelasticity, or interfacial slip during vibration or debonding.^5,6 At smaller scales, coupled micro- and nano-beam resonators and multi-walled carbon nanotubes can be accurately modeled as double-beam systems with nonlinear coupling due to electrostatic, van der Waals, or interlayer forces.^7,8

In civil and transport engineering, twin-girder bridge decks and railway tracks have been modeled as double beams connected by nonlinear or viscoelastic layers to capture the response of connectors, sleepers, and fasteners under moving or harmonic loads.⁷ Likewise, composite or functionally graded two-layer beams with imperfect interfaces can be represented as double beams with nonlinear interfacial stiffness to account for interlayer slip, hysteresis, and energy dissipation.^4,6 Studies have further shown that initial imperfections or geometric defects—such as curvature, residual stresses, or localized deviations—strongly influence the nonlinear vibration response, often leading to bifurcations and jump phenomena when the upper beam is subjected to harmonic excitation.^2,3,8

It is known that to achieve high-performance and efficient structures, and to enhance the applicability and reliability of proposed designs, it is essential to understand their vibrational behavior thoroughly. Many researchers investigated the dynamic characteristics of double beam systems. For example, Fei et al.⁹ examined the dynamic behavior of a double tilted and tensioned beam system. The two different beams were linked by elastic springs, and the governing equations of motion were formulated by considering the impacts of the flexural rigidity, sag, inclined angle, and boundary conditions. The double beam system was reduced to four degrees of freedom system, and a semi analytical semi numerical method was presented to obtain the frequency equation of the structure.

Furthermore, Murmu and Adhikari¹⁰ explored the stability and vibration analyses of coupled nanobeams under initial compressive pre-stressed conditions and linked by elastic springs. Eringen’s non-local elasticity theory was utilized to find the relations for the transverse vibrations of the double nanobeams. The nanoscale and connecting spring impacts were addressed, and it was observed that they have a major effect on the natural frequencies of the nonlocal double-nanobeams. The findings have shown that the non-local impacts and the rise of the pre-stressed force tend to reduce the natural frequencies. Szmidt¹¹ used experimental and analytical methods to study the dynamics of twin cantilever beams linked by a viscoelastic layer undergoing shear deformation. A continuous system was presented to describe the model, and the stiffness and the damping properties were found.

Li et al.¹² focused on the free and forced oscillations of a double-beam system joint by a viscoelastic member. A semi-analytical approach was used to determine the natural frequencies and the mode shapes. The forced vibration responses were found by using a modal-expansion iterated technique to obtain an orthogonality condition that is needed to decouple the differential equations of motion of the double beam. The reliability of the suggested approach was validated by various examples. Additionally, Nguyen¹³ studied the influence of a point mass location on the natural vibrations of cracked double beams linked by an elastic member. Numerical simulations were carried out, and it was revealed that the frequency of the cracked double-beams intermittently alters when the point mass is placed at the crack location, and the crack position can be identified by the location of peaks of the wavelet transform of the frequency mass location. A normal mode method was applied by Kelly and Srinivas¹⁴ to determine the natural frequencies and the mode shapes for a set of Euler-Bernoulli beams acted upon axial forces and linked by continuous elastic spring layer.

Mirzabeigy et al.¹⁵ used the differential transform method to find natural frequencies and mode shapes of two parallel beams linked by a layer with variable stiffness. It was presumed that the edges of the beams were elastically restrained against translation and rotation by connecting them to translational and rotational springs. The effects of the stiffness layer, flexural rigidity of the beams, and the boundary conditions on the vibration performance were analyzed. Oniszczuk¹⁶ used the Bernoulli-Fourier approach to calculate the natural frequencies of two parallel simply supported beams connected by a Winkler elastic member. Two partial differential equations were used to study the dynamics of the system, and the initial-value problem was considered to obtain the natural frequencies and the corresponding mode shapes. Hadian Jazi¹⁷ used the modified couple stress theory to report the nonlinear forced vibrations of double nanobeams under translating particle load. The system was presumed to be lying on an elastic layer. Hamilton’s principle and Timoshenko beam theory were adopted to derive the partial differential equations of motion. The Galerkin approach was used to obtain the reduced nonlinear temporal equations of motion. A parametric study was performed, and it was illustrated that the material length scale, the stiffness of the elastic layers, and the velocity of the moving particle have significant impacts on the displacements of the double nanobeams.

Rezaiee-Pajand and Hozhabrossadati¹⁸ performed the free vibration analysis of a double-beam system with elastic restraints at one end and free at the other end. The beams were linked by a mass-spring device. The Fourier transform was employed to determine the frequency parameters and the mode shapes of the system under study. Zhou et al.¹⁹ introduced a new Hamiltonian-based method to determine the natural frequencies, the mode shapes, and the displacements of double nanobeams resting on an elastic layer. Besides, Mohanty et al.²⁰ used computational techniques to explore the natural vibrations of a two-layer elastic beam lying on a variable Pasternak layer. Rahman and Lee²¹ used a modified version of the harmonic balance method to study the forced nonlinear vibrations of axially loaded single and double beams. Beams with simply supported and clamped were considered. The Galerkin method was adopted to formulate the ordinary differential equations of motion. The modal amplitudes were expressed as the summation of the zero, first, and second level responses. Hence, a set of decoupled nonlinear algebraic equations is generated at each solution level. The impacts of the damping coefficient, the axial force, and the excitation level on the dynamic behavior of the single and double beam systems were conducted.

Eshtehardiha et al.²² presented a nonlinear energy harvester employing magnetically coupled double cantilever beams with piezoelectric layers. Unlike typical linear harvesters with limited bandwidth, the suggested system exploits magnetic nonlinearity and internal resonance to improve performance. The assumed modes method was used to formulate the vibration equations that were solved numerically. It was revealed that magnetic coupling broadens the harvesting frequency range. Experimental validation confirmed the model’s accuracy. Fang et al.²³ performed the natural and forced oscillations of double Euler-Bernoulli beams. It was assumed that the primary beam carries a tip mass with rotary inertia, linked to a secondary beam through an elastic layer. A parametric study was conducted to examine the impacts of the tip mass, the rotary inertia, the elastic layer stiffness, and the beam rigidity on the system frequencies and dynamic behavior.

Based on the above literature, one can observe that no efforts have been dedicated to conduct the nonlinear forced behavior of double simply supported Euler-Bernoulli beams with initial deflections. Hence, the purpose of this research is to fill this gap. To provide a more structured comparison of the literature, Table 1 has been established to clarify the research gap. In this table, a comparison of 11 key studies based on their beam models, coupling laws, imperfection modeling, and solution methods have been considered. The table explicitly demonstrates that the simultaneous treatment of geometric imperfections and nonlinear interlayers in double-beam systems under forced excitation remains a distinctive contribution of this work.

Table 1.

Comprehensive summary of key studies on the dynamic behavior of double beams.

Reference	Beam model	Coupling law	Imperfection modeling	Excitation type	Solution method	Multi-mode/Internal resonance
Oniszczuk (2000)	Euler-Bernoulli	Linear (Winkler)	None	Free vibration	Bernoulli-Fourier	Yes
Kelly & Srinivas (2009)	Euler-Bernoulli	Linear	None	Free vibration	Normal mode	Yes
Murmu & Adhikari (2012)	Nonlocal Euler-Bernoulli	Linear	None	Free vibration	Analytical	Yes
Li et al. (2016)	Euler-Bernoulli	Linear (Viscoelastic)	None	Free & Forced (concentrated harmonic force acting on the midspan of the upper beam)	Semi-Analytical	Yes
Rahman & Lee (2017)	Euler-Bernoulli	linear	None	Harmonic force acting on the upper beam	Multi-level Residue HB	Yes
Fei et al. (2018)	Tensioned Euler-Bernoulli	Linear	None	Free vibration	Semi-analytical	Yes
Liu & Yang (2019)	Euler-Bernoulli	Linear (Viscoelastic)	None	Free & Forced (harmonic forces)	Closed-form Analytical	Yes
Mirzabeigy & Madoliat (2019)	Euler-Bernoulli	Linear (variable stiffness Winkler-type elastic layer)	None	Free vibration	Semi-Analytical	Yes
Jazi (2020)	Timoshenko	Linear	None	Moving Particle	Galerkin/Numerical	No
Eshtehardiha et al. (2021)	Euler-Bernoulli	Linear	None	Forced (base excitation)	Analytical /Experimental	Yes (Internal Res.)
Fang et al. (2022)	Euler-Bernoulli	Linear	None	Forced	Analytical	Yes
Present Work	Euler-Bernoulli	Nonlinear (Quad/Cubic)	Initial Deflection	Forced (concentrated harmonic force acting on the midspan of the upper beam)	Harmonic Balance / Continuation	No (Fundamental)

Although this study focuses on Euler-Bernoulli beams, the suggested formulation may serve as a basis for several significant extensions. Future advancements may include Timoshenko beam theory to account for shear deformation and rotating inertia, which are essential for the analysis of thicker beams or high-frequency modes. Additionally, the model might be used to investigate the performance of multi-layered systems (more than two beams) or structures that are subjected to moving loads and particles, simulating high-speed transport infrastructure. Besides, considering non-local elasticity effects would allow the formulation to be adapted for nano-electromechanical systems (NEMS), where small-scale impact becomes more significant. As observed in Table 1, most studies on double-beam systems assume perfectly straight members. While some researchers have considered initial pre-stress to analyze stability and frequency shifts, the dynamic performance under the effect of initial geometric imperfections remains under-explored, particularly in the presence of quadratic and cubic interlayer nonlinearities.

The rest of this article is organized as follows. In Section 2, the details of the model are presented, and the governing equations of motion of double simply supported beam with initial deflections are shown. In Section 3, the Galerkin approach is used to obtain the temporal nonlinear equations of motion. Then, the harmonic balance method and the arc length continuation scheme are applied. The nonlinear forced vibrations are investigated. In section 4, the proposed model is validated, and the findings are presented and examined. In section 5, conclusions are drawn, and further developments are suggested in Section 5.

Mathematical modeling

Figure 1 shows double Euler-Bernoulli beams with initial deflections connected through a distributed linear and nonlinear elastic layer. It is assumed that the upper beam is subjected to a harmonic force located at its midspan. In this study, it is assumed that the nonlinear restoring force of the elastic layer is expressed as Huang et al.²⁴

F_{spring} = k_{1} (w_{1} - w_{2}) + k_{2} {(w_{1} - w_{2})}^{2} + k_{3} {(w_{1} - w_{2})}^{3}

(1)

where $k_{1}$ , $k_{2}$ , and $k_{3}$ are the linear stiffness, the quadratic nonlinear stiffness, and the cubic nonlinear stiffness of the elastic layer, respectively.

Figure 1.

Schematic diagram of double Euler–Bernoulli beams with initial deflections connected by a nonlinear elastic layer and subjected to a harmonic force.

The nonlinear coupling layer, defined by the restoring force given in equation (1), is a generalized model designed to represent several physical configurations depending on the engineering application. For example, in aerospace and automotive sandwich panels, the layer models a soft or viscoelastic core or adhesive bonded joints where nonlinearities arise from material behavior such as adhesive softening or hyperelasticity. Moreover, in transport engineering, specifically for railway tracks and twin-girder bridges, the distributed nonlinear layer represents the combined response of connectors, sleepers, and fasteners under dynamic loads. Besides, the formulation is also consistent with double-beam systems situated on nonlinear elastic or viscoelastic foundations, such as pipelines or railway tracks on a Winkler-type base. The linear stiffness is always positive, representing the primary elastic resistance of the coupling layer. An increase in $k_{1}$ generally reduces vibration amplitudes and shifts resonance regions to higher frequencies. Moreover, quadratic nonlinearity is often associated with structural or geometric asymmetry in the system’s deformation or boundary conditions. Physically, it can lead to a “softening” behavior, where the system’s effective stiffness decreases as the oscillation amplitude rises. Crucially, the net quadratic nonlinearity of the system is not determined by $k_{2}$ alone. It depends on the interaction between the coupling stiffness and the initial displacements of the beams. Depending on the orientation and magnitude of these initial curvatures, the net quadratic effect may be positive or negative, shifting the equilibrium position and dictating the direction of the mean component shift in the vibration signal. Additionally, cubic nonlinearity is typically linked to mid-plane stretching and geometric nonlinearities. A positive leads to “hardening” behavior, where the restoring force increases more than linearly with displacement, causing the resonance frequency to rise at higher amplitudes.

With the aid of Zhai et al.,²⁵ the governing partial differential equations of motion for double Euler-Bernoulli beams with initial deflections and subjected to a harmonic force are given as

\begin{matrix} ρ A {\overset{\cdot\cdot}{w}}_{1} + η {\overset{\cdot}{w}}_{1} + E I w_{1}'''' \\ - \frac{E A}{2 L} (w_{1} ″ + w_{1, 0} ″ (x)) \int_{0}^{L} ({' w}_{1}^{2} + 2 w_{1}' w_{1, 0}') dx \\ + k_{1} (w_{1} - w_{2}) + k_{2} {(w_{1} - w_{2})}^{2} + k_{3} {(w_{1} - w_{2})}^{3} \\ = F_{0} \cos (Ω t) δ (x - \frac{L}{2}) \end{matrix}

(2)

\begin{matrix} ρ A {\overset{\cdot\cdot}{w}}_{2} + η {\overset{\cdot}{w}}_{2} + E I w_{2}'''' \\ - \frac{E A}{2 L} (w_{2} ″ + w_{2, 0} ″ (x)) \int_{0}^{L} ({' w}_{2}^{2} + 2 w_{2}' w_{2, 0}') dx \\ + k_{1} (w_{2} - w_{1}) + k_{2} {(w_{2} - w_{1})}^{2} + k_{3} {(w_{2} - w_{1})}^{3} = 0 \end{matrix}

(3)

where $ρ$ is the density of the beams’ material, E is the modulus of elasticity, I is the second moment of area, A is the cross-sectional area of the beams, L is the beams’ length, $w_{i}$ is the transverse displacement of beam i, $w_{i, 0}$ is the initial deflection of beam, $k_{1}$ , $k_{2}$ , $k_{3}$ are the linear, quadratic nonlinear, and cubic nonlinear stiffness parameters of the elastic layer, respectively, $δ ()$ is the Dirac delta function, $F_{0}$ is magnitude of the transverse force, $Ω$ is the excitation frequency, t is time, and x is the axial direction of the beam. It is worthwhile stating that the double beams considered in the present work are identical. The prime indicates the partial derivative with respect to $x$ whereas the dot indicates the derivative with respect to time. For interested readers, the derivation of equations (2) and (3) is shown in Appendix 1.

To simplify the equations of motion, the following dimensionless parameters are introduced as

\begin{matrix} x^{*} = \frac{x}{L}, w^{*} = \frac{w}{L}, t^{*} = t \sqrt{\frac{EI}{ρ A L^{4}}}, Ω^{*} = Ω \sqrt{\frac{ρ A L^{4}}{EI}}, \\ F_{0}^{*} = \frac{F_{0} L^{3}}{EI}, c^{*} = \frac{η L^{2}}{\sqrt{ρ AEI}}, K_{1} = \frac{k_{1} L^{4}}{EI} \\ K_{2} = \frac{k_{2} L^{5}}{EI}, K_{3} = \frac{k_{3} L^{6}}{EI} \end{matrix}

(4)

These dimensionless parameters are substituted into equations (2) and (3) to convert the governing equations of motion into a non-dimensional form. It should be stated that the asterisk notations are removed for simplicity. The governing equations of motion are expressed as

\begin{matrix} {\overset{\cdot\cdot}{w}}_{1} + w_{1}^{''''} + c {\overset{\cdot}{w}}_{1} \\ - \frac{1}{2} (w_{1}^{''} + w_{1, 0}^{''}) \int_{0}^{1} ({' w}_{1}^{2} + 2 w_{1}' w_{1, 0}') dx \\ + K_{1} (w_{1} - w_{2}) + K_{2} {(w_{1} - w_{2})}^{2} + K_{3} {(w_{1} - w_{2})}^{3} \\ = F_{0} \cos (Ω t) δ (x - \frac{1}{2}) \end{matrix}

(5)

\begin{matrix} {\overset{\cdot\cdot}{w}}_{2} + w_{2}'''' + c {\overset{\cdot}{w}}_{2} \\ - \frac{1}{2} (w_{2}'' + w_{2, 0}'') \int_{0}^{1} ({' w}_{2}^{2} + 2 w_{2}' w_{2, 0}') dx \\ + K_{1} (w_{2} - w_{1}) + K_{2} {(w_{2} - w_{1})}^{2} + K_{3} {(w_{2} - w_{1})}^{3} = 0 \end{matrix}

(6)

The authors would like to emphasize that the mathematical background presented here for double-beam systems with initial deflections can be adapted for emerging technologies. For example, the nonlinear coupling model can be used to design nonlinear energy harvesters, where the interaction between parallel beams is tuned to broaden the frequency bandwidth for improved power scavenging. Moreover, the inclusion of the initial deflections provides a fundamental metric for structural health monitoring (SHM). By understanding how these geometric changes shift resonance peaks and trigger bifurcations, engineers can find early signs of damage or delamination in composite sandwich panels used in aerospace and marine hulls. These systems also serve well models for micro-electromechanical systems (MEMS) resonators, where nonlinear interlayer forces like van der Waals or electrostatic interactions control the sensitivity and the stability of these devices.

Solution procedure

The nonlinear partial differential equations (5) and (6) can be discretized in space and time utilizing the Galerkin’s approach as²⁶

w (x, t) = \sum_{j = 1}^{N} ϕ_{j} (x) q_{j} (t)

(7)

where N denotes the number of retained modes, $ϕ_{j} (x)$ are the linear mode shapes for a simply supported beam, and $q_{j} (t)$ are the generalized coordinates. In the current study, it is assumed that the concentrated harmonic force will only excite the fundamental mode.^27,28 Moreover, it is supposed that this is the dominant mode for the double beams system, and higher modes can be ignored. Furthermore, it is presumed that this mode is not involved in an internal resonance with any other higher modes.²⁹ Therefore, other modes not being directly or indirectly excited, will decay to zero with time due to the presence of damping. Thus, the beams displacements are given as

w_{1} (x, t) = ϕ_{1} (x) q_{1} (t)

(8)

w_{2} (x, t) = ϕ_{1} (x) q_{2} (t)

(9)

where the first bending mode shape of the beams is given as

ϕ_{1} (x) = \sin (π x)

(10)

The initial deflections for the simply supported beams are given as³⁰

w_{i, 0} (x) = g_{i} \sin (π x)

(11)

where $g_{i}$ is the value of the initial deflection at the midsection of beam i, where i = 1 for the upper beam, and 2 for the lower beam. Equations (8)–(11), are substituted into equations (5) and (6), then the result is multiplied by $ϕ_{1} (x)$ and integrated over the length of the beam to obtain the nonlinear ordinary differential equations that are given as

\begin{matrix} {\overset{\cdot\cdot}{q}}_{1} + c {\overset{\cdot}{q}}_{1} + (π^{4} + \frac{1}{2} g_{1}^{2} π^{4}) q_{1} \\ + \frac{3}{4} g_{1} π^{4} q_{1}^{2} + \frac{1}{4} π^{4} q_{1}^{3} + K_{1} (q_{1} - q_{2}) \\ + \frac{8}{3 π} K_{2} {(q_{1} - q_{2})}^{2} + \frac{3}{4} K_{3} {(q_{1} - q_{2})}^{3} = 2 F_{0} \cos (Ω t) \end{matrix}

(12)

\begin{matrix} {\overset{\cdot\cdot}{q}}_{2} + c {\overset{\cdot}{q}}_{2} + (π^{4} + \frac{1}{2} g_{2}^{2} π^{4}) q_{2} + \frac{3}{4} g_{2} π^{4} q_{2}^{2} \\ + \frac{1}{4} π^{4} q_{2}^{3} + K_{1} (q_{2} - q_{1}) + \frac{8}{3 π} K_{2} {(q_{2} - q_{1})}^{2} \\ + \frac{3}{4} K_{3} {(q_{2} - q_{1})}^{3} = 0 \end{matrix}

(13)

Harmonic balance

The periodic solutions of the equations of motion can be obtained using the harmonic balance method. The periodic solutions are assumed as

q_{i} (t) = a_{i 0} + \sum_{j = 1}^{N} (a_{ij} \cos (j Ω t) + b_{ij} \sin (j Ω t))

(14)

where $a_{i 0}$ is a constant that represents the solution mean of the generalized coordinate $q_{i} (t)$ , N is the number of harmonics retained in the approximation, $a_{ij}$ are the coefficients of the cosine terms, and $b_{ij}$ are the coefficients of the sine terms. The suggested solution given in equation (14) is substituted into the equations of motion, multiplied by each harmonic, then the result is integrated with respect to the time t from 0 to $2 π / Ω$ , leading to a system of $2 N + 1$ nonlinear algebraic equations, expressed as

f (X, Ω) = 0

(15)

where

X = {\begin{matrix} \begin{matrix} a_{i 0} & a_{i 1} & \dots \end{matrix} & \begin{matrix} a_{ik} & b_{i 1} & \begin{matrix} \dots & b_{ik} \end{matrix} \end{matrix} \end{matrix}}^{T}

(16)

Arc length continuation method

In nonlinear systems, it is essential to obtain periodic steady-state responses and to identify jump phenomena and nonlinear resonances.^31,32 Several numerical methods can be employed, including Incremental Harmonic Balance (IHB), Newton–Raphson continuation on frequency, path-following algorithms, and frequency stepping with stability switching. In this article, the pseudo arc-length continuation scheme is employed to find all possible stable and unstable periodic solutions of X (at a particular value of $Ω$ ).^33,34 Introducing the arc-length parameter $η$ , the augmented equation is written as

f (X, Ω) - η = 0

(17)

Using a quadratic form of the pseudo arc length, equation (17) is rewritten as

f (Y) - η = Y^{' T} (Y - Y_{k - 1}) - η = 0

(18)

where $Y = {[X, Ω]}^{T}$ .

The slope of the equilibrium path ( $Y'$ ) is expressed in terms of two adjacent previous points $Y_{k - 2}$ and $Y_{k - 1}$ as

Y' = \frac{Y_{k - 1} - Y_{k - 2}}{‖ Y_{k - 1} - Y_{k - 2} ‖}

(19)

The initial estimation of the next point is determined as

Y_{u} = Y_{k - 1} + Δ η Y'

(20)

To obtain the periodic solution, the prediction obtained in equation (20) is corrected iteratively utilizing the Newton-Raphson method until the exact solution on the equilibrium path is found. An adaptive pseudo-arc-length continuation method is utilized to trace the solution branches of the nonlinear system. To ensure robust convergence across high-curvature regions and turning points, the algorithm employs a linear predictor that extrapolates the subsequent state from the previous two solution points. The linear predictor is coupled with a Newton-Raphson corrector that is configured with a residual tolerance of 10⁻¹¹ and a maximum limit of 60 iterations per step. The system Jacobian is evaluated numerically through a finite difference scheme. The step size is managed dynamically as the arc-length increment is automatically reduced via bisection upon a convergence failure in the Newton loop and increased by a certain factor following successful iterations to maintain computational efficiency. Turning points and bifurcations are detected by monitoring sign changes in the determinant of the augmented Jacobian matrix. Upon detection of such a transition, the scheme employs a recursive bisection improvement to determine the critical coordinate where the determinant of the Jacobian is approximately zero. This adaptive logic further enables the solver to navigate “jumps” or gaps in the solution manifold by locally enhancing the arc-length resolution until convergence is achieved.

Results

In this section, the frequency response curves of double simply supported beams are presented. The upper beam is presumed to experience a centrally concentrated harmonic force. It is known that a convergence test is essential to ascertain the adequate number of harmonics utilized in the study. Consequently, the steady-state amplitudes are ascertained and graphed against the excitation frequency for varying harmonic counts for a double simply supported beam with parameters $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , $c = 0.1$ , and $F_{0} = 1$ . Figure 2 indicates that three harmonics are enough for convergence, and they can capture the vibrational behavior of double beams. Consequently, the generalized coordinates $q_{1} (t)$ and $q_{2} (t)$ are expressed as

\begin{matrix} q_{1} (t) = a_{10} + a_{11} \cos (Ω t) + b_{11} \sin (Ω t) \\ + a_{12} \cos (2 Ω t) + b_{12} \sin (2 Ω t) \\ + a_{13} \cos (3 Ω t) + b_{13} \sin (3 Ω t) \end{matrix}

(21)

\begin{matrix} q_{2} (t) = a_{20} + a_{21} \cos (Ω t) + b_{21} \sin (Ω t) \\ + a_{22} \cos (2 Ω t) + b_{22} \sin (2 Ω t) \\ + a_{23} \cos (3 Ω t) + b_{23} \sin (3 Ω t) \end{matrix}

(22)

Figure 2.

The frequency response of $q_{1} (t)$ and $q_{2} (t)$ obtained using the harmonic balance method (HB), for different numbers of harmonics, of a double simply supported beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , $c = 0.1$ , and $F_{0} = 1$ .

In the current research, the expressions shown in equations (21) and (22) will be utilized to examine the dynamic behavior of the double beams system at different parameters.

A comparative study has been conducted to validate the accuracy of the proposed solution derived from the harmonic balance method (HB). Figures 3 to 6 compare the frequency amplitude curve derived from HB with those produced using numerical integration (NI) for a double simply supported beam, with the parameters used in Figure 1. The solid line represents the steady state amplitudes derived from the HB, whereas the circles denote the results produced by the NI approach. The efficacy of the proposed method can be corroborated by observing that the outcomes of both approaches align closely. In all presented figures, the force and the damping coefficients are taken as $F_{0} = 1$ , and $c = 0.1 .$ The reduced-order nonlinear model obtained via the Galerkin discretization was solved using the harmonic balance method combined with the pseudo arc-length continuation scheme. To verify the accuracy of this semi-analytical approach, the same governing nonlinear ordinary differential equations were independently solved using direct numerical integration in the time domain. Numerical validation was conducted using the MATLAB ode45 integrator, based on the explicit Runge-Kutta method. A time step of $Δ τ = 0.01$ was utilized, with relative and absolute tolerances maintained at $10^{- 8}$ and $10^{- 10}$ , respectively, to ensure high precision. To isolate the steady-state response, the first 60% of the 600-unit dimensionless time history was discarded to remove all transient effects. The steady-state amplitude was then extracted using the maximum absolute values of the steady states. This extraction method was selected because it accurately captures the peak physical reach of the beams, accounting for the asymmetric oscillations (mean component shift) induced by the initial geometric imperfections.

Figure 3.

The frequency responses of $q_{1} (t)$ and $q_{2} (t)$ obtained using the harmonic balance method (HB) compared with that obtained using direct numerical integration (NI) of a double simply supported beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , $c = 0.1$ , and $F_{0} = 1$ .

Figure 4.

The frequency responses of (a) $q_{1} (t)$ and (b) $q_{2} (t)$ obtained using the harmonic balance method (HB) compared with that obtained using direct numerical integration (NI) of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 10$ , $g_{1} =$ $g_{2} = 2.0$ , $c = 0.1$ , and $F_{0} = 1$ .

Figure 5.

The frequency responses of (a) $q_{1} (t)$ and (b) $q_{2} (t)$ obtained using the harmonic balance method (HB) compared with that obtained using direct numerical integration (NI) of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 2$ , $g_{1} =$ $g_{2} = 0.5$ , $c = 0.1$ , and $F_{0} = 1$ .

Figure 6.

The frequency responses of (a) $q_{1} (t)$ (b) and $q_{2} (t)$ obtained using the harmonic balance method (HB) compared with that obtained using direct numerical integration (NI) of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 10$ , $g_{1} =$ $g_{2} = 0.5$ , $c = 0.1$ , and $F_{0} = 1 .$ .

Then, the steady-state amplitudes derived from numerical integration were subsequently compared with those obtained from harmonic balance solutions. The comparison exhibited excellent agreement over the entire range of excitation frequencies considered, including near resonance peaks and turning points where nonlinear effects such as jump phenomena and multi-valued solutions occur. This agreement confirms that the proposed solution strategy accurately captures both stable and unstable branches of the nonlinear frequency–response curves. By investigating Figures 3 to 6, one can observe that validation has been done for multiple cases spanning different values of $K_{1}, K_{2}, K_{3}$ , and at various levels of initial deflections $g_{1}$ , $g_{2}$ .

In the present study, the authors first calculate a periodic steady-state solution using the Harmonic Balance (HB) method and then use the values of this solution at the initial time ( $t = 0$ ) as the starting conditions for the time-domain numerical simulation. In other words, numerical integration does not start from arbitrary initial conditions (such as zero displacement and velocity), but instead from a state that already lies on the periodic response predicted by the HB method. This approach ensures that the numerical time integration starts exactly on the periodic orbit where no long transient decay is needed, and where the numerical solution directly reproduces the steady-state response. This is especially important in nonlinear systems, where multiple stable and unstable periodic solutions may exist, and where poor initial conditions may lead to convergence to the wrong attractor. From a physical perspective, the structure is assumed to already be vibrating in its steady periodic regime when the numerical simulation begins.

In the examination of nonlinear vibrations in beams, quadratic and cubic nonlinearities frequently arise as a result of geometric effects and material properties. These nonlinearities are essential in determining the system’s dynamic behavior. Quadratic nonlinearities, frequently associated with asymmetry in deformation or boundary conditions, can induce softening behavior, resulting in a decline in the beam’s effective stiffness as the amplitude increases. Cubic nonlinearities are typically linked to mid-plane stretching and geometric nonlinearity. This often induces a hardening response, wherein stiffness grows with amplitude. The coexistence of both types of nonlinearities leads to competition, resulting in varied dynamic behaviors due to their combined influence. The system may exhibit a net hardening effect, a net softening effect, or, in certain instances, almost linear behavior where nonlinear terms cancel each other out. This is contingent upon their respective magnitudes and signs. The competition between quadratic and cubic nonlinearities is crucial for practical applications of beam structures, as it directly influences resonance characteristics, stability, and the amplitude-frequency relationship in nonlinear vibration analysis.

Figure 7 illustrates the impact of several values of the linear stiffness parameter $K_{1} = {1, 10, 20}$ on the amplitude-frequency responses for double slightly curved simply supported beams with $K_{2} = 5, K_{3} = 1,$ and $g_{1} = g_{2} = 0.5$ . The vertical axis denotes the amplitudes $A_{1}$ and $A_{2}$ , indicating the maximum displacements of beam 1 and beam 2, respectively, while the abscissa shows the excitation frequency $Ω$ . As $K_{1}$ grows, the amplitudes $A_{1}$ and $A_{2}$ diminish, and the second resonance region shifts to the right as a result of the elevated natural frequencies. Furthermore, the hardening nonlinearity reduces as $K_{1}$ improves. The hardening nonlinearity signifies a nonlinear response in a system marked by an increase in stiffness as displacement rises, leading to an upward curvature of the response curve at higher frequencies. This signifies that the restoring force escalates more than linearly with displacement, leading to an increase in the system’s natural frequency at larger amplitudes. Moreover, Figure 7 illustrates that the $A_{1}$ response curves exhibit loops for $K_{1}$ values of 1 and 10. The amplitude $A_{1}$ attains a minimal value of 0.023 at Ω ≈ 10.96 for $K_{1}$ = 10, and a lesser value of 0.006 at Ω ≈ 11.4 for $K_{1}$ = 20. Furthermore, a region exhibiting coexisting periodic solutions of the response amplitude $A_{2}$ emerges for $K_{1}$ = 1.

Figure 7.

The Influence of the linear stiffness parameter on the frequency-response curves of a double simply supported beam system with $K_{2} = 5, K_{3} = 1,$ and $g_{1} = g_{2} = 0.5$ .

Figure 8 illustrates the influence of the linear stiffness parameter $K_{1}$ on the steady-state amplitudes of double curved simply supported beams, with $K_{2} = K_{3} = 10,$ and $g_{1} = g_{2} = 1.0$ . The maximum amplitudes of the beams (the peaks) diminish as $K_{1}$ rises. The first resonance regions nearly overlap as $K_{1}$ enhances, indicating that $K_{1}$ exerts a negligible effect on the fundamental natural frequency of the double curved simply supported beams. Conversely, the second resonance regions are displaced to the right as $K_{1}$ propagates, revealing that the second natural frequency increases with the rising of $K_{1}$ . The linear stiffness parameter exerts a softening influence on the steady-state amplitude of the first peak, while it imposes a hardening effect on the steady-state amplitude of the second peak. This asymmetry reveals that higher-order modal interactions are more sensitive to coupling stiffness than the fundamental mode. Moreover, the coexistence of softening behavior at the first peak and hardening behavior at the second peak emphasizes that linear stiffness does not uniformly alter system dynamics but instead redistributes nonlinear effects across resonance regions. This observation has practical implications for design strategies targeting specific frequency bands. Softening nonlinearity denotes a nonlinear behavior in which the system’s stiffness diminishes with increasing displacement, resulting in a downward curvature of the response curve at elevated amplitudes. This implies that the restoring force increases at a rate slower than linear with displacement, resulting in a decrease in the system’s natural frequency at greater amplitudes.

Figure 8.

The effect of the linear stiffness parameter on the frequency-response curves of a double simply supported beam system with $K_{2} = K_{3} = 10,$ and $g_{1} = g_{2} = 1.0$ .

Figure 9 depicts the effect of varying nonlinear quadratic stiffness parameter values $K_{2} = {1, 5, 10}$ on the amplitude-frequency responses of double slightly curved simply supported beams, with $K_{1} = K_{3} = 1,$ and $g_{1} = g_{2} = 0.5$ . It is demonstrated that as $K_{2}$ raises, the amplitudes $A_{1}$ and $A_{2}$ augment. Moreover, the efficacy of the hardening nonlinearity diminishes as $K_{2}$ increases. Furthermore, at elevated excitation frequencies, it is shown that the $A_{2}$ response curves exhibit loops, indicating the presence of multiple periodic solutions. This figure confirms the pronounced effect of quadratic nonlinear stiffness on the system response. Increasing this parameter leads to higher vibration amplitudes and progressively weaker hardening behavior. The quadratic nonlinearity introduces asymmetry into the restoring force, which counteracts the cubic hardening mechanism. At higher excitation frequencies, the formation of closed loops in the response curves reveals strong nonlinear coupling and multiple coexisting periodic solutions. This behavior is indicative of increased susceptibility to jump phenomena and dynamic instability, highlighting the destabilizing role of quadratic nonlinearity in double-beam systems.

Figure 9.

The impact of the nonlinear quadratic stiffness parameter on the frequency-response curves of a double simply supported beam system with $K_{1} = K_{3} = 1,$ and $g_{1} = g_{2} = 0.5$ .

Figure 10 illustrates the influence of the nonlinear quadratic stiffness parameter $K_{2}$ on the steady-state amplitudes of double slightly curved simply supported beams, with $K_{1} = K_{3} = 10,$ and $g_{1} = g_{2} = 0.5$ . The figure indicates that as the nonlinear quadratic stiffness parameter improves, the maximum amplitudes of the second resonance region (the second peak) exhibit a tiny increase, but the first peaks remain unchanged. Furthermore, the primary resonance regions demonstrate (almost) linear characteristics for both beams. This may be ascribed to the cancellation of the quadratic and cubic nonlinearities. One can notice that the increase in $K_{2}$ yields lower hardening effect (the frequency-amplitude curve is less bent to the right) on the second resonance region.

Figure 10.

The influence of the nonlinear quadratic stiffness parameter on the frequency-response curves of a double simply supported beam system with $K_{1} = K_{3} = 10,$ and $g_{1} = g_{2} = 0.5$ .

Figure 11 shows the influence of the nonlinear cubic stiffness parameter $K_{3}$ on the steady-state amplitudes of the double slightly curved simply supported beams, where $K_{1} = K_{2} = g_{1} = g_{2} = 1.0$ . It is noticed that the maximum amplitudes of the beams (the peaks) marginally decline as the magnitude of $K_{3}$ enhances. The first resonance regions for both beams demonstrate linear characteristics. The second resonance regions exhibit nonlinear behavior, transitioning from softening to hardening as the value of $K_{3}$ increases. The transition from softening to hardening behavior in the second resonance region demonstrates the dominant role of cubic nonlinearity at higher amplitudes. Unlike quadratic nonlinearity, the cubic term does not significantly affect the fundamental resonance, emphasizing its amplitude-dependent nature. This selective influence suggests that cubic stiffness can be used as an effective tuning parameter for controlling nonlinear resonance behavior.

Figure 11.

Variations of the frequency-response curves for different nonlinear cubic stiffness parameters of a double simply supported beam system with $K_{1} = K_{2} = 1,$ and $g_{1} = g_{2} = 1.0$ .

In Figure 12, the primary resonance regions for both beams are marginally displaced to the left, indicating softening behavior. The behaviors of the steady-state amplitudes of the second peaks for both beams exhibit linearity for $K_{3} = 5$ . This observation confirms that cubic nonlinearity alone does not guarantee hardening behavior; rather, its interaction with other system parameters governs the overall response. The near-linear behavior at higher stiffness values further suggests that nonlinear effects can be deliberately neutralized through parameter balancing.

Figure 12.

Variations of the frequency-response curves for different nonlinear cubic stiffness parameters of a double simply supported beam system with $K_{1} = K_{2} = 10,$ and $g_{1} = g_{2} = 1.0$ .

Figure 13 demonstrates the amplitude-response curves for amplitudes $A_{1}$ and $A_{2}$ , calculated for double curved simply supported beams with $K_{1} = K_{2} = K_{3} = g_{2} = 1$ , across various values of the initial rise of the upper beam, $g_{1}$ . The top beam exhibits hardening behavior at $g_{1}$ =0 and softening behavior at $g_{1}$ =2. Nonetheless, the beam has linear behavior when $g_{1}$ =1. This non-monotonic behavior demonstrates that initial curvature acts as an effective nonlinear stiffness parameter. The rightward shift of resonance regions reflects increased natural frequencies due to geometric stiffening. Furthermore, the resonance regions are displaced to the right as $g_{1}$ enhances, attributable to the elevation of the fundamental natural frequencies. The amplitude $A_{1}$ remains unaffected by the increase in $g_{1}$ . The figure indicates that the amplitude $A_{2}$ grows as $g_{1}$ rises from 0 to 1, then decreasing as $g_{1}$ ascends from 1 to 2.

Figure 13.

The effect of the initial rise of the upper beam on the frequency-response curves of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 1,$ and $g_{2} = 1.0$ .

Figure 14 further shows that the vibration amplitude of the upper beam initially increases and then decreases with rising imperfection magnitude. This trend suggests the existence of an optimal imperfection level that maximizes vibration response. Such sensitivity underscores the necessity of accounting for manufacturing tolerances and geometric defects in practical designs, as small geometric imperfections can significantly alter dynamic performance.

Figure 14.

The influence of the initial rise of the upper beam on the frequency-response curves of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 5,$ and $g_{2} = 1.0$ .

Figure 15 elucidates the impact of the initial rise of the lower beam $g_{2}$ on the frequency-response diagrams of double slightly curved simply supported beams, where $K_{1} = K_{2} = K_{3} = 1$ , and $g_{1} = 0.5$ . The frequency response curves have a rightward bend, indicating hardening behavior. The peak of the $A_{1}$ curve is marginally influenced by the increase of $g_{2}$ . Examining the first peak of the $A_{2}$ curve, the behavior of the frequency-amplitude curve is linear at $g_{2} = 0$ , while it demonstrates a rightward curvature (hardening behavior) at $g_{2} = 1$ and $g_{2} = 2$ , with this nonlinearity becoming increasingly pronounced as $g_{2}$ enhances. The second peak of the $A_{2}$ curve exhibits hardening nonlinear behavior at $g_{2} = 0$ , while the curve appears linear at $g_{2} = 1$ and $g_{2} = 2$ . The second apex of the $A_{2}$ curve diminishes as $g_{2}$ propagates. Figure 15 demonstrates that increasing the initial rise of the lower beam consistently enhances hardening behavior. The progressive bending of the frequency–response curves indicates that geometric imperfections in the non-excited beam still exert strong influence through nonlinear coupling. The reduction in resonance amplitudes at higher imperfection levels suggests enhanced stiffness and reduced energy transfer efficiency.

Figure 15.

The effect of the initial rise of the lower beam on the frequency-response curves for (a) q₁(t) and (b) q₂(t) of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 1,$ and $g_{1} = 0.5$ .

In Figure 16, at elevated values of the linear and nonlinear stiffness parameters ( $K_{1} = K_{2} = K_{3} = 5$ ), the maximum displacement of the upper beam increases, and the hardening nonlinearity intensifies as $g_{2}$ escalates. Additionally, it is demonstrated that the $A_{1}$ response curve creates a loop when $g_{2} = 0$ . Conversely, the maximum displacements of the lower beam (the first and second peaks) diminish as $g_{2}$ rises. The figure shows that at higher stiffness levels, the influence of lower beam imperfection becomes more pronounced, leading to loop formation and multi-valued responses. This behavior confirms that geometric imperfections and stiffness parameters act synergistically rather than independently. The contrasting amplitude trends of the upper and lower beams further highlight complex energy redistribution mechanisms within the coupled system.

Figure 16.

The effect of the initial rise of the lower beam on the frequency-response curves of a double simply supported beam system with $K_{1} = K_{2} = K_{3} = 5,$ and $g_{1} = 0.5$ .

Figure 17 presents the time history, the fast Fourier transform, and the phase portrait of the double-beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , and a pre-jump frequency of $Ω = 11.82$ (the jump occurs at $Ω = 11.838$ ). The steady-state time histories for both the upper beam $q_{1} (t)$ and the lower beam $q_{2} (t)$ reveal a stable periodic oscillation. However, a distinct vertical shift is observed; the beams do not oscillate symmetrically about the zero-equilibrium position. This is a direct physical manifestation of the quadratic nonlinearities introduced by the initial displacements $g_{1}$ and $g_{2}$ . The stability of the envelope over the time span confirms that the system has successfully shed all transient noise and reached a true limit cycle. In the phase plane, the trajectory forms a well-defined, closed limit cycle. The “egg-like” or distorted shape of this cycle is a direct visual signature of the quadratic and cubic nonlinearities; specifically, the quadratic terms induce asymmetry in the orbit, while the cubic terms distort the curvature relative to a linear system’s ellipse. Due to the initial displacements of the beams, the center of this limit cycle is not at the origin. This spatial offset in the phase plane is the visual counterpart to the mean component identified in the frequency domain. The FFT analysis confirms that the response is dominated by the operating frequency of $Ω = 11.82$ , while the significant “mean component” at $Ω = 0$ represents the non-zero average position of the beams during oscillation. The rapid attenuation of the higher harmonics reinforces the validity of the analytical model, proving that the three-harmonic Harmonic Balance expansion is both mathematically sufficient and physically accurate for characterizing the system’s dynamics.

Figure 17.

(a) and (d) The steady state forced responses, (b) and (e) the FFT, (c) and (f) the phase portraits obtained for a double simply supported beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , and a pre-jump frequency of $Ω = 11.82$ .

Figure 18 illustrates the time history, the fast Fourier transform, and the phase portrait of the double-beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , and a post-jump frequency of $Ω = 11.804$ . By initializing the numerical solver on the predicted unstable branch, the simulation successfully captures the physical departure from this saddle-point equilibrium and the subsequent convergence to the stable lower-amplitude attractor, effectively illustrating the “jump” to the lower-energy state. The resulting steady-state response is characterized by a low-amplitude, asymmetric oscillation about a non-zero equilibrium, represented in the frequency spectrum by a significant mean component alongside a dominant fundamental harmonic. As in Figure 17, the periodic steady-state response forms a closed “egg-like” shape which is a direct visual signature of the quadratic and cubic nonlinearities. Because of the initial displacements, the center of the limit cycle is not at the origin. This spatial offset in the phase plane is the visual counterpart to the mean component in the frequency domain. The FFT confirms that the response is dominated by the operating frequency $Ω$ , and the “mean component” or “static offset” at $Ω = 0$ represents the non-zero average position of the beam during oscillation.

Figure 18.

(a) and (d) The steady state forced responses, (b) and (e) the FFT, and (c) and (f) the phase portraits obtained for a double simply supported beam system with $K_{1} = 1, K_{2} = 5, K_{3} = 1$ , $g_{1} = 1$ , $g_{2} = 0.5$ , and a post-jump frequency of $Ω = 11.804$ .

Conclusions

The forced nonlinear response of double simply supported Euler-Bernoulli beam system with initial deflections connected with an elastic nonlinear elastic layer was conducted. In this work, a central harmonic transverse force was acting on the upper beam. The governing integro-partial differential equations of motion were presented, and the system was effectively converted into nonlinear ordinary differential equations through the application of the Galerkin method. The HB approach, along with the pseudo-arclength continuation scheme, was employed to produce frequency-amplitude curves. Certain results were corroborated by comparing them with those obtained from direct numerical integration of the ordinary differential equations of motion.

The findings indicate that increasing the linear stiffness of the connecting layer reduces steady-state amplitudes and shifts resonance regions to higher frequencies while simultaneously reducing hardening nonlinear behavior. In contrast, advancing the quadratic nonlinear stiffness leads to a growth in steady-state vibration amplitudes and a potential diminution of hardening nonlinearity. The influence of cubic nonlinear stiffness was found to be more intricate and dependent on other system parameters; however, its progression clearly facilitates a shift from softening to hardening behavior within the second resonance region. Furthermore, the advancement of both linear and quadratic nonlinear stiffness parameters promotes the formation of multi-valued solutions, including jumps and hysteresis loops, in specific frequency ranges. Ultimately, the initial deflections of the upper and lower beams play a significant role in influencing the system’s linear and nonlinear characteristics, natural frequencies, and steady-state amplitudes.

This work offers a comprehensive analytical framework and significant insights into the nonlinear dynamics of double-beam systems with initial imperfections, emphasizing the crucial effect of linear and nonlinear coupling stiffness, along with geometric initial deflections, on the overall vibrational characteristics. Additionally, the utilization of double beam systems may assist sustainable advancement in industrial applications. This study advances new design and methodology that can improve the resilience and efficiency of infrastructure. It is recommended to analyze the dynamic behavior of other double structures, including plates and rods.

Footnotes

Appendix 1

After neglecting the axial displacement, the longitudinal strain of the mid-plane, accounting for stretching and initial geometric imperfection, is given by³⁵

(23)

ϵ_{0} = w_{0}' w + \frac{1}{2} w^{' 2}

For a double beam system, we define these strains for each beam (i = 1, 2) as

(24)

ϵ_{0 i} = w_{0 i}' w_{i} + \frac{1}{2} w_{i}^{' 2}

The total energy of the system is the sum of the energies of the two beams, and is given as

(25)

T = \sum_{i = 1}^{2} \frac{1}{2} \int_{0}^{L} ρ A {\overset{\cdot}{w}}_{i}^{2} dx

The total bending potential energy $V_{b}$ is expressed as

(26)

V_{b} = \sum_{i = 1}^{2} \frac{1}{2} \int_{0}^{L} EI w_{i}^{' 2} dx

The stretching potential energy $V_{s}$ is given as

(27)

V_{s} = \frac{EA}{2 L} \int_{0}^{L} {(w_{0 i}' w_{i} + \frac{1}{2} w_{i}^{' 2})}^{2} dx

The potential energy of the continuous elastic $V_{k}$ layer is given as

(28)

V_{k} = \int_{0}^{L} (k_{1} {(w_{1} - w_{2})}^{2} + \frac{1}{3} k_{2} {(w_{1} - w_{2})}^{3} + \frac{1}{4} k_{3} {(w_{1} - w_{2})}^{4}) dx

The total potential energy is expressed as

(29)

V = V_{b} + V_{s} + V_{k}

Extended Hamilton’s principle can be written as

(30)

\int_{t_{1}}^{t_{2}} (δ T - δ V) dt + \int_{t_{1}}^{t_{2}} \int_{0}^{L} (F_{i} + η {\overset{\cdot}{w}}_{i}) δ w_{i} dt = 0

Where

(31)

F_{1} = F_{0} \cos (Ω t) δ (x - \frac{L}{2})

(32)

F_{2} = 0

Substituting equations (25)–(28) into equation (30) and performing the integrations yields the following two nonlinear coupled integro-partial differential equations:

(33)

\begin{matrix} ρ A {\overset{\cdot\cdot}{w}}_{1} + η {\overset{\cdot}{w}}_{1} + E I w_{1}'''' \\ - \frac{E A}{2 L} (w_{1} ″ + w_{1, 0} ″ (x)) \int_{0}^{L} (w_{1}^{2} + 2 w_{1}' w_{1, 0}') dx \\ + k_{1} (w_{1} - w_{2}) + k_{2} {(w_{1} - w_{2})}^{2} + k_{3} {(w_{1} - w_{2})}^{3} \\ = F_{0} \cos (Ω t) δ (x - \frac{L}{2}) \end{matrix}

(34)

\begin{matrix} ρ A {\overset{\cdot\cdot}{w}}_{2} + η {\overset{\cdot}{w}}_{2} + E I w_{2}'''' \\ - \frac{E A}{2 L} (w_{2} ″ + w_{2, 0} ″ (x)) \int_{0}^{L} (w_{2}^{2} + 2 w_{2}' w_{2, 0}') dx \\ + k_{1} (w_{2} - w_{1}) + k_{2} {(w_{2} - w_{1})}^{2} + k_{3} {(w_{2} - w_{1})}^{3} = 0 \end{matrix}

Handling Editor: Divyam Semwal

ORCID iD

Ma’en S. Sari

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Liu

Yang

. A closed-form analytical solution method for vibration analysis of elastically connected double-beam systems. Compos Struct 2019; 212: 598–608.

Lee

H-H

. Non-linear vibration of a multilayer sandwich beam with viscoelastic layers. J Sound Vib 1998; 216(4): 601–621.

Kovac

Anderson

Scott

. Forced non-linear vibrations of a damped sandwich beam. J Sound Vib 1971; 17(1): 25–39.

Adali

. Variational principles for a double Rayleigh beam system undergoing vibrations and connected by a nonlinear Winkler–Pasternak elastic layer. Nonlinear Eng 2023; 12(1): 20220259.

Cabello

Zurbitu

Renart

, et al. A non-linear hyperelastic foundation beam theory model for double cantilever beam tests with thick flexible adhesive. Int J Solids Struct 2016; 80: 19–27.

Foraboschi

. Analytical solution of two-layer beam taking into account nonlinear interlayer slip. J Eng Mech 2009; 135(10): 1129–1146.

Mao

. Free vibration analysis of elastically connected multiple-beams by using the Adomian modified decomposition method. J Sound Vib 2012; 331(11): 2532–2542.

Mirzabeigy

Madoliat

. A note on free vibration of a double-beam system with nonlinear elastic inner layer. J Appl Comput Mech 2019; 5(1): 174–180.

Fei

Danhui

Cheng

, et al. Analysis on the dynamic characteristic of a tensioned double-beam system with a semi theoretical semi numerical method. Compos Struct 2018; 185: 584–599.

10.

Murmu

Adhikari

. Nonlocal elasticity based vibration of initially pre-stressed coupled nanobeam systems. Eur J Mech 2012; 34: 52–62.

11.

Szmidt

. Shear deformation damping of a double-beam structure. J Sound Vib 2016; 370: 163–175.

12.

Sun

. Dynamical behavior of a double-beam system interconnected by a viscoelastic layer. Int J Mech Sci 2016; 105: 291–303.

13.

Nguyen

. Crack detection of a double-beam carrying a concentrated mass. Mech Res Commun 2016; 75: 20–28.

14.

Kelly

Srinivas

. Free vibrations of elastically connected stretched beams. J Sound Vib 2009; 326(3-5): 883–893.

15.

Mirzabeigy

Madoliat

Vahabi

. Free vibration analysis of two parallel beams connected together through variable stiffness elastic layer with elastically restrained ends. Adv Struct Eng 2017; 20(3): 275–287.

16.

Oniszczuk

. Free transverse vibrations of elastically connected simply supported double-beam complex system. J Sound Vib 2000; 232(2): 387–403.

17.

Hadian Jazi

. Nonlinear vibration of an elastically connected double Timoshenko nanobeam system carrying a moving particle based on modified couple stress theory. Arch Appl Mech 2020; 90(12): 2739–2754.

18.

Rezaiee-Pajand

Hozhabrossadati

. Free vibration analysis of a double-beam system joined by a mass-spring device. J Vib Control 2016; 22(13): 3004–3017.

19.

Zhou

Fan

, et al. Exact vibration analysis of a double-nanobeam-systems embedded in an elastic medium by a Hamiltonian-based method. Physica E Low Dimens Syst Nanostruct 2018; 99: 220–235.

20.

Mohanty

Bhattacharya

Yadav

, et al. Free vibration analysis of a two-layer elastic beam resting on a variable Pasternak foundation. Mater Today Proc 2018; 5(9): 19242–19249.

21.

Rahman

Lee

. New modified multi-level residue harmonic balance method for solving nonlinearly vibrating double-beam problem. J Sound Vib 2017; 406: 295–327.

22.

Eshtehardiha

Tikani

Ziaei-Rad

. Experimental and numerical investigation of energy harvesting from double cantilever beams with internal resonance. J Sound Vib 2021; 500: 116022.

23.

Fang

Hao

. Free and forced vibrations of an undamped double-beam system carrying a tip mass with rotary inertia. J Eng Mech 2022; 148(2): 04021141.

24.

Huang

Xiao

, et al. A subharmonic modulation-based approach for debonding damage detection of composite stiffened plates. Nondestructive Test Eval 2024; 39(4): 802–831.

25.

Zhai

Ding

, et al. Nonlinear transverse vibrations of a slightly curved beam with hinged–hinged boundaries subject to axial loads. Arch Appl Mech 2022; 92(7): 2081–2094.

26.

Alimoradzadeh

Tornabene

Dimitri

. Nonlinear axial-lateral coupled vibration of functionally graded-fiber reinforced composite laminated (FG-FRCL) beams subjected to aero-thermal loads. Int J Non Linear Mech 2024; 159: 104612.

27.

Tian

Yang

Zhang

, et al. Vibration reduction in beam bridge under moving loads using nonlinear smooth and discontinuous oscillator. Adv Mech Eng 2016; 8(6): 1687814016652564.

28.

Alfahmi

Erturk

. Programmable hardening and softening cubic inductive shunts for piezoelectric structures: Harmonic balance analysis and experiments. J Sound Vib 2024; 571: 118029.

29.

Nayfeh

. Linear and nonlinear structural mechanics. John Wiley & Sons, 2024.

30.

Zhu

Chen

X-C

Guo

, et al. Static and dynamic characteristics of the post-buckling of fluid-conveying porous functionally graded pipes with geometric imperfections. Int J Mech Sci 2021; 189: 105947.

31.

Altan

Hacıoğlu

. Model predictive control of three-axis gimbal system mounted on UAV for real-time target tracking under external disturbances. Mech Syst Signal Process 2020; 138: 106548.

32.

Yanarateş

Okur

Altan

. Performance analysis of digitally controlled nonlinear systems considering time delay issues. Heliyon 2023; 9(10): e20994.

33.

Karličić

Cajić

Paunović

, et al. Periodic response of a nonlinear axially moving beam with a nonlinear energy sink and piezoelectric attachment. Int J Mech Sci 2021; 195: 106230.

34.

Zhang

. Nonlinear vibrations of porous-hyperelastic cylindrical shell under harmonic force using harmonic balance and pseudo-arc length continuation methods. Thin-Walled Struct 2024; 198: 111767.

35.

Lacarbonara

. A theoretical and experimental investigation of nonlinear vibrations of buckled beams. Master’s thesis, Virginia Polytechnic Institute and State University, 1997.