1. Introduction

The buckling and postbuckling response of beams/struts (of either infinite or finite length) on elastic foundation has been the subject of a huge number of publications over the last 50 years, which were based on various linear as well as nonlinear subgrade models [1-3]. Even a limited reference of all these investigations lies beyond the scope of the present work.

This simple structural component can rather easily simulate the stability response related to practical problems of various disciplines. Among these, one may quote the lateral thermal buckling of railway tracks [4], the stability of the top chord of low (pony) truss bridges [5], the mechanical response of metallic aortic stents [6] etc. Other theoretical approaches have also been reported [7-9], which however did not consider any underlying physical problem.

The overwhelming majority of the works cited above (and of numerous others not referenced here for brevity) was based on the assumption that the subgrade is always in contact with the supported structure (bilateral), i.e. that the foundation always reacts, either in tension or in compression.

However, the simple ideal beam only in axial compression resting on an elastic foundation may also model practical applications, in which the subgrade is of the so-called one-way, i.e. tensionless, of unilateral type. In such situations, a separation between parts of the beam and the foundation may occur, and hence the relevant problem is a discontinuous – unbonded contact one. Such a phenomenon has been reported in the literature, but in a limited number of works. For instance, in an earlier paper [10], Yun and Kyriakides simulated the “beam-mode buckling” exhibited by buried pipelines in compression (where a section of the line lifts through the ground), adopting a tensionless Winkler foundation. Later, Silveira et al. [11] studied – among other types of structures – the response of a simply supported beam acted upon simultaneously by end moments and a compressive axial force, resting on a Winkler tensionless foundation. The authors used a semi-analytical methodology via a Ritz type approach, which led to the determination of the contact and no contact region between beam and subgrade.

Moreover, Battiprolu et al. [12, 13] in order to simulate the response of a pinned-pinned beam interacting with polyurethane foam foundation, adopted a beam-model on a nonlinear tensionless viscoelastic subgrade under axial compression and distributed as well as concentrated transverse loading. They also investigated the effect of various parameters, such as the relative foundation-beam stiffness, nonlinearities etc. Separation between beam and subgrade was evident also in this case.

Finally, in a very recent paper [14], Ioakimidis derived the conditions of complete contact for a beam on a tensionless Winkler foundation, using symbolic computations available in modern quantifier elimination software [Mathematica, [15]]. In this study, the case of pure axial compression was considered, and as stated, the adopted subgrade model is frequently actually the case in practice.

To this end, when the above described simple beam-model lacks imperfections and lateral restrictions, if it is loaded by a pure axial compressive force, only four (4) postbuckling deformation patterns may occur, as the load gradually increases from zero. Namely, the 1^st one (abbreviated in what follows by PB1) is associated with the in-plane buckling “away” from the foundation (full separation – no contact at all), that is realized by the 1^st Euler stable buckling mode in a single semi-wave form, with well-known properties. The 2^nd pattern (PB2), out-of-plane Euler buckling in any mode, also stable, may as well occur but is not affected by the presence of the subgrade (friction neglected). The corresponding critical loads may be greater of small than the ones of the 1^st pattern, depending on the beam’s properties and cross-section orientation. The 3^rd simplest possible postbuckling configuration (PB3) is related to a single wave-form-like deformed shape in full contact with the subgrade, i.e. the whole beam is deflected inside the foundation. This will be possible only for certain combinations of the beam and foundation parameters [5, 14, 16-17] and the critical load may also be smaller or greater than the ones of PB1 and PB2. The last pattern (PB4) includes discontinuous contact, with yet unknown in number and length consecutive regions in contact and no contact. However, energy considerations and application of the principle of Mathematical Induction [18] have shown [19] that, in the foregoing case, the buckled deformed configuration is associated with only one part of the beam in contact with the foundation and the remaining not in contact. This finding was based on a Pasternak tensionless subgrade model, and hence is a fortiori valid for the Winkler model as well.

Focusing on establishing the nonlinear stability of PB4, which is an unstable case with profound symmetry, the present work tackles the problem by constructing the total potential energy function of the system and exploring its nature via Catastrophe Theory [20-22]. The three parameter Singularity arising was found to be a special reduced case of the Butterfly, namely the symmetric one [23], with its Bifurcation Set strongly resembling Swallowtail-like sections of the original Singularity [24]. This set, evaluated by means of successive eliminations, as demonstrated by Deng [25], produced critical buckling loads and corresponding deformed configurations. The results obtained were compared with the ones of two more alternative approaches and good agreement was established. It was found that the unstable configuration dealt with herein may appear either prior to or between the in- and-out-of-plane Euler buckling, which may lead to sudden jumps and unexpected bifurcations, a fact of importance for practical design and simulations.

This could be avoided by carefully selecting the beam and subgrade parameters, which is not an easy task, since especially the stiffness of the foundation cannot be directly determined [26]. Similar undesirable phenomena have also been reported in the literature, but for bilateral contact only [7, 8, 27, 28].

For the completeness and full validation of the present theoretical findings, in their entirety, experiments should be carried out. This task is ongoing by the authors, and hopefully its outcome will be soon available.

2. Problem Description

We consider a simply supported inextensional Euler-Bernoulli beam of flexural rigidity EI and of uniform cross section, resting on a linear (Winkler type) tensionless foundation, characterized by a stiffness parameter k_f, and acted upon by an axial compressive force P. Let W(y) be the beam’s in-plane flexural deflection. The material of the beam is assumed linearly elastic and there are no lateral restrictions. The undeformed pre-buckling configuration is shown in Figure 1, that follows.

Figure 1. Axially compressed simply supported beam resting on a Winkler type unilateral elastic foundation

We then focus on the 4^th possible postbuckling behavior described in the Introduction, which is related to one portion of the beam in contact with the foundation and the remaining one with no contact. Due to the nature of the foundation (tensionless and linear) the problem is symmetric and thus one may arbitrarily choose to investigate the case of either the 1^st or the 2^nd region being in contact. In this work we consider the 1^st one in contact, depicted in Figure 2.

Figure 2. Postbuckling configuration of the beam dealt with

Hence, it is convenient to express this shape in a Fourier series expansion, retaining only the two leading coefficients, which suffice for the problem at hand, i.e.

(1)

If more than two terms are used, no additional accuracy is expected, only obstacles in the analysis that follows.

This smooth function satisfies the boundary conditions, and moreover, in order that there is free rotation at the two ends of the beam and one global extremum in the deflection of each region, it is easily found that the ratio of the two state variables q₁ and q₂, designated as β = q₁/q₂, must satisfy the inequality . Additionally, for the function given in Eq. (1) to possess only one more root within the interval , symbolic computations using the powerful Reduce command embedded in Mathematica [15] yield the following conditions, which represent a logical expression of a clear dissolution – conjunction type,

(2)

resulting to the evaluation of the length a depicted in Figure 2, which is found equal to

(3)

3. Mathematical Formulation

We then proceed with the evaluation of the system’s total potential energy function V_T being in fact the sum of the strain energy due to bending U₁, the strain energy due to subgrade reaction U₂, and the work of the external force Ω. In doing this, we follow standard procedures of Mechanics and we adopt an approximate nonlinear curvature expression [29, 30] of the form

(4)

The various components of V_T are thereafter formulated as:

(5)

The last expression is the so-called fixed length condition, adopted from Gilmore, page 258 [21].

We then introduce the following dimensionless parameters:

(6)

Moreover, we approximate function arcsec[p₁/p₂], which appears in the evaluated expression of U, according to the formula

(7)

up to the 6^th order. Higher order approximants lead to identical results, since z, i.e. -2p₁/p₂ is (in absolute value) restricted to vary between the boundaries of inequality (2).

Radicals appearing in the various expressions of the components of the total potential energy function are also approximated via series expansions up to the 6^th order. The reason for the choice of this specific order of truncation, and not of a higher one, is explained at the end of this Section.

The final outcome of the dimensionless total potential energy function of the beam yields:

(8)

where

(9)

This function is obviously smooth (with continuous derivatives up to the 6^th order), since it involves only simple algebraic manipulations of the parameters and no discontinuities, due to the absence of denominators that may be equal to zero. Hence, the potential and the families of functions that it represents may be treated via Catastrophe Theory.

Evidently, the strongly nonlinear total potential energy function of the original (continuous) structural system depends on three parameters: the magnitude of the external compressive load λ, the stiffness of the foundation μ and the length of the beam’s portion in contact with the subgrade a. The truncated potential also depends on three parameters, as in (9), which in turn depend on the original ones (β silently refers to a). Since the proposed analysis leads to a potential dependent only on one generalized coordinate, and the problem at hand is an elastic stability one, the 6^th order truncation appears a reasonable choice. Note that the swallowtail and butterfly singularities associated with such a choice are scarce in structural engineering; additionally, the results produced will be verified in a satisfactory manner, through the application of two more different approaches, as presented in Section 4.

4. Catastrophe Theory Analysis

In this Section we explore the nature of the truncated total potential energy function of the system, according to the Theory of Catastrophes. Apparently, the function given in Eq. (8), although it can be seen as a reduced – special case of the Butterfly Singularity [21, 22], it still represents a symmetric butterfly function, exhibiting in general a symmetric catastrophe manifold [23]. This feature will be sought by exploring the nature of the foregoing Singularity via an alternative scheme. In doing this, we employ an elimination method, reported by Deng [25], and using simple Algebra, we reach to the Bifurcation Set (B_s) of the Catastrophe at hand. Dropping the subscript of p₁ in Eq. (8), the corresponding equilibrium equations at the critical states and the final expression of the B_s are given by:

(10a,b)

Using the expressions of ξ_i, i=1,2,3 given in (9), and after cumbersome symbolic manipulations, Eq. (10b) yields a second order polynomial equation with respect to λ, the coefficients of which are lengthy nonlinear functions of β and μ, not presented herein for brevity. However, for arbitrary positive values of μ and for β varying according to inequality (2), it is found that this equation has always two distinct real positive roots, namely λ₁ and λ₂. Their values are given in Eqs. (11a,b) as follows:

(11a)

(11b)

The parameterized families of surfaces defined by Eq. (10b) could be represented as three-dimensional contour plots in the (λ, β, μ) space and then compared with standard bifurcation sets of the Elementary Catastrophes in a geometrical manner [24]. This will be demonstrated in the next Section.

For real applications, i.e. for a specific value of one may engage expressions (11a,b) in a minimization procedure via Mathematica [15], and evaluate the critical (buckling) load and the corresponding critical value of β = β_cr. The latter is the quotient of the state variables at the instance of buckling.

5. Results and Discussion

Initially the Bifurcation Set defined in Eq. (10b) is presented graphically in terms of three-dimensional contour plots. Due to its strong nonlinearity, a full range reconstruction requires immense computer power, and hence partial range reconstructions are produced for increasing values of λ and μ. These and shown in Figures 3 and 4.

Figure 3. Partial 3D contour plots of the Bifurcation Set for small values of λ and μ, seen from different viewpoints

Figure 4. As in Figure 3, but for large values of λ and μ

The surfaces contained in these Figures clearly resemble to Swallowtail-like sections of the Symmetric Butterfly Catastrophe, as given in [24]; this is rather expected, due to the nature of the truncated potential. These surfaces however cannot be produced directly from this potential using advanced symbolic manipulations and sophisticated graphics (due to unavoidable memory overflow), but only via the Bifurcation Set evaluated according to the above scheme. Moreover, and according to the preceding analysis, numerical results are presented in what follows, considering seven (7) combinations (Combos) of beam and subgrade, adopted from the relevant literature.

More specifically, Combo 1 is taken from the work of Kounadis et al. [29], while the remaining ones are those used in [11] and [31]. The critical (buckling) load λ_cr, found in this investigation, is compared with the corresponding buckling loads of in plane buckling (1^st mode – no contact), of out of plane lateral buckling (1^st mode – no influence of subgrade reaction) and of the case of full contact in a simple semi-wave, if existent. The outcome of the analysis is shown in Table 1, from where one may observe that – for the combinations considered – postbuckling response PB3 (full contact) is not possible.

Table 1. Numerical results for the beam – subgrade combinations considered, based on the proposed Catastrophe Theory approach

Before discussing these results, it is considered appropriate to compare them with ones obtained from different methods. In doing this, two additional approaches are utilized for the problem at hand. The 1^st is based on a flat Galerkin scheme and the 2^nd on modeling the system in SAP 2000 software [32].

The former approach is the one described in [19] properly adjusted to the Winkler subgrade model, by vanishing the parameter associated with the shear layer in the Pasternak model. The latter takes advantage of the modeling capabilities of the aforementioned software according to the following steps: (a) for each Combo number the beam is designed as a single straight elastic frame element, simply supported inextensional and restricted to deform only in-plane; in cases where the libraries of the software do not include the cross-section corresponding to the specific Combo, it is drawn externally and imported afterwards as a dxf file, (b) the beam is discretized – divided densely into frame sub-assemblages with free to rotate internal joints in all directions, but also restricted to translate only vertically, (c) at each intermediate joint a gap – compression only – nonlinear link is inserted, with its stiffness defined from each individual Combo accordingly, and (d) a nonlinear buckling analysis is performed. The self-weight of the beam is not accounted for.

Both alternative approaches yield results in good agreement with the ones obtained from the method proposed herein. The 1^st (Galerkin) approach produces upper-bounds, while the SAP model lower-bounds. Interestingly enough, the FE model failed to converge in more than two buckling modes; this may be considered as a simple – but existent – indication of absence of a 3^rd mode, according to the findings of the Mathematical Induction mentioned in the Introduction, also for the simplest subgrade model. The results of the two approaches are summarized in Table 2.

Table 2. Numerical results for the beam – subgrade combinations considered, based on the proposed Catastrophe Theory approach

From all the above it is readily perceivable that, in some cases, the critical load corresponding to the unbonded contact configuration is smaller than both the in-plane and out-of-plane classical Euler load. This implies that an unstable bifurcation may occur prior to stable states, a finding important for design and simulations. In other cases, this specific load lies between the Euler loads. As the external compressive force increases from zero, this last finding may be interpreted as a possibility of sudden jumps between unstable and stable states. This is a well-known property of the butterfly singularity, and may occur even in the case of the simplest subgrade model (dealt with in this work). This undesirable behavior should be accounted for, when addressing the buckling of beams on elastic foundations. Similar jumping phenomena have also been reported in the literature [7, 8, 27, 28] for nonlinear bilateral subgrade models. As mentioned in the Introduction, to avoid such a behavior, a proper choice of beam and subgrade parameters should be made.

Characteristic quantitative plots of the postbuckling deformed shape of the beam for PB4 are presented in Figure 5, for different values of μ. As this parameter increases, the corresponding length of the part of the beam in contact decreases, as well expected from the whole analysis.

Figure 5. Qualitative buckled configurations for (a) large, (b) medium, and (c) small values of μ

Figure 6. Sections – contour plots of the Bifurcation Set for four (4) of the combinations of Table 1

Finally, it is also expedient to show contour (β, λ) plots – sections of the Bifurcation Set for some of the above combinations of the parameters, as illustrated in Figure 6. From these plots one may see the numbers and types of the stationary points in the chosen territories of the parameter space, as well as the evaluated critical value of β (as given in the above Table), which provides the smallest critical load. The reconstruction of these kind of plots is a standard procedure accompanying Singularity Theory, and this is the main reason of their presentation.

To the knowledge of the authors, no experimental works have been yet reported, concerning the problem at hand. It is recommended that such experiments should be carried out, in order to fully validate the theoretical findings presented so far.

6. Conclusions

The most important conclusions drawn from the present study are the following:

(a) The consideration of discontinuous (unbonded - unilateral) contact between beam and foundation, being an unstable postbuckling response, leads to a total potential energy function which is a special reduced form of the Butterfly singularity, namely the symmetric one. The resulting three-parameter bifurcation set, an outcome of successive eliminations, resembles Swallowtail-like sections of this particular Catastrophe reported in the literature. Results obtained from the present approach, in terms of critical loads and unbonded beam lengths were found in satisfactory agreement with ones of two alternative approaches.

(b) For the cases considered, unstable bifurcations prior to Euler buckling may occur, while also sudden jumps may be possible (a salient feature of the Butterfly Catastrophe), while full contact in a single semi-wave form is not possible. For structural design, these findings are considered of importance, and should be accounted for both in practice and simulations. Similar phenomena have also been reported in the literature for nonlinear bilateral subgrade models.