1 Introduction

Buckling of an Euler–Bernoulli beam (E–B) resting on an elastic Winkler foundation has been thoroughly studied from various points of view in many engineering fields for more than 80 years. A comprehensive review on different theoretical elastic and viscoelastic foundation models in oscillatory systems can be found in [1]. The review covers the simplest foundation models to the most complicated ones, and fully describes the recent theories on the topic of mechanical foundations. Special attention in [1] is paid to publications which consider the dynamics of an E–B beam resting on a nonlinear elastic foundation. The dynamics and buckling loads for an E–B beam resting on an inhomogeneous elastic Winkler foundation also plays an importing role in the study of problems of soil–solid interaction [2]. Also E–B beams buckling is of interest for models of energy harvesting E–B beams, in which instability is an important issue [3], and for controllable artificial devices designing [4]. The buckling of beams on an elastic foundation has been discussed in the seminal work of [5]. Other references can be found in books such as [6, 7], and in papers such as [8,9,10]. In [11] two methods for solving the eigenvalue problems of vibrations and stability of a beam on a variable Winkler elastic foundation are presented and compared. The first is based on using the exact stiffness, consistent mass, and geometric stiffness matrices for a beam on a variable Winkler elastic foundation. The second method is based on adding an element foundation stiffness matrix to the regular beam stiffness matrix, for vibrations and stability analysis. In [12] free vibrations of an Euler–Bernoulli beam resting on a variable Winkler foundation is considered. Constant, linear and parabolic variations are considered. The problem is handled for three different boundary conditions: simply supported-simply supported, clamped–clamped and cantilever (clamped-free) beams. The governing differential equations of the beam are solved by using differential transform method. In [13] the buckling and free vibrations of Timoshenko beams resting on variable elastic foundation are analyzed by means of a new finite element formulation. In [14] vibration characteristics of axially functionally graded nanobeams resting on variable elastic foundations are investigated based on a nonlocal strain gradient theory. The authors of [14] considered linear, parabolic, and sinusoidal variations of the Winkler foundation in longitudinal direction. The governing equations in [14] were solved by applying a Galerkin-based solution for different boundary edges. The eigenvalue problem for the buckling loads and natural frequencies of a braced beam on an elastic foundation were studied in [15, 16]. It was found that the location of the translational springs, which were attached to the beam, has significant impact on the buckling loads, on the buckling shapes, and on the eigenfrequencies of the structure. In particular it was shown that under special conditions, an ideal spring stiffness exists, such that the elastic supports do not deflect when the beam buckles. Also, in [16] the eigenvalue problems for the buckling loads and natural frequencies of a braced beam on an elastic foundation were investigated. The conclusion made in [16] is that the study of the eigenvalues variation patterns can offer a design guidance for using a lateral brace of translational springs to strengthen the structure. In [17] a tire model with a flexible belt on an elastic multi-stiffness foundation is investigated via theoretical modeling and an experimental modal analysis. In [18] the vibrations of an isotropic beam on a variable Winkler foundation were investigated by using a modified differential quadrature method. In [19] the natural frequencies and the buckling stresses of a deep beam-column resting on elastic foundations were obtained by using the method of power series expansions of the displacement components and by using a numerical method. In [20] the general solution for a problem describing the vibrations of a beam on a variable Winkler elastic foundation is presented. The exact solution of the dynamic response of the beam is obtained by considering the reaction force of the foundation on the beam as an external force acting on the beam, which is an integral equation including the displacement of the beam. This integral equation was solved approximately and numerically. In [21] the finite difference and the finite element methods are applied to determine the natural frequencies of non-prismatic and non-homogeneous beams, subjected different boundary conditions and resting on a variable Winkler foundation. In this paper, we consider the problem on how to determine the buckling load of an Euler–Bernoulli beam resting on an inhomogeneous elastic Winkler foundation taking into account different damping models (including space hysteresis type of models). Damping effects were considered in many studies (see, for example, [22,23,24,25]). In contrast to previous papers, we investigate a beam on an inhomogeneous elastic foundation. Localization phenomena in one-dimensional imperfect continuous structures were analyzed, both for dynamic cases and for buckling cases in [26,27,28,29]. By using an asymptotic approach we find analytical formulas for the buckling load for some particular types of inhomogeneities in the elastic foundation of the beam. The buckling load depends on the inhomogeneities in the elastic foundation, and this effect exhibits a “space resonance,” i.e., the magnitude of the critical force depends on the inclusion location. This result generalizes the results previously obtained in [15, 16]. Also analytically it will be shown that it is possible that a so-called added mass instability can occur under the influence of damping. The analytical expressions for these “added masses” for different damping models (including the space hysteresis one) will be given. These effects arise as a result of an interaction between the main mode, which is close to instability, and all the other stable modes. This interaction is induced by damping, and we will discuss how it depends on the type of damping model.

2 Statement of the problem

The equation describing the beam dynamics is given by:

$$\begin{aligned} EI u_{xxxx} + a u_{xx} + b(x) u + \epsilon D [u] + m u_{tt} =0 \ , \end{aligned}$$
(1)

where u(x, t) is the beam displacement, \(t \ge 0\) is the time, \(x \in [0, L]\) is the space coordinate in the axial direction of the beam, \(m=A\rho \) is the mass of the beam per unit length, A is the beam’s cross sectional area, \(\rho \) is the beam material density, E is the Young’s modulus, I is the moment of the cross-section inertia, \(\epsilon >0\) is a small parameter, a is a longitudinal force, which can have a positive (compression), or a negative (expansion) sign, and EI is the bending rigidity of the beam. In (1) D denotes a linear operator acting on u(x, t) and its derivatives, and defines damping. The function \(b(x)>0\) is an elastic foundation coefficient. Notice that the differential equation (1), and that given boundary and initial conditions can be transformed to a dimensionless form when we rescale the variables. For the rescaling, the following choice is made: \(x={\bar{x}}L\), \(u={\bar{u}}L\), \(e_0=\frac{I}{AL^2}\), \(t={\bar{t}}\frac{L}{c_0}\), \(c_0^2=\frac{E}{A\rho }\), \({\bar{a}}=\frac{a}{AE}\), \({\bar{b}}=\frac{L^2{b}}{Ac_0^2\rho }\). For simplification, the bar is omitted and the final equation then takes the form:

$$\begin{aligned} e_0 u_{xxxx} + a u_{xx} + b(x) u + \epsilon D [u] + u_{tt} =0 \ , \end{aligned}$$
(2)

where \(x \in [0,1]\). The following initial conditions are considered:

$$\begin{aligned} u(x, 0)=v_0(x), \quad u_{t}(x, 0)=v_1(x), \quad x \in [0, 1], \end{aligned}$$
(3)

where \( || {v_0}_{xx}|| + ||v_1|| < \infty . \) Here, we use the standard notation ||v|| for the norm, \(|| f||^2=\langle f, f \rangle \) , and \(\langle f, g \rangle \) is the scalar product in \(L_2[0, 1]\): \(\langle f, g \rangle =\int _{0}^1 f(x) g(x) \mathrm{d}x\).

The boundary conditions are assumed to be simply supported ones:

$$\begin{aligned} u(0, t)=u(1, t)=0, \quad u_{xx}(0, t)=u_{xx}(1, t)=0. \end{aligned}$$
(4)

However, the applied perturbation approach, which is used in this paper, is applicable to problems with other boundary conditions.

3 Eigenfunctions

Our first step is to consider the unperturbed equation (2) with \(\epsilon =0\). The unperturbed equation (2) can be solved by using Fourier’s method, i.e., by substitution of \(u(x, t)=Re \ \psi (x) \exp (i\omega t)\) into (2) with \(\epsilon =0\). Then, for \(\psi \) we obtain the following eigenfunction problem

$$\begin{aligned} {{\mathcal {L}}} \psi _n =\lambda _n\psi _n, \end{aligned}$$
(5)

where the orthonormal eigenfunctions \(\psi _n\) satisfy the boundary conditions (4).

For \(b(x)=const\) the eigenfunctions are well known (see [5]). To handle the case \(b\ne const\) we apply perturbation theory [30, 31] by assuming that \(b(x)=b_0 + \delta _b b_1(x)\), where \(\delta _b>0\) is a dimensionless small parameter. For positive integers n and \(\lambda _n >0\) let us denote by \( \omega _n\) the natural frequencies, which are defined by \(\omega _n^2 =\lambda _n\), where \(n=1, 2, \ldots \). We set formally \(\omega _{-k}=-\omega _k\). For \(a \le 0\) and the given boundary conditions (4) it is well-known that the eigenvalues \(\lambda _n\) of the operator \({{\mathcal {L}}}\) are always positive. In this case with \(\epsilon =0\) we have stable solutions, where the natural frequencies \(\omega _n\) are positive and \( \omega _n^2 =\lambda _n\). For the opposite case \(a > 0\) an instability is possible, when \(\lambda _n <0\) for a certain n. There exists a critical value \(a_c\) for the parameter a such that for \(a> a_c\) at least one eigenvalue \(\lambda _n <0\) occurs for some n. This can be shown, for example, by the Raleigh variation principle for eigenvalues.

3.1 Unperturbed eigenfunctions

Let \(\delta _b=0\). We denote the corresponding unperturbed eigenfunctions and eigenvalues by \(\psi _n^{(0)}\) and \(\lambda _n^{(0)}\), respectively. In the case of hinged ends we have

$$\begin{aligned} \lambda _n^{(0)}=e_0(\pi n )^4 - a (\pi n)^2 +b_0, \quad n=1, 2, \ldots \end{aligned}$$
(6)

and \(\psi _n^{(0)}=\sqrt{2} \sin (\pi n x)\).

3.2 Perturbed eigenfunctions: non-resonance case

For small positive \(\delta _b\) we can apply the well-known perturbation theory [30, 31]. Under the assumptions

$$\begin{aligned} \varDelta _{kn}=|\lambda _n^{(0)} -\lambda _k^{(0)}| \gg \delta _b, \quad \forall k \ne n \end{aligned}$$
(7)

we have the asymptotical formulas for the eigenvalues

$$\begin{aligned} \lambda _n=\lambda _n^{(0)} + \delta _b \langle b_1 \psi _n^{(0)}, \psi _n^{(0)} \rangle +O(\delta _b^2) , \end{aligned}$$
(8)

and for the perturbed eigenfunctions:

$$\begin{aligned} \psi _n=\psi _n^{(0)} + \delta _b \psi _n^{(1)} + O(\delta _b^2), \end{aligned}$$
(9)

where

$$\begin{aligned} \psi _n^{(1)} = - \sum _{k \ne n} \frac{\langle b_1 \psi _n^{(0)} , \psi _k^{(0)} \rangle \psi _k^{(0)}}{\lambda _k^{(0)} - \lambda _n^{(0)}} . \end{aligned}$$
(10)

Relation (8) can be used to find an approximation for the buckling load \(a_c\) (see below). Note that in the case \(\varDelta _{kn} =O(\delta _b)\) for certain \(k \ne n\) the relations (8), (9), (10) have to be modified. We consider this case in the next subsections.

3.3 On condition (7)

A sufficient condition to satisfy (7), can be found as follows. Using the relation (6), we obtain

$$\begin{aligned} \varDelta _{kn}=e_0\big ((\pi n)^4 -(\pi k)^4\big ) - a \big ((\pi n)^2 -(\pi k)^2\big ), \end{aligned}$$

which implies

$$\begin{aligned} \varDelta _{kn}=\big ((\pi n)^2 -(\pi k)^2\big ) \big (e_0\big ((\pi n)^2 + (\pi k)^2\big ) - a \big ). \end{aligned}$$

For \(k \ne n\) one has \(|(\pi n)^2 -(\pi k)^2|\ge 3\pi ^2\) and \((\pi n)^2 + (\pi k)^2> 4\pi ^2\). Therefore, if, for example,

$$\begin{aligned} 4 \pi ^2 e_0 > a \end{aligned}$$

then condition (7) is satisfied. The resonance condition reads

$$\begin{aligned} \varDelta _{kn}=O(\delta _b), \end{aligned}$$
(11)

for some \(k \ne n\). It is satisfies if

$$\begin{aligned} |(\pi n)^2 + (\pi k)^2 - a e_0^{-1}| =O(\delta _b) \end{aligned}$$
(12)

for some positive integers n and k, \(n \ne k\). Note that we are dealing here with an internal resonance. The number of modes involved in the resonance depends on \(a e_0^{-1}\pi ^{-2}\) and on the right-hand side of Eq. (12). It might be possible that not 2 but 4 or even more modes satisfy Eq. (12), that is, the sum of two squares can be in more than one way equal to the same number. In the next subsection we will consider an internal resonance involving only two modes.

3.4 Perturbed eigenfunctions in case of an internal resonance

Let us represent the eigenfunctions \(\psi _n\) and \(\psi _k\) as (see [31])

$$\begin{aligned}&\psi _n=c_{1n} \psi _n^{(0)} + c_{1k} \psi _k^{(0)} + \delta _b \phi _n,\\&\quad \psi _k=c_{2n} \psi _n^{(0)} + c_{2k} \psi _k^{(0)} + \delta _b \phi _k, \end{aligned}$$

where \(c_{1n}, c_{1k}, c_{2n}, c_{2k}\) are unknown coefficients and the correction functions \(\phi _{n}\) and \(\phi _{k}\) are orthogonal in \(L_2[0,1]\) to both \(\psi _n^{(0)}\), and \( \psi _k^{(0)}\). Then, the unknown coefficients and corrections to the eigenvalues can be found from the following \(2\times 2\) linear algebraic eigenvalue problem:

$$\begin{aligned} R C= {\bar{\lambda }} C, \end{aligned}$$
(13)

where

$$\begin{aligned} R_{11}= & {} \int _0^1 b_1 (x)|\psi _n^{(0)}(x) |^2\mathrm{d}x, \\ R_{22}= & {} \int _0^1 b_1 (x) |\psi _k^{(0)}(x)|^2 (x)\mathrm{d}x, \\ R_{12}= & {} \int _0^1 b_1(x) \psi _n^{(0)}(x) \psi _k^{(0)}(x)\mathrm{d}x, \end{aligned}$$

and \(C=(C_1, C_2)^{tr}\), where tr stands for the transposed. The first eigenvalue \({\bar{\lambda }}_{+}\) of (13) gives us the perturbation of \(\lambda _n^{(0)}\), and the second one \({\bar{\lambda }}_{-}\) is the perturbation of \(\lambda _k^{(0)}\):

$$\begin{aligned} \lambda _n= & {} \lambda _n^{(0)} + \delta _b {\bar{\lambda }}_{+} +O(\delta _b^2), \end{aligned}$$
(14)
$$\begin{aligned} \lambda _k= & {} \lambda _k^{(0)} + \delta _b {\bar{\lambda }}_{-} +O(\delta _b^2). \end{aligned}$$
(15)

The \({\bar{\lambda }}_{\pm }\) are given by

$$\begin{aligned} {\bar{\lambda }}_{\pm }=\frac{ R_{11} + R_{22} \pm \sqrt{(R_{11} -R_{22})^2 + 4R_{12}^2}}{2}. \end{aligned}$$
(16)

The two independent eigenvectors \((C_{1}^{(l)}, C_{2}^{(l)})^{tr}\), \(l=1,2\) are related to the coefficients \(c_{ln}\) in the following way:

$$\begin{aligned} c_{ln}= C_{1}^{(l)}, \quad c_{lk}= C_{2}^{(l)}, \quad l=1,2. \end{aligned}$$
(17)

Note that the eigenfunctions are orthonormal, therefore the coefficients \(c_{ln}\) are also orthonormal. It is convenient to write down \(c_{ln}\) as

$$\begin{aligned} c_{1n}= & {} \sin (\xi _{nk}), \quad c_{2n}=\cos (\xi _{nk}),\\ c_{1k}= & {} -\cos (\xi _{nk}), \quad c_{2k}=\sin (\xi _{nk}), \end{aligned}$$

where \(\xi _{nk}\) is a resonance angle. Under condition \({{\bar{R}}}_{12}\ne 0\) for that angle one has the following relation:

$$\begin{aligned} \tan \xi _{nk}= & {} \frac{ R_{12}}{{\bar{R}}_{12}},\nonumber \\ 2{{\bar{R}}}_{12}= & {} (R_{11} - R_{22}) - \sqrt{(R_{11} -R_{22})^2 + 4R_{12}^2}. \end{aligned}$$
(18)

If \({{\bar{R}}}_{12}=0\), then \(R_{12}=0\) and we have no resonance, because system (13) is diagonal and the interaction between the modes n and k is of order \(O(\delta _b^2)\). The parameter \(\xi _{nk}\) plays an important role, as will be shown in the coming subsection.

3.5 A new resonance effect: harmonics disappear for a long time

We suppose that damping is absent, i.e., \(\epsilon =0\). Suppose that the initial data are given by a simple harmonic function, for example

$$\begin{aligned} u(x, 0)=A\sin (\pi n x), \quad u_t(x,0)=0, \end{aligned}$$

where A is an amplitude (for a more general case, where \(u_t(x,0) =B \sin (\pi nx)\), our results are similar). We have

$$\begin{aligned} u(x,t)= A \sum _{l=1}^{\infty } C_{nl} \cos (\omega _l t) \psi _l(x), \end{aligned}$$
(19)

where

$$\begin{aligned} C_{nl}=\int _0^1 \sin (\pi n x) \psi _l(x) \mathrm{d}x. \end{aligned}$$

The Fourier coefficient \( u_n(t)=\int _0^1 u(x,t) \sin (\pi n x) \mathrm{d}x \) evolves in time according to

$$\begin{aligned} u_n(t)= A \sum _{l=1}^{\infty } C_{nl}^2 \cos (\omega _l t). \end{aligned}$$
(20)

In the non-resonance case the main contribution in the sum (20) is given by the single term \(l=n\) and \( u_n(t)\) oscillates with frequency \(\omega _n\). The effect of the elastic foundation produces small contributions of order \(O(\delta _b)\) in all other terms in (20) for \( l \ne n\). In the resonance case [that is , when (12) is satisfied for \(a e_0^{-1}=5\pi ^{2}\)] we obtain that the two modes \(\psi _n\) and \(\psi _k\) give rise to significant contribution in the sum of (20), that is,

$$\begin{aligned} u_n(t)= A \big ( C_{nn}^2 \cos (\omega _n t) + C_{nk}^2 \cos (\omega _k t) +O(\delta _b)\big ), \end{aligned}$$
(21)

and

$$\begin{aligned} u_k(t)= A \big ( C_{kk}^2 \cos (\omega _k t) + C_{nk}^2 \cos (\omega _n t) +O(\delta _b)\big ). \end{aligned}$$
(22)

From the previous subsection it follows that

$$\begin{aligned} C_{nn}=c_{1n}=\sin \xi _{nk}, \quad C_{nk}=c_{1k}=-\cos \xi _{nk} \end{aligned}$$

and, as a result, we have

$$\begin{aligned} u_n(t)= & {} A \Big ( (\sin \xi _{nk})^2 \cos (\omega _n t) + (\cos \xi _{nk})^2 \cos (\omega _k t) \nonumber \\&+O(\delta _b)\Big ). \end{aligned}$$
(23)

Note that in this resonance case \(\omega _n -\omega _k=O(\delta _b)\). By introducing the mean frequency and the deviation by

$$\begin{aligned} {\bar{\omega }} =\frac{\omega _n + \omega _k}{2}, \quad {\tilde{\omega }} =\frac{\omega _n - \omega _k}{2} \end{aligned}$$

respectively, we can transform Eq. (23) as follows:

$$\begin{aligned} u_n(t)= & {} A \Big (\cos ({\bar{\omega }} t) \sin ({\tilde{\omega }} t) \nonumber \\&+ \big (\sin ^2 \xi _{nk} - \cos ^2 \xi _{nk}\big )\cos ({\bar{\omega }} t) \cos ({\tilde{\omega }} t) \Big ) +O(\delta _b). \end{aligned}$$
(24)

Using (18) one obtains

$$\begin{aligned} u_n(t)= & {} A \Big (\cos ({\bar{\omega }} t) \sin ({\tilde{\omega }} t) \nonumber \\&+ \frac{{{\bar{R}}}_{12}^2 - R_{12}^2}{R_{12}^2 + {{\bar{R}}}_{12}^2} \cos ({\bar{\omega }} t) \cos ({\tilde{\omega }} t) \Big ) +O(\delta _b). \end{aligned}$$
(25)

Relation (25) shows a new internal resonance effect only possible for a beam on a Winkler foundation. It is possible that \( u_n(t)\) is close to \(O(\delta _b)\) for large times and we observe beats (see the plot in Fig. 1). The maximal beat effect, when \(u_n(t)\) is close to 0 within a time interval, of \(O(\delta _b^{-1})\), arises when \(\xi _{nk} \approx \pi /4\). Note that \(\xi _{nk} \approx \pi /4\) for \(|R_{11}-R_{22}|=O(\delta _b)\), i.e., for

$$\begin{aligned} \int _0^1 b_1(x) \sin ^2(\pi n x) \mathrm{d}x-\int _0^1 b_1(x) \sin ^2(\pi k x) \mathrm{d}x =O(\delta _b). \end{aligned}$$
(26)

We will use this condition (26) to have maximal beat effect in Sect. 4.

Fig. 1
figure 1

This plot shows the time evolution of the Fourier coefficient \(u_n(t)\) in the solution u(x, t) when \(u(x,0)=A\sin (\pi n x)\) and \(u_t(x,0)=0\) in the resonance case. The plot is computed by using formula (24), where \(n=1, k=2\) and the parameters are \({\bar{\omega }}=1, \xi _{nk}=\pi /4, {\tilde{\omega }}=0.1\) and \(A=2\). The correction \(O(\delta _b)\) is approximated by the first three nonzero terms in the Fourier series with small coefficients. We see beats induced by an internal resonance.

Full size image

3.6 Strong Winkler foundation

In the case where the perturbation \(b_1\) is not small, we can apply the Rayleigh variation principle: the minimal eigenvalue \(\lambda \) can be found as the minimum of the Rayleigh quotient Q over all test functions, that is,

$$\begin{aligned} \lambda _{\min }= \min _{\psi } {Q[\psi ]}, \quad Q=\frac{V[\psi ]}{||\psi ||^2}, \quad \end{aligned}$$
(27)

where the test functions \(\psi \in C^2[0, 1]\) should satisfy the boundary conditions (4), and V is given by

$$\begin{aligned} V[\psi ]=e_0 ||\psi _{xx}||^2 - a ||\psi _{x}||^2 + \int _0^1 b(x) \psi (x)^2 \mathrm{d}x. \end{aligned}$$

We will use this approach with test functions to study so-called localized modes, in Sect. 4. In [28] a similar problem is studied, but with a different type of method. The advantage of the use of test functions is that we can avoid complicated integral equations.

4 Beam buckling load for localized and non-localized modes

In this section we will consider different particular types of inhomogeneities in the elastic foundation of the beam, and we will approximately determine for which parameters values buckling occurs.

Example A

a regular harmonic perturbation in the foundation. Consider a smooth perturbation \(b_1(x)= \sin (\gamma x)\), where \(\gamma >0\). Then, (8) gives

$$\begin{aligned} \varDelta \lambda _n= \lambda _n-\lambda _n^{(0)} \approx - \delta _b \frac{(1- \cos \gamma ) (2\pi n)^2}{\gamma (\gamma ^2 - ( 2\pi n)^2)} \end{aligned}$$
(28)

for \(\gamma \ne 2\pi n\). For \(\gamma =2\pi n+ {\tilde{\gamma }}\), where \({\tilde{\gamma }} =O(\delta _b)\), we can simplify the expression in the right-hand side of (28) by using \(1- \cos (\gamma )= {\tilde{\gamma }} ^2/2 +O({\tilde{\gamma }}^4)\) , yielding

$$\begin{aligned} \varDelta \lambda _n \approx - \delta _b {\tilde{\gamma }}/4. \end{aligned}$$
(29)

For the buckling load \(a_c(n)\) corresponding to the stability loss for the n-th mode we approximately obtain

$$\begin{aligned} a_c(n)= (\pi n)^{-2} (b_0 + \varDelta \lambda _n) + e_0 (\pi n)^{2}, \end{aligned}$$
(30)

where \(\varDelta \lambda _n\) is defined by (28) for \(\gamma \ne 2\pi n\) or by (29) for \(\gamma =2\pi n+O(\delta _b)\).

Example B

the elastic foundation has point inhomogeneities in the form of concentrated springs with negative or positive stiffness.

Let us consider a case which corresponds to point inhomogeneities in the form of concentrated springs with negative or positive stiffness. We will consider an irregular perturbation \(b_1(x)=\sum _{j=1}^{n_d} \beta _j \delta (x -x_j)\), where \(\delta \) stands for the Dirac delta function and \(\beta _j\) are coefficients. For this case we will also observe a “space resonance” effect. Suppose that \(x_j=(j-1/2)r+{\tilde{x}}\), where \(j=1,\ldots , n_d\) and \(r=1/n_d\). The parameter \(\beta _j\) may have negative, or positive values. So, we assume that all point inclusions of the elastic foundation are equidistant and the parameter \({\tilde{x}} \in (0, 1/2n_d)\) describes a shift of the inclusion coordinates with respect to the beam edges. Under these assumptions, all inclusions lie within (0, 1). When \({\tilde{x}}=0\) all inclusions are located symmetrically with respect to the beam edges. Moreover, we suppose that \(n_d \ll O(\delta _b^{-1})\) .

Let \(\beta _j=\beta \) for all j. Then, there two cases have to be considered. In the first case, where \(k=n/n_d\) is not an integer, the perturbation of \(\lambda _n\) depends on the number of inclusions only and does not depend on the mode number n and the shift \({\tilde{x}}\) (see “Appendix”):

$$\begin{aligned} \varDelta \lambda _n =\lambda _n-\lambda _n^{(0)} = \beta \delta _b n_d + O({\delta _b}^2). \end{aligned}$$
(31)

This means that the contributions of the inclusions mutually cancel each other.

In the second case, when \(k=n/n_d\) is an integer, we obtain (see “Appendix”)

$$\begin{aligned} \varDelta \lambda _n = \beta \delta _b n_d \big (1- (-1)^k \cos (2\pi n {\tilde{x}})\big ) + O({\delta _b}^2). \end{aligned}$$
(32)

The last relation shows that a “space resonance” effect occurs, which depends on \({\tilde{x}}\). This effect is stronger if \({\tilde{x}}=0\) (see Fig. 2). Then, for the odd k we have \(\varDelta \lambda _n=2\beta \delta _b n_d\), i.e., the effect is doubled with respect to the case of a non-integer k, while for even k one obtains \(\varDelta \lambda _n=0\). In both cases we can use (30) for \(a_c(n)\). So, we obtain a higher value for \(\varDelta \lambda _n\) , if \(n/n_d\) is an odd integer and \({\tilde{x}}=0\). The results are illustrated in Figs. 2, 3, 4 and 5. Figure 2 shows the effect of the “space resonances,” where we see that the number of peaks decreases with the number of inclusions \(n_d\). Figure 3 shows the dependence of the buckling load on \(b_0\) for a given \(\delta _b\) and different inclusion numbers. In Fig. 4 we observe a weak effect of irregularity in the dashed and solid curves for the dependence of \(a_c\) on \(b_0\). This effect can be explained if we take into account that the number of the critical modes is an integer depending on \(b_0\), and which makes jumps for some \(b_0\). The effect is induced by the inclusions (\(n_d=3)\) and the “space resonances.” The dashed curve corresponds to the inhomogeneous elastic foundation (\(\delta _b=0.1\)), the star curve shows the case of the foundation without inclusions, and the solid one describes the case for a negative \(\delta _b\) (\(\delta _b=-0.1\)). We can observe those jumps in Fig. 5, and we see that the jumps exactly correspond to the small irregularities in Fig. 4. In Fig. 5 we also see the most essential effect of the inclusions in the elastic foundation: the number of the critical modes is different for the cases with inclusions and without inclusions, and it sharply depends on the sign of \(\delta _b\). For negative \(\delta _b\) we obtain more jumps. Note that, as it follows from the Central Limit Theorem that for random \(x_j\) and \(n_d \gg 1\), when we are dealing with a random perturbation, the effect of the inclusions is almost absent and \({\delta _b}^{-1}(\lambda _n-\lambda _n^{(0)})\) is small. In fact, then the contributions of different inclusions mutually cancel each other. By the relations obtained above we can analyze condition (26) leading to the maximal beats in the examples A and B. In example A we use the results of “Appendix,” and the expressions (28) and (29). Then, after elementary algebra, we obtain that (26) is satisfied for \(\gamma =2\pi l +O(\delta _b)\), where l is a positive integer (this integer can be equal to n or k).

Fig. 2
figure 2

The dependence of \(\varDelta \lambda _n\) up to \(O(\delta _b^2)\) on the mode number n. The parameters are \(e_0=10^{-6}\), \(L=1\), the maximal mode number \(N_{\max }=50\), \(b_0=1\), \(\delta _b=0.01\), and \(\beta =1\). The shift \({\tilde{x}}=0\). The curves correspond to the cases, where the number of point inclusions are \(n_d=3\) and \(n_d=10\)

Full size image
Fig. 3
figure 3

The dependence of the buckling load \(a_c\) on the parameter \(b_0\). The value \(\delta _b=0.1\) is fixed. The parameters are \(e_0=10^{-6}\), \(L=1\), the maximal mode number \(N_{\max }=50\), \({\tilde{x}}=0\), and the number of inclusions \(n_d=1, 3\), 10

Full size image
Fig. 4
figure 4

The dependence of the buckling load \(a_c\) on the parameter \(b_0\). We see a weak irregularities for the dashed and solid curves. The effect is induced by the inclusions (\(n_d=3)\) and the “space resonances.” The dashed curve corresponds to the inhomogeneous elastic foundation (\(\delta _b=0.1\)), the star curve shows the case of the foundation without inclusions, and the solid one describes the case of the negative \(\delta _b\) (\(\delta _b=-0.1\)). The parameters are \(e_0=10^{-6}\), \(L=1\), the maximal mode number \(N_{\max }=50\), \(b_0=1\), \({\tilde{x}}=0\)

Full size image
Fig. 5
figure 5

The dependence of the critical mode index \(n_c\) on the parameter \(b_0\) for three different values of \(\delta _b\): \(\delta _b=0.1\), \(\delta _b=0\) and \(\delta _b=-0.1\). The parameters are \(e_0=10^{-6}\), \(L=1\), the maximal possible mode number \(N_{\max }=50\), (i.e., we take the mode indices n from the set \(\{1, \ldots , N_{\max } \}\)), \(b_0=1\), \({\tilde{x}}=0\) and the number of inclusions \(n_d=3\)

Full size image

In example B, expression (31) shows that (26) is satisfied for all non-integers \(n/n_d\) and \(k/n_d\), i.e., periodic delta-like inclusions practically always induces maximal beats for internal resonances.

Finally, in case of a point inhomogeneity in the form of a concentrated spring with a negative stiffness, we note that localized modes can exist [29]. We discuss these effects related to localized modes in the following example.

Example C

Rectangular inclusions in the elastic foundation.

If b(x) is piecewise constant then so-called localized modes can occur. These modes are concentrated at inhomogeneities of the elastic foundation. These modes can also lead to beam buckling. To simplify the problem we consider the following, simple coefficient b(x), which depends on the positive parameters \(x_0, b_0, b_{\min } < b_0\), and d:

$$\begin{aligned} b(x)= & {} b_0, \quad |x-x_0| > d,\\ b(x)= & {} b_{\min }, \quad |x-x_0| \le d, \end{aligned}$$

where \(x_0 > d\) and \(x_0+d <1\). Let us denote by \(\varDelta b=b_0 -b_{\min }\) the “depth” of the inclusion. The quantity 2d is its width. If \(d \ll 1\) then so-called localized modes concentrated at \(x=x_0\) and exponentially decreasing in \(|x-x_0|\) exist [28]. These eigenfunctions are studied in [28], where also their influence on the Euler instability is considered. For \(b=const\), \(e_0 \ll 1\), and for compression force a(x) depending on x, the localized modes can be investigated by the WKB method, see [26]. Note that in the limit \(d \rightarrow 0\) we obtain a point inclusion as considered in [29]. For \(d \ll 1\) the unperturbed non-localized eigenfunctions \(\psi _{n}\) have the form \(\sqrt{2} \sin (\pi n x)\). Then, we have the asymptotics, which are valid up to corrections of the order \(d^2\):

$$\begin{aligned} V[\psi _n] \approx \varPhi (z_n)=e_0 z_n^4 - a z_n^2 + b_0 - 2d \varDelta b \sin ^2(\pi n x_0) , \end{aligned}$$
(33)

where \(\quad z_n= \pi n\), and \( V[\psi _n] \) is defined in (27). If a buckling mode has number n, then

$$\begin{aligned} a_{c, nloc}(n) \approx e_0 (\pi n)^2 + (\pi n)^{-2}\big (b_0 - 2d \varDelta b \sin ^2(\pi n x_0) \big ). \end{aligned}$$
(34)

For small \(e_0\) we can consider \(z_n\) as a real valued parameter (since the minimum of \(\varPhi (z_n)\) is obtained for large n). Then we minimize \(\varPhi (z)\) with respect to z, which implies:

$$\begin{aligned} \pi n \approx \sqrt{a/2e_0}, \end{aligned}$$

and

$$\begin{aligned} \quad \min V \approx - \frac{a^2}{4e_0} + b_0 - 2d \varDelta b \sin ^2(\pi nx_0). \end{aligned}$$

The value of \(a_c\) can be found by the condition that \(\min V =0\) , which implies:

$$\begin{aligned} a_{c, nloc} \approx 2 \sqrt{ e_0 \big ( b_0 - 2 d\varDelta b \sin ^2(\pi n_c x_0) \big )}, \end{aligned}$$
(35)

where the critical mode number is

$$\begin{aligned} n_c = [(b_0/e_0)^{1/4} \pi ^{-1}] , \end{aligned}$$

and where [z] stands for the integer nearest to z.

Note that under our assumptions \(\varDelta b < b_0\). Relation (35) is the expression for the buckling load for an infinite beam (see [29]). We see that the onset of the Euler instability can give rise to a high-frequency mode. Moreover, we see an effect of “space resonance”: that is, the magnitude of \(a_c\), induced by the inclusion, depends on the inclusion location. In the “space resonance” case, the change of \(a_c\) is proportional to the inclusion number \(n_c\). Let us now consider the case of the localized modes. We suppose that \(e_0\) is a small parameter. In the Rayleigh quotient [see (27)] we substitute the test functions, which are well localized at \(x_0\), for example,

$$\begin{aligned} \psi _\text {test, loc} =r^{-1} \theta ((x-x_0)/r), \end{aligned}$$

where \(\theta (y)\) is a well localized function at \(y=0\) (for instance, we can take \(\theta (y)=\exp (- y^{\gamma }), \quad \gamma \ge 2\)). Here, the parameter r defines the characteristic radius of the test mode localization. For small \(r\ll d\) one has:

$$\begin{aligned} V [\psi _\text {test, loc}]\approx & {} \varPhi (r)=c_0 e_0 r^{-5} - c_1 a r^{-3} + c_2 b_{\min }r^{-1} \nonumber \\&+ O(1) \end{aligned}$$
(36)

for \(r\rightarrow 0\), where the positive constants \(c_1, c_2\), and \(c_0\) are independent of r, and are given by

$$\begin{aligned} c_0= & {} \int _{-\infty }^{\infty } \Bigl (\frac{\mathrm{d}z^2\theta (z)}{\mathrm{d}z^2}\Bigr )^2 \mathrm{d}z, \quad c_1=\int _{-\infty }^{\infty }\Bigl (\frac{\mathrm{d}z\theta (z)}{\mathrm{d}z}\Bigr )^2 \mathrm{d}z, \\ c_2= & {} \int _{-\infty }^{\infty } \theta ^2 \mathrm{d}z. \end{aligned}$$

The same expression (36) can be obtained for an infinite beam. Thus

$$\begin{aligned} a_{c, loc} =2 \sqrt{ e_0 \ b_{\min } } (1 + o(1)), \end{aligned}$$
(37)

where \(o(1) \rightarrow 0\) as \(e_0 \rightarrow 0\). Note that our approach does not need a solution of a complicated integral equation and numerical simulations, and it can be extended to study more general forms for the functions b(x) . By comparing the relations (35) and (37), we observe that for \( b_{\min } \ll b_0\) the onset of the Euler instability starts at the localized modes. If a localized mode exists then the buckling load \(a_c\), which corresponds to this mode is always less than the buckling load for a non-localized mode. To see this, we have to compare relations (35) and (37) and take into account that \(b_0 > b_{\min }\). Note that the Euler beam buckling in presence of localized modes is also studied in [28] by another method.

5 The influence of damping on buckling

Let us consider the influence of damping on the behavior of the beam near buckling. Solutions of Eq. (2) can be expressed in an eigenfunction expansion by using the orthonormal eigenfunctions \(\psi _n\):

$$\begin{aligned} u(x,t)=\sum _{n \in {{\mathbb {N}}}} X_n(t)\psi _n(x), \end{aligned}$$
(38)

where \( X_n(t)\) are unknown functions. Then, the coefficients \( X_n(t)\) have to satisfy:

$$\begin{aligned} \frac{\mathrm{d}z^2 X_n}{\mathrm{d}t^2} + \omega _n^2 X_n = - \epsilon \left\langle D\left[ \sum _{k=1}^{\infty } X_k(t) \psi _k(x)\right] , \psi _n(x) \right\rangle . \end{aligned}$$
(39)

Assume that the Euler–Bernoulli beam becomes unstable for a certain \(a=a_c\). Let us consider values of a close to this critical value \(a_c\). Let us denote by Y the magnitude of the mode corresponding to \(\psi _{n_*}=\varPsi \), which is close to instability, and \(n_*\) denotes its number. Note that it is possible that \(n_* \ne 1\), as was shown in the previous section. For small \(|a-a_c|\) the frequency \(\omega _{n^*}\) is small and close to zero for a certain \(n_*\). For that reason we introduce a new small parameter \(\omega _{n^*}=\mu >0\). Obviously the amplitude \(Y= X_{n_*}\) of the mode for which the beam lost its stability, is slowly varying in time.

We suppose that \(\delta _b, \epsilon \) and \(\mu \) are small parameters which are related in the following way:

$$\begin{aligned} \delta _b=\epsilon {\tilde{\delta }}_b, \quad \mu =\epsilon \mu _0, \end{aligned}$$
(40)

where \(\mu _0\) and \({\tilde{\delta }}_b\) are positive parameters, which are independent of the small, positive parameter \(\epsilon \). Further, to solve the equations describing the time evolution of the stable (fast) modes \(X_n\) with \(n \ne n_*\), we use a multiple time-scale perturbation method with fast time t and slow time \(T=\epsilon t\) . Note that the parameter \(\delta _b\) is not involved directly in the suggested multiple time-scale approach, but it is present in the equation for \(a_c(\delta _b)\) and in the orthonormal eigenfunctions \(\psi _n(x)\), which depend on \(\delta _b\) according to the expressions as obtained in the previous section. To solve (39), we assume that \(Y=Y(t, T)\) and \(X_n=X_n(t, T)\) for \(n \ne n_*\). We obtain then the following system of equations for \(Y={X_n}_*\) and the fast modes \(X_n(t)\) with \(n \ne n_*\) :

$$\begin{aligned}&\frac{\partial ^2 Y}{\partial t^2} + {{\mathcal {B}}}Y + \epsilon ^2 \mu _0^2 Y + \epsilon \langle D[Y \varPsi ], \varPsi \rangle \nonumber \\&\quad =-\epsilon \sum _{n=1,|n| \ne n_*}^{\infty } \langle D[X_n(t,T) \psi _n ], \varPsi \rangle , \end{aligned}$$
(41)
$$\begin{aligned}&\frac{\partial ^2 X_n(t,T)}{\partial t^2} + {{\mathcal {B}}} X_n + \omega _n^2 X_n \nonumber \\&\quad = -\epsilon \left\langle D\Biggl [Y \varPsi + \sum _{k=1,|k| \ne n_* }^{\infty } X_k(t,T) \psi _k \Biggr ], \psi _n\right\rangle , \end{aligned}$$
(42)

where \({{\mathcal {B}}} \) is the following differential operator:

$$\begin{aligned} {{\mathcal {B}}} =2\epsilon \frac{\partial ^2 }{\partial t \partial T} + \epsilon ^2 \frac{\partial ^2 }{\partial T^2}. \end{aligned}$$

The \(+\) sign in (41) before the term \(\epsilon ^2 \mu _0^2 Y\) in the left-hand side of (41) corresponds to a weakly stable situation, when \(a < a_c\). To study the weakly unstable case, when \(a > a_c\), we should put the − sign in the left-hand side of (41) before the term \(\epsilon ^2 \mu _0^2 Y\). To fix ideas, we choose the \(+\) sign.

To find asymptotic approximations of the solutions of (41) and (42), we should define the form of the operator D. To this end, let us discuss briefly some damping models. We will use damping models which are for instance suggested in [25]. In this paper, two cases will be distinguished: external and internal damping. For the external damping we have

$$\begin{aligned} D[u]= \int _0^1 \int _0^{t} h_{td} (x-\xi , t-\tau ) u_t(\xi , \tau ) \mathrm{d}z\xi \mathrm{d}z\tau , \end{aligned}$$
(43)

and for internal damping

$$\begin{aligned} D[u]= \int _0^1 \int _0^{t} h_{td} (x-\xi , t-\tau ) L_s u_t(\xi , \tau ) \mathrm{d}z\xi \mathrm{d}z\tau , \end{aligned}$$
(44)

where \(L_s\) is a linear operator, \( L_s u = \frac{\partial ^4 u}{\partial ^4 \xi }. \) The simple damping (air damping) and Kelvin–Voigt damping (referred to as the SD and KV cases) occur if \(h_{td}=\delta (t-\tau ) \delta (x-\xi )\) in the relations (43) and (44), respectively. For space hysteresis induced damping (that we will be referred to as the SH case) we will take the kernel \(h_{sd}\) [25] as a function which is well localized in \(|x-\xi |\), for example, it may be a Gaussian function:

$$\begin{aligned} h_{td}=h_{sd}= \delta (t-\tau ) (\sigma \sqrt{\pi })^{-1} \exp \big (- (x-\xi )^2 /2\sigma ^2 \big ), \end{aligned}$$

where \(\sigma \) a positive parameter with \(\sigma \ll 1\). From (41) it follows that the first order part of Y(t, T) depends only on T. And so Eq. (42) has the following asymptotic solution:

$$\begin{aligned} X_n (t, T) = {{\bar{X}}}_n(T) + {\tilde{X}}_n(t, T), \end{aligned}$$
(45)

where

$$\begin{aligned} {\tilde{X}}_n(t, T)=A_n(T) \sin (\omega _n t) + B_n(T) \cos (\omega _n t) \end{aligned}$$

and

$$\begin{aligned} {{\bar{X}}}_n (T) = -\omega _n^{-2} \langle D[ Y \varPsi ], \psi _n\rangle , \end{aligned}$$
(46)

where \(A_n, B_n\) are functions of the slow time T, which can be found by the standard procedure [32,33,34,35,36]. We suppose that \(A_n, B_n=O(1)\). Then, the terms \({\tilde{X}}_n(t, T)\) produce small contributions in Y with respect to the terms \({{\bar{X}}}_n (T)\). To see this, let us consider Eq. (41) with \(X_n={\tilde{X}}_n(t, T)\) in the right-hand side. The corresponding asymptotic solution \({\tilde{Y}}(t, T)\) involves the harmonics \(\sin (\omega _n t), \ \cos (\omega _n t)\) with \(n \ne n_*\) and coefficients depending on T. This part of the solution has order \(O(\epsilon )\) because the right-hand side of Eq. (41) is proportional to \(\epsilon \). Now let us consider (41) with \(X_n={{\bar{X}}}_n( T)\) in the right-hand side, which is proportional to \(\epsilon ^2\) and depends only on T. The corresponding solution \({{\bar{Y}}}\) is defined by

$$\begin{aligned}&\frac{\mathrm{d}z^2 {{\bar{Y}}}}{dT^2} + \mu _0^2 {{\bar{Y}}} + \epsilon ^{-1} \langle D[{{\bar{Y}}} \varPsi ], \varPsi \rangle \nonumber \\&\quad =- \epsilon ^{-1} \sum _{|n| \ne n_*}^{\infty } \langle D[{{\bar{X}}}_n(T) \psi _n ], \varPsi \rangle , \end{aligned}$$
(47)

and, in general, has order O(1). Therefore, we conclude that \(X_n ={{\bar{X}}}_n(T)\) and \({{\bar{Y}}}(T)\) satisfy Eqs. (41) and (42) up to terms of higher orders in \(\epsilon \). In fact, let us substitute \(Y={{\bar{Y}}}\) into Eq. (42), giving

$$\begin{aligned}&\frac{\partial ^2 X_n(t)}{\partial t^2} + {{\mathcal {B}}} X_n + \omega _n^2 X_n \\&\quad = -\epsilon \left\langle D\Biggl [{{\bar{Y}}} \varPsi + \sum _{|k| \ne n_* }^{\infty } X_k(t) \psi _k \Biggr ], \psi _n\right\rangle . \end{aligned}$$

The function \(X_n(t, T)={{\bar{X}}}_n(T)\) satisfies the last equation up to

$$\begin{aligned} \epsilon \left\langle D\Biggl [ \sum _{|k| \ne n_* }^{\infty } {\bar{X}}_k(T) \psi _k \Biggr ], \psi _n \right\rangle =O(\epsilon ^2), \end{aligned}$$

and since \({{\bar{Y}}}\) is a function of T only, it follows that the function \({{\bar{X}}}_n(T)\) is of order \(\epsilon \).

We substitute the relations (45) and (46) into Eq. (41). Since Y(T) is slowly varying in time, it follows for the cases SD, KV and SH that

$$\begin{aligned} {{\bar{X}}}_n (T) =\alpha _n \frac{\mathrm{d}z {{\bar{Y}}}}{dT}, \end{aligned}$$

where the constants \(\alpha _n\) determine an interaction force between the unstable mode and the n-th stable one, and have the form:

$$\begin{aligned} \alpha _n= & {} \int _0^1 \varPsi (x) \psi _n(x) \mathrm{dn}x, \end{aligned}$$
(48)
$$\begin{aligned} \alpha _n= & {} \int _0^1 \frac{\mathrm{d}^4\varPsi }{\mathrm{d}x^4} \psi _n(x) \mathrm{d}x, \end{aligned}$$
(49)
$$\begin{aligned} \alpha _n= & {} \int _0^1 \int _0^1 \frac{\mathrm{d}z^4 \varPsi (\xi )}{\mathrm{d}z\xi ^4} \psi _n(x) h_{sd}(x-\xi ) \mathrm{d}x \mathrm{d}z\xi , \end{aligned}$$
(50)

and

$$\begin{aligned} \alpha _n = \int _0^1 \int _0^1 \varPsi (\xi ) h_{sd}(x-\xi ) \psi _n(x) \mathrm{d}x \end{aligned}$$
(51)

for SD, KV damping, internal space hysteresis damping, and for external space hysteresis damping, respectively. Taking into account that in (47)

$$\begin{aligned} \langle D[{{\bar{X}}}_n(T) \psi _n ], \varPsi \rangle \approx - \omega _n^{-2} \alpha _n^2 \frac{\mathrm{d}z^2{{\bar{Y}}}}{dT^2} \end{aligned}$$

which follows from \({{\bar{X}}}_n (T) =\alpha _n \frac{\mathrm{d}z {{\bar{Y}}}}{dT}\), and from the arguments given after (47), and so we finally obtain::

$$\begin{aligned} (1+ m_\text {add}) \frac{\mathrm{d}z^2 {{\bar{Y}}}}{dT^2} + \mu _0^2 {{\bar{Y}}} + \kappa \frac{\mathrm{d}z {{\bar{Y}}}}{dT} =0, \end{aligned}$$
(52)

where the relative “added mass” \(m_\text {add}\) is

$$\begin{aligned} m_\text {add} = - \sum _{n \ne n_*}^{\infty } \omega _n^{-2} \alpha _n^2, \end{aligned}$$
(53)

where

$$\begin{aligned} \kappa= & {} 1, \end{aligned}$$
(54)
$$\begin{aligned} \kappa= & {} \int _0^1 \frac{\mathrm{d}^4\varPsi }{\mathrm{d}x^4} \varPsi (x) \mathrm{d}x, \end{aligned}$$
(55)
$$\begin{aligned} \kappa= & {} \int _0^1 \int _0^1 \frac{\mathrm{d}z^4 \varPsi (\xi )}{\mathrm{d}\xi ^4} \varPsi (x) h_{sd}(x-\xi ) \mathrm{d}x \mathrm{d}z\xi , \end{aligned}$$
(56)

and

$$\begin{aligned} \kappa = \int _0^1 \int _0^1 \varPsi (\xi ) h_{sd}(x-\xi ) \varPsi (x) \mathrm{d}x \end{aligned}$$
(57)

for SD, KV damping, internal space hysteresis damping, and for external space hysteresis damping, respectively.

So, we conclude that a weak interaction between modes produces an “added mass” effect. This “added mass” effect diminishes when the natural frequencies increase, and is proportional to the square of the damping magnitude. However, it does not depend on the parameter \(\epsilon _0\), which determines the deviation between a and the critical value \(a_c\). In fact, this effect is proportional to the coefficient \(\alpha _n\), which is defined by the relations (48), (49) (50) and (51). Now we can conclude the following:

  1. (a)

    in the SD case \(\alpha _n=0\) and the added mass effect is absent at the order 1 level;

  2. (b)

    in the KV damping case \(\alpha _n=0\) for beams on homogeneous elastic foundations, and \(\alpha _n \) is proportional to \(\delta _b\) for inhomogeneous elastic foundations and hinged supported beams;

  3. (c)

    similar results can be obtained for the external and internal SH cases. Let \(\sigma \ll 1\). Then, we conclude that for simply supported beams and external or internal space hysteresis damping \(\alpha _n=O(\exp (-1/\sigma ))\) for beams on homogeneous elastic foundations, and \(\alpha _n = O(\delta _b)\) for inhomogeneous elastic foundations.

Note that the modal interaction does not change (up to order \(\epsilon \)) the value of the critical Euler force \(a_c\). Indeed, the value of \(a_c\) depends on parameters involved in the self-adjoint operator \({{\mathcal {L}}}\) since the instability arises when the spectrum of that operator contains 0.

To study damping effects, let us consider Eq. (52). By solving (52) we find that

$$\begin{aligned} {{\bar{Y}}}= c_1 \exp ( \theta _{1} T) + c_1 \exp ( \theta _{2} T), \end{aligned}$$

where

$$\begin{aligned} \theta _{1}= & {} \frac{-\kappa + \sqrt{ \kappa ^2 - 4\mu _0^2(1+ m_\text {add})}}{2(1+ m_\text {add})}, \nonumber \\ \theta _2= & {} \frac{-\kappa - \sqrt{ \kappa ^2 - 4\mu _0^2(1+ m_\text {add})}}{2(1+ m_\text {add})}. \end{aligned}$$
(58)

Here, we have two regimes, the first one corresponds to the case for which the system is relatively far away from an instability, and damping forces dominate (\(\kappa \gg 4\mu _0^2(1+m_\text {add})\)). When \(a < a_c\) it follows from (58) there exists a weakly decreasing mode with \({\text {Re}}\theta _{1} \approx - \mu _0^2\kappa ^{-1}\), and for \(a > a_c\) we also obtain a weakly decreasing mode with \({\text {Re}}\theta _{1} \approx -\mu _0^2 \kappa ^{-1}\). We conclude that in both cases the “added mass” does not affect stability. The second regime corresponds to the case for which the damping forces are relatively small (\(\kappa \ll 4\mu _0^2(1+m_\text {add})\)). For \(a < a_c\) we have \({\text {Re}}(\theta _1) = -\kappa /2(1+ m_\text {add})\) and \({\text {Im}}(\theta _1) \approx \mu _0 \sqrt{1/(1+ m_\text {add})}\). In this case we see that the negative “added mass” effect leads to an increase in the oscillation frequency and moreover that effect decreases the damping. For \(a > a_c\) we see that the negative “added mass” effect also reinforces the instability.

6 Interaction of resonant modes and damping

6.1 Main equation

This section is “served” as an example to show additional existing and interesting phenomena. We restrict ourselves to the case of space hysteresis damping, and we only indicate what might be expected. For sure, there is much to be discovered, but this is outside the scope of this paper. Let \(\psi _n\), \(\psi _k\) be two resonant modes (see Sect. 3.4). For these modes, we introduce the detuning parameter \(\omega _{nk}\) by

$$\begin{aligned} \omega _{nk}=\epsilon ^{-1}(\omega _n - \omega _k). \end{aligned}$$
(59)

Let us consider the most interesting case of space hysteresis damping:

$$\begin{aligned} D u (x,t)=\int _0^1 \int _0^1 G(x, x') u_t(x', t) \mathrm{d}x'. \end{aligned}$$

To describe the mutual evolution of two resonant modes, we are restricting ourselves to the following solution form:

$$\begin{aligned} u(x,t)= X_n (t) \psi _{n} + X_k(t) \psi _k, \end{aligned}$$

where \(X_n, X_k\) are unknown functions. We neglect here effects of added mass and buckling, and we assume that \(\omega _n, \omega _k=O(1)\) (these effects will be interesting to study in a future research). Then, for \(X_n, X_k\) we obtain the following system of equations

$$\begin{aligned} \frac{\mathrm{d}^2 X_n}{\mathrm{d}t^2} + \omega _n^2 X_n= & {} - \epsilon \left( d_{nn} \frac{\mathrm{d}X_n}{\mathrm{d}t} + d_{nk} \frac{\mathrm{d}X_k}{\mathrm{d}t} \right) , \end{aligned}$$
(60)
$$\begin{aligned} \frac{\mathrm{d}^2 X_k}{\mathrm{d}t^2} + \omega _k^2 X_k= & {} - \epsilon \left( d_{kk} \frac{\mathrm{d}X_k}{\mathrm{d}t} + d_{kn} \frac{\mathrm{d}X_n}{\mathrm{d}t} \right) , \end{aligned}$$
(61)

where

$$\begin{aligned} d_{l_1 l_2}= \int _0^1 \int _0^1 G(x, x') \psi _{l_1}(x) \psi _{l_2}(x') \mathrm{d}x \mathrm{d}x', \end{aligned}$$

where \(l_1, l_2\) take values in the set \(\{n, k\}\). This system can be solved analytically, however, to find oscillation magnitudes, it is simpler to seek solutions in the following asymptotic form:

$$\begin{aligned} X_l(t)= A_l(T) \exp (\mathrm {i}\omega _l t) + O(\epsilon ), \ldots , \quad l=n, k, \end{aligned}$$
(62)

where \(T=\epsilon t\) is a slow time and \(\mathrm {i}=\sqrt{-1}\). We focus our attention on the oscillation amplitude behavior. The phase angle behavior can be studied in a similar way. Then, by applying the standard two time-scales perturbation method one obtains for \(A_n\), and \(A_k\):

$$\begin{aligned} 2\omega _n \frac{dA_n}{dT}= & {} - \left( d_{nn} \omega _n {A_n} + d_{nk} \omega _k {A_k} \exp (-\mathrm {i}\omega _{nk}T)\right) , \end{aligned}$$
(63)
$$\begin{aligned} 2\omega _k \frac{dA_k}{dT}= & {} - \epsilon \left( d_{kn} \omega _n {A_n} \exp (\mathrm {i}\omega _{nk}T) + d_{kk} \omega _k {A_k} \right) . \end{aligned}$$
(64)

By using (63) we express \(A_k\) in \(A_n\), giving

$$\begin{aligned} A_k=- d_{nk}^{-1} \omega _n\omega _k^{-1} \exp (\mathrm {i}\omega _{nk}T) \left( 2\frac{dA_n}{dT} + d_{nn} A_n\right) , \end{aligned}$$

and then, it follows from (64) that \(A_n\) satisfies

$$\begin{aligned} 4\frac{d^2A_n}{dT^2} + b \frac{dA_n}{dT} + c A_n=0, \end{aligned}$$
(65)

where

$$\begin{aligned} b= & {} 2(d_{nn} + d_{kk} +2 \mathrm {i}\omega _{nk}), \\ c= & {} 2\mathrm {i}d_{nn}\omega _{nk} +d_{kk}d_{nn}-d_{nk}^2. \end{aligned}$$

By using (66) we obtain

$$\begin{aligned} A_l(T)= C_{+,l} \exp (\theta _{+} T) + C_{-, l} \exp (\theta _{-} T), \quad l=n, k, \end{aligned}$$

where

$$\begin{aligned} \theta _{\pm }= \frac{1}{8}\big (- b \pm \sqrt{D_0} \big ), \end{aligned}$$
(66)

and where \(D_0=b^2- 16c\). For a given damping model the coefficients \(d_{nn}\),\(d_{kk}\), and \(d_{nk}\) can be determined, and so, the interaction of these resonant modes can be studied in detail. The analysis performed reveals that depending on the \(d_{l_1 l_2}\) parameters, the modal interaction can lead either to an increase or a decrease in the oscillation amplitudes. The reader should note, that the imaginary parts of the expression obtained for \(\theta _{\pm }\) corresponds to a small variation of the oscillation frequency induced by damping.

7 Conclusions

In this paper, the dynamics of and the buckling load for an Euler–Bernoulli beam resting on an inhomogeneous elastic Winkler foundation are studied. An analytical, asymptotic method is proposed to study the stability of the Euler–Bernoulli beam for various types of inhomogeneities in the elastic foundation and for different types of damping. Based on the Rayleigh variation principle, beam buckling loads are determined for cases of harmonically perturbed types of inhomogeneities in the elastic foundation, and for cases of point inhomogeneities in the form of concentrated springs in the elastic foundation, and for cases with rectangular inclusions in the elastic foundation. Buckling loads are influenced by inhomogeneities in the elastic foundation, and this effect exhibits a “space resonance”: that is the magnitude of the critical load depends on the inclusion location. We can control the magnitude of the buckling load by using this “space resonance” effect, and by taking a particular number of inclusions. The investigation of the beam dynamics shows the possibility of internal resonances for particular values of the beam rigidity and longitudinal force. Such types of resonances, which are usually typical for nonlinear systems, are only possible for the beam due to its inhomogeneous foundation. For large times also a beat effect can be observed. The maximal displacement during this beating was observed for a specific relation between the beam rigidity and the magnitude of the longitudinal force. Instead of a beat caused by an external force, the beat effect in the considered system is caused by modal interactions. Also analytically it was shown that damping can give rise to “added mass” effects for beams near buckling. The analytical expressions of this “added mass” effect for different damping models (including space hysteresis types) have been obtained. This effect arises as a result of an interaction between the main mode, which is close to instability, and all the other stable modes. This interaction is induced by damping, and we discussed how it depends on the type of damping model.