Dynamics and buckling loads for a vibrating damped Euler–Bernoulli beam connected to an inhomogeneous foundation

In this paper, the dynamics and the buckling loads for an Euler–Bernoulli beam resting on an inhomogeneous elastic, Winkler foundation are studied. An analytical, asymptotic method is proposed to determine the stability of the Euler–Bernoulli beam for various types of inhomogeneities in the elastic foundation taking into account different types of damping models. Based on the Rayleigh variation principle, beam buckling loads are computed for cases of harmonically perturbed types of inhomogeneities in the elastic foundation, 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. 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. The occurrence of so-called added mass effects near buckling instabilities under the influence of damping have been found. The analytical expressions for this “added mass” effect have been obtained for different damping models including space hysteresis types. This effect arises as a result of an interaction between the main mode, which is close to instability, and all the other stable modes of vibration.

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.

Statement of the problem
The equation describing the beam dynamics is given by: where u(x, t) is the beam displacement, t ≥ 0 is the time, x ∈ [0, L] is the space coordinate in the axial direction of the beam, m = Aρ is the mass of the beam per unit length, A is the beam's cross sectional area, ρ is the beam material density, E is the Young's modulus, I is the moment of the cross-section inertia, > 0 is a small parameter, a is a longitudinal force, which can have a positive (compression), or a negative (expansion) sign, and E I 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: Ac 2 0 ρ . For simplification, the bar is omitted and the final equation then takes the form: where x ∈ [0, 1]. The following initial conditions are considered: where ||v 0x x || + ||v 1 || < ∞. Here, we use the standard notation ||v|| for the norm, || f || 2 = f, f , and f, g is the scalar product in The boundary conditions are assumed to be simply supported ones: However, the applied perturbation approach, which is used in this paper, is applicable to problems with other boundary conditions.

Eigenfunctions
Our first step is to consider the unperturbed equation (2) with = 0. The unperturbed equation (2) can be solved by using Fourier's method, i.e., by substitution of u(x, t) = Re ψ(x) exp(iωt) into (2) with = 0. Then, for ψ we obtain the following eigenfunction problem where the orthonormal eigenfunctions ψ n satisfy the boundary conditions (4). For b(x) = const the eigenfunctions are well known (see [5]). To handle the case b = const we apply perturbation theory [30,31] For positive integers n and λ n > 0 let us denote by ω n the natural frequencies, which are defined by ω 2 n = λ n , where n = 1, 2, . . .. We set formally ω −k = −ω k . For a ≤ 0 and the given boundary conditions (4) it is well-known that the eigenvalues λ n of the operator L are always positive. In this case with = 0 we have stable solutions, where the natural frequencies ω n are positive and ω 2 n = λ n . For the opposite case a > 0 an instability is possible, when λ 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 λ n < 0 occurs for some n. This can be shown, for example, by the Raleigh variation principle for eigenvalues.

Perturbed eigenfunctions: non-resonance case
For small positive δ b we can apply the well-known perturbation theory [30,31]. Under the assumptions we have the asymptotical formulas for the eigenvalues and for the perturbed eigenfunctions: where Relation (8) can be used to find an approximation for the buckling load a c (see below). Note that in the case Δ kn = O(δ b ) for certain k = n the relations (8), (9), (10) have to be modified. We consider this case in the next subsections.

On condition (7)
A sufficient condition to satisfy (7), can be found as follows. Using the relation (6), we obtain For k = n one has |(πn) 2 − (πk) 2 | ≥ 3π 2 and (πn) 2 + (πk) 2 > 4π 2 . Therefore, if, for example, 4π 2 e 0 > a then condition (7) is satisfied. The resonance condition reads for some k = n. It is satisfies if for some positive integers n and k, n = k. Note that we are dealing here with an internal resonance. The number of modes involved in the resonance depends on ae −1 0 π −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.

Perturbed eigenfunctions in case of an internal resonance
Let us represent the eigenfunctions ψ n and ψ k as (see [31]) where c 1n , c 1k , c 2n , c 2k are unknown coefficients and the correction functions φ n and φ k are orthogonal in k . Then, the unknown coefficients and corrections to the eigenvalues can be found from the following 2 × 2 linear algebraic eigenvalue problem: where where tr stands for the transposed. The first eigenvalueλ + of (13) gives us the perturbation of λ (0) n , and the second oneλ − is the perturbation of λ (0) k : Theλ ± are given byλ The two independent eigenvectors (C (l) 2 ) tr , l = 1, 2 are related to the coefficients c ln in the following way: Note that the eigenfunctions are orthonormal, therefore the coefficients c ln are also orthonormal. It is convenient to write down c ln as where ξ nk is a resonance angle. Under conditionR 12 = 0 for that angle one has the following relation: IfR 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(δ 2 b ). The parameter ξ nk plays an important role, as will be shown in the coming subsection.

A new resonance effect: harmonics disappear for a long time
We suppose that damping is absent, i.e., = 0. Suppose that the initial data are given by a simple harmonic function, for example where A is an amplitude (for a more general case, where u t (x, 0) = B sin(πnx), our results are similar). We have where 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 ω n . The effect of the elastic foundation produces small contributions of order O(δ b ) in all other terms in (20) for l = n. In the resonance case [that is , when (12) is satisfied for ae −1 0 = 5π 2 ] we obtain that the two modes ψ n and ψ k give rise to significant contribution in the sum of (20), that is, and From the previous subsection it follows that C nn = c 1n = sin ξ nk , C nk = c 1k = − cos ξ nk and, as a result, we have Note that in this resonance case ω n − ω k = O(δ b ). By introducing the mean frequency and the deviation bȳ respectively, we can transform Eq. (23) as follows: Using (18) one obtains u n (t) = A cos(ωt) sin(ωt) 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(δ b ) for large times and we observe beats (see the plot in Fig. 1). The maximal beat effect, when We will use this condition (26) to have maximal beat effect in Sect. 4.

Strong Winkler foundation
In the case where the perturbation b 1 is not small, we can apply the Rayleigh variation principle: the minimal eigenvalue λ can be found as the minimum of the Rayleigh quotient Q over all test functions, that is, where the test functions ψ ∈ C 2 [0, 1] should satisfy the boundary conditions (4), and V is given by 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.

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(γ x), where γ > 0. Then, (8) gives for γ = 2πn. For γ = 2πn +γ , whereγ = O(δ b ), we can simplify the expression in the right-hand side of (28) by using 1 − cos(γ ) For the buckling load a c (n) corresponding to the stability loss for the n-th mode we approximately obtain where Δλ n is defined by (28) for γ = 2πn or by (29) 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 where δ stands for the Dirac delta function and β j are coefficients. For this case we will also observe a "space resonance" effect. Suppose that x j = ( j − 1/2)r +x, where j = 1, . . . , n d and r = 1/n d . The parameter β j may have negative, or positive values. So, we assume that all point inclusions of the elastic foundation are equidistant and the parameterx ∈ (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). Whenx = 0 all inclusions are located symmetrically with respect to the beam edges. Moreover, we suppose that n d O(δ −1 b ) . Let β j = β 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 λ n depends on the number of inclusions only and does not depend on the mode number n and the shiftx (see "Appendix"): 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") The last relation shows that a "space resonance" effect occurs, which depends onx. This effect is stronger if x = 0 (see Fig. 2). Then, for the odd k we have Δλ n = 2βδ 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 Δλ n = 0. In both cases we can use (30) for a c (n). So, we obtain a higher value for Δλ n , if n/n d is an odd integer andx = 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 δ 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 (δ b = 0.1), the star curve shows the case of the foundation without inclusions, and the solid one describes the case for a negative δ b (δ 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 δ b . For negative δ b we obtain more jumps. Note that, as it follows from the Central Limit Theorem that for random x j and n d 1, when we are dealing with a random perturbation, the effect of the inclusions is almost absent and n ) 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 γ = 2πl + O(δ b ), where l is a positive integer (this integer can be equal to n or k).
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. 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: where x 0 > d and x 0 + d < 1. Let us denote by Δb = b 0 − b min the "depth" of the inclusion. The quantity 2d is its width. If d 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 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 → 0 we obtain a point inclusion as considered in [29]. For d 1 the unperturbed non-localized eigenfunctions ψ n have the form √ 2 sin(πnx). Then, we have the asymptotics, which are valid up to corrections of the order d 2 : where z n = πn, and V [ψ n ] is defined in (27). If a buckling mode has number n, then a c,nloc (n) ≈ e 0 (πn) 2 + (πn) −2 b 0 − 2dΔb sin 2 (πnx 0 ) .
For small e 0 we can consider z n as a real valued parameter (since the minimum of Φ(z n ) is obtained for large n). Then we minimize Φ(z) with respect to z, which implies: The value of a c can be found by the condition that min V = 0 , which implies: where the critical mode number is and where [z] stands for the integer nearest to z.
Note that under our assumptions Δ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, where θ(y) is a well localized function at y = 0 (for instance, we can take θ(y) = exp(−y γ ), γ ≥ 2). Here, the parameter r defines the characteristic radius of the test mode localization. For small r d one has: for r → 0, where the positive constants c 1 , c 2 , and c 0 are independent of r , and are given by The same expression (36) can be obtained for an infinite beam. Thus where o(1) → 0 as e 0 → 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 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.

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 ψ n : where X n (t) are unknown functions. Then, the coefficients X n (t) have to satisfy: 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 ψ n * = Ψ , which is close to instability, and n * denotes its number. Note that it is possible that n * = 1, as was shown in the previous section. For small |a − a c | the frequency ω n * is small and close to zero for a certain n * . For that reason we introduce a new small parameter ω n * = μ > 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 δ b , and μ are small parameters which are related in the following way: where μ 0 andδ b are positive parameters, which are independent of the small, positive parameter . Further, to solve the equations describing the time evolution of the stable (fast) modes X n with n = n * , we use a multiple time-scale perturbation method with fast time t and slow time T = t . Note that the parameter δ b is not involved directly in the suggested multiple time-scale approach, but it is present in the equation for a c (δ b ) and in the orthonormal eigenfunctions ψ n (x), which depend on δ 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 = n * . We obtain then the following system of equations for Y = X n * and the fast modes X n (t) with n = n * : where B is the following differential operator: The + sign in (41) before the term 2 μ 2 0 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 2 μ 2 0 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 and for internal damping where L s is a linear operator, L s u = ∂ 4 u ∂ 4 ξ . The simple damping (air damping) and Kelvin-Voigt damping (referred to as the SD and KV cases) occur if h td = δ(t −τ )δ(x −ξ) 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 − ξ |, for example, it may be a Gaussian function: where σ a positive parameter with σ 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: whereX 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 termsX n (t, T ) produce small contributions in Y with respect to the termsX n (T ). To see this, let us consider Eq. (41) with X n =X n (t, T ) in the right-hand side. The corresponding asymptotic solutionỸ (t, T ) involves the harmonics sin(ω n t), cos(ω n t) with n = n * and coefficients depending on T . This part of the solution has order O( ) because the right-hand side of Eq. (41) is proportional to . Now let us consider (41) with X n =X n (T ) in the right-hand side, which is proportional to 2 and depends only on T . The corresponding solutionȲ is defined by and, in general, has order O(1). Therefore, we conclude that X n =X n (T ) andȲ (T ) satisfy Eqs. (41) and (42) up to terms of higher orders in . In fact, let us substitute Y =Ȳ into Eq. (42), giving The function X n (t, T ) =X n (T ) satisfies the last equation up to and sinceȲ is a function of T only, it follows that the functionX n (T ) is of order . 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 thatX where the constants α n determine an interaction force between the unstable mode and the n-th stable one, and have the form: and for SD, KV damping, internal space hysteresis damping, and for external space hysteresis damping, respectively. Taking into account that in (47) which follows fromX n (T ) = α n dzȲ dT , and from the arguments given after (47), and so we finally obtain:: where the relative "added mass" m add is where and 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 0 , which determines the deviation between a and the critical value a c . In fact, this effect is proportional to the coefficient α n , which is defined by the relations (48), (49) (50) and (51). Now we can conclude the following: (a) in the SD case α n = 0 and the added mass effect is absent at the order 1 level; (b) in the KV damping case α n = 0 for beams on homogeneous elastic foundations, and α n is proportional to δ b for inhomogeneous elastic foundations and hinged supported beams; (c) similar results can be obtained for the external and internal SH cases. Let σ 1. Then, we conclude that for simply supported beams and external or internal space hysteresis damping α n = O(exp(−1/σ )) for beams on homogeneous elastic foundations, and α n = O(δ b ) for inhomogeneous elastic foundations.
Note that the modal interaction does not change (up to order ) the value of the critical Euler force a c . Indeed, the value of a c depends on parameters involved in the self-adjoint operator 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 where , 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 (κ 4μ 2 0 (1+m add )). When a < a c it follows from (58) there exists a weakly decreasing mode with Re θ 1 ≈ −μ 2 0 κ −1 , and for a > a c we also obtain a weakly decreasing mode with Re θ 1 ≈ −μ 2 0 κ −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 (κ 4μ 2 0 (1 + m add )). For a < a c we have Re(θ 1 ) = −κ/2(1 + m add ) and Im(θ 1 ) ≈ μ 0 √ 1/(1 + m 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.

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 ψ n , ψ k be two resonant modes (see Sect. 3.4). For these modes, we introduce the detuning parameter ω nk by Let us consider the most interesting case of space hysteresis damping: To describe the mutual evolution of two resonant modes, we are restricting ourselves to the following solution form: where X n , X k are unknown functions. We neglect here effects of added mass and buckling, and we assume that ω n , ω 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 d 2 X n dt 2 + ω 2 n X n = − d nn where 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: where T = t is a slow time and i = √ −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 : By using (63) we express A k in A n , giving and then, it follows from (64) that A n satisfies By using (66) we obtain 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 θ ± corresponds to a small variation of the oscillation frequency induced by damping.

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.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. We observe that the quantityS n is the real part of sum in the geometric sequence with initial value exp √ −12π(x + r/2) and common ratio q = exp( √ −1 2πnr), leading directly to (31) and (32).