Introduction

The concept of beams on elastic soil foundation has been widely used in different fields of engineering, such as strip foundation, railroads tracks, building, dams and airport runway. The soil mechanic exhibits a very complex behavior of foundations due to the heterogeneity, physical composition, and presence of imperfections and pores of soils. Concepts, status and various analysis methods of the soil-structure interaction researches have been illustrated in the recent art state review (Prakash et al. 2016).

The quantification of the soil-structure interaction is a challenge until the present moment. The modeling of the contact between the structure and the soil foundation is a primordial task of this analysis. Analytical solutions are restricted compared to numerical studies using two-parameter soil foundation model of beams resting on isotropic or anisotropic elastic half-space (Johnson 1985; Kachanov et al. 2003). In numerical domain, various models have been largely developed for modeling the soil foundation, which are classified into three categories: (1) continuum models, (2) mixed models and (3) spring models.

In continuum mechanic concept, the medium is defined by a continuously distributed matter through the half-space that the constitutive law can be described by a linear elastic isotropic behavior (Irgens 1980). The solution for a simplified continuum using the finite element idealization was developed by Reissner (1967). The continuum can be analyzed with many numerical methods, such as the finite element method (FEM), the boundary element method (BEM) or combined methods between FEM and BEM, which are suitable to the soil-structure interaction analysis. Particularly, FEM is well known and widely used in many approaches to study soil-structure interaction behavior and BEM shows many advantages in the modeling field showing a high accord with infinite and semi-infinite spaces (Bolteus 1984; Tezzon et al. 2015).

Due to the complexity of the interaction problems, analytical solutions are rarely used and desired alternatives would be a numerical approach (Dinev 2012; Hassan and Doha 2015; Ai and Cai 2016). The FEM is still popular method in this domain (Bourgeois et al. 2012; Su and Li 2013) but it has disadvantages compared to BEM and spectral element method (SEM) (Omolofe 2013; Mokhatari et al. 2016) using absorbent frontier modeling. The FEM requires the discretization of the domain to high number of finite elements but the problem can be solved with BEM where only the boundary of the domains involved can be discretized (Padron et al. 2011; Ai and Cheng 2013; Ribeiro and Paiva 2014).

Beam using Winkler foundation model (1867) is used in various practical problems. In this approach, the foundation flexibility is considered as a set of continuous springs and has employed to model soil-structure interaction (Kim and Yang 2010; Chore et al. 2010; Sapountzakis and Kanpitsis 2011; Raychowdhury 2011; Limkatanyu et al. 2012). Springs introduced provide a resistance in the vertical direction only confining the deformation of the soil foundation. Evidently, the one-parameter soil foundation model suffers a handicap due to the discontinuity in the supporting medium. To improve the Winkler model, two-parameter and three-parameter soil foundation models have been developed taking into consideration shear deformations of the soil. To overcome the one-parameter foundation model deficiencies, many researchers proposed various foundation models to describe rigorously soil response insuring the interconnection between vertical springs.

In the first approach, the interaction between springs has been established by an elastic tensioned membrane (Filonenko 1968) and the second one uses beam elements or a plate to interact between them (Hetényi 1966). In this case, the tensioned membrane is quantified by a shear parameter. Third, to model the mutually effect between springs, Kerr (1964) integrated a shear layer dividing the soil medium of foundation to two different layers.

The two-parameter soil foundation has shown considerable developments during the last decade. A new finite element formulation was developed to study the shallow and the raft foundation response eliminating the limits of Winkler model (Mullapudi 2010). In the way to study large deflections of functionally graded beam resting on two-parameter elastic foundation, a finite element procedure was developed (Gan and Kien 2014). Finally, the effect of material non-homogeneity and the two-parameter elastic foundation were used to quantify the response of simply supported beams. The foundation medium behavior is assumed to be linear elastic, homogeneous and isotropic with two parameters describing the reaction of the elastic foundation on the beam (Avcar 2016).

In this analysis, the soil-structure interaction problem has been studied using effectively two-parameter model of the layered soil. Both the beam and the substrate are described by means of FEM integrating shear deformations. However, to ensure vanishing displacement at the frontier of the substrate, mesh has to be extended far away from the loaded region. To improve the computational efficiency, two-dimensional finite elements are refined in the loaded area. The modeling uses plane stress state for the shear beam in adhesive contact with plane strain state of the soil foundation. Adding, a parametric study is elaborated to show the influence of (1) the horizontal behavior of the soil foundation, (2) mechanical properties of the soil foundation and (3) ballasted layer.

Physical modeling

Transverse deflections of a thin elastic beam can be produced due to the reciprocal effect of the loading and the foundation reactions. Really, the reaction forces of layered soils are proportional to the deflection at each point at the interface medium. Among deficiencies of the Winkler model is that the displacement discontinuity which appears between the loaded and unloaded regions of the foundation surface. The displacement field must be a continuous function along the all surface (Fig. 1).

Fig. 1
figure 1

a Practical soil foundation and b Winkler’s foundation

In the Winkler model, vertical soil sections of the foundation are modeled by identical, discrete elastic springs characterized by their mechanical property, Ks (Fig. 1b). To improve prediction of the entire system behavior, the contribution of the interface medium between the foundation and the beam must be integrated. The literature shows many methodologies interpreting the soil-beam interaction response. This interaction was established by means of flexural element (beam or plate) (Fig. 2), shear layers (Fig. 3) or pre-tensioned membranes (Fig. 4).

Fig. 2
figure 2

Interaction using beam or plate

Fig. 3
figure 3

Shear layer model

Fig. 4
figure 4

Tensioned membrane model

On the other hand, to take into account the lateral interaction due to the soil cohesion, a plate is used introducing the shear effect (\(K_{\text{G}}^{{}} \left( {\frac{{{\text{d}}_{{}}^{2} w(x)}}{{{\text{d}}x_{{}}^{2} }}} \right)\)) (Fig. 5). Further, a set of springs are introduced in the longitudinal direction of the soil foundation guaranteeing the interaction between longitudinal and transverse behaviors (Fig. 6).

Fig. 5
figure 5

Shear effect model

Fig. 6
figure 6

Shear spring elements

Finite element formulation

To improve the soil-structure interaction behavior, the incorporation of the shear phenomenon seems so necessary. The change of reactions at the interface between the beam and the soil can be expressed as a combination of:

$$R(x) = K_{\text{s}}^{{}} w(x) - K_{\text{G}}^{{}} \frac{{{\text{d}}_{{}}^{2} w(x)}}{{{\text{d}}x_{{}}^{2} }},$$
(1)

where w(x) is the vertical deflection of centroidal line, Ks and KG are the sub-grade and the shear foundation modulus, respectively.

The theory of beams leads to the fourth order of the differential equation for the deformed central line of the beam resting on elastic soil:

$$E_{\text{b}}^{{}} I_{\text{b}}^{{}} \frac{{{\text{d}}_{{}}^{4} w(x)}}{{{\text{d}}x_{{}}^{4} }} + K_{\text{s}}^{{}} bw(x) - K_{\text{G}}^{{}} b\frac{{{\text{d}}_{{}}^{2} w(x)}}{{{\text{d}}x_{{}}^{2} }} - q(x) = 0,$$
(2)

where \(E_{\text{b}}^{{}} I_{\text{b}}^{{}}\) is the flexional rigidity of the beam, q(x) is the external loading on the structure and b is the width of the beam.

The finite element formulation described is the same used in various developments (Cen et al. 2009; Logan 2012; Zienkiewicz et al. 2005). The beam is divided into a number of elements of the length L and flexural rigidity \(E_{\text{b}}^{{}} I_{\text{b}}^{{}}\). The interpolation of the displacement field, \(w_{\text{e}}^{{}} (x)\), for a linear finite element can be written as:

$$w_{\text{e}}^{{}} (x) = \left\langle {N(x)} \right\rangle \left\{ {q_{\text{e}}^{{}} } \right\}.$$
(3)

The shape functions are \(N_{1}^{{}} (x) = 1 - 3\left( {\frac{x}{L}} \right)_{{}}^{2} + 2\left( {\frac{x}{L}} \right)_{{}}^{3} ,\) \(N_{2}^{{}} (x) = x - 2\frac{{x_{{}}^{2} }}{L} + \frac{{x_{{}}^{3} }}{{L_{{}}^{2} }},\) \(N_{3}^{{}} (x) = 3\left( {\frac{x}{L}} \right)_{{}}^{2} - 2\left( {\frac{x}{L}} \right)_{{}}^{3}\) and \(N_{4}^{{}} (x) = - \frac{{x_{{}}^{2} }}{L} + \frac{{x_{{}}^{3} }}{{L_{{}}^{2} }}\), and \(\left\{ {q_{\text{e}}^{{}} } \right\}\) is the degree of freedom vector. The strain energy for the beam element, U, is combined of:

$$U = U_{\text{B}}^{{}} + U_{\text{F}}^{{}} ,$$
(4)

\(U_{\text{B}}^{{}}\) and \(U_{\text{F}}^{{}}\) are the strain energies for the bending beam and the soil foundation.

In this case, the strain energy of the soil foundation can be divided in two components as:

$$U_{\text{F}}^{{}} = U_{\text{s}}^{{}} + U_{\text{G}}^{{}} .$$
(5)

For the two-parameter foundation model, Shen (2012) gave the strain energy expression, UF, as:

$$U_{\text{F}}^{{}} = \frac{1}{2}\int\limits_{L}^{{}} {K_{\text{s}}^{{}} } w_{{}}^{2} (x){\text{d}}x + \frac{1}{2}\int\limits_{L}^{{}} {K_{\text{G}}^{{}} } \left( {\frac{\partial w(x)}{\partial x}} \right)_{{}}^{2} x.$$
(6)

Substituting the relation (3) into the Eq. (6), we obtain:

$$U_{\text{F}}^{{}} = \frac{1}{2}\left\{ {q_{\text{e}}^{{}} } \right\}_{{}}^{t} \left[ {\int\limits_{L}^{{}} {K_{\text{s}}^{{}} \left[ {N(x)} \right]_{{}}^{t} \left[ {N(x)} \right]} dx} \right]\left\{ {q_{\text{e}}^{{}} } \right\} + \frac{1}{2}\left\{ {q_{\text{e}}^{{}} } \right\}_{{}}^{t} \left[ {\int\limits_{L}^{{}} {K_{\text{G}}^{{}} } \left[ {\frac{{{\text{d}}N(x)}}{{{\text{d}}x}}} \right]_{{}}^{t} \,\left[ {\frac{{{\text{d}}N(x)}}{{{\text{d}}x}}} \right]{\text{d}}x} \right]\left\{ {q_{\text{e}}^{{}} } \right\}.$$
(7)

Or,

$$U_{\text{F}}^{{}} = \frac{1}{2}\left\{ {q_{\text{e}}^{{}} } \right\}_{{}}^{t} \left[ {MK_{\text{s}}^{{}} } \right]\left\{ {q_{\text{e}}^{{}} } \right\} + \frac{1}{2}\left\{ {q_{\text{e}}^{{}} } \right\}_{{}}^{t} \left[ {MK_{\text{G}}^{{}} } \right]\left\{ {q_{\text{e}}^{{}} } \right\}.$$
(8)

The corresponding stiffness matrices are

$$\left[ {KM_{\text{s}}^{{}} } \right] = \int\limits_{L}^{{}} {K_{s}^{{}} \left[ {N(x)} \right]_{{}}^{t} \left[ {N(x)} \right]} dx = \frac{{K_{1}^{{}} }}{420}\left[ {\begin{array}{*{20}c} {156} & {22L} & {54} & { - 13L_{{}}^{2} } \\ {} & {4L_{{}}^{2} } & {13L} & { - 13L_{{}}^{2} } \\ {} & {} & {156} & { - 22L_{{}}^{2} } \\ {\text{Sym}} & {} & {} & {4L_{{}}^{2} } \\ \end{array} } \right]$$
(9)
$$\left[ {KM_{\text{G}}^{{}} } \right] = \int\limits_{L}^{{}} {K_{\text{G}}^{{}} } \left[ {\frac{{{\text{d}}N(x)}}{{{\text{d}}x}}} \right]_{{}}^{t} \left[ {\frac{{{\text{d}}N(x)}}{{{\text{d}}x}}} \right]{\text{d}}x = \frac{{K_{2}^{{}} }}{30}\left[ {\begin{array}{*{20}c} {36} & {3L} & { - 36} & {3L} \\ {} & {4L_{{}}^{2} } & { - 3L} & { - L_{{}}^{2} } \\ {} & {} & {36} & { - 3L_{{}}^{2} } \\ {\text{Sym}} & {} & {} & {4L_{{}}^{2} } \\ \end{array} } \right],$$
(10)

with \(K_{1}^{{}} = K_{\text{s}}^{{}} bL\) and \(K_{2}^{{}} = \frac{{K_{\text{G}}^{{}} b}}{L}\)

In the same manner, the strain energy of the beam bending is

$$U_{\text{B}}^{{}} = \frac{1}{2}E_{\text{b}}^{{}} I_{\text{b}}^{{}} \left\{ {q_{\text{e}}^{{}} } \right\}_{{}}^{t} \left[ {\int\limits_{L}^{{}} {\left[ {\frac{{{\text{d}}_{{}}^{2} N(x)}}{{{\text{d}}x_{{}}^{2} }}} \right]_{{}}^{t} } \left[ {\frac{{{\text{d}}_{{}}^{2} N(x)}}{{{\text{d}}x_{{}}^{2} }}} \right]dx} \right]\left\{ {q_{\text{e}}^{{}} } \right\}.$$
(11)

Substituting the relation (3) into the Eq. (11), the stiffness matrix can be deduced as:

$$\left[ {K_{\text{b}}^{{}} } \right] = \frac{{E_{\text{b}}^{{}} I_{\text{b}}^{{}} }}{{L_{{}}^{3} }}\left[ {\begin{array}{*{20}c} {12} & {6L} & { - 12} & {6L} \\ {} & {4L_{{}}^{2} } & { - 6L} & {2L_{{}}^{2} } \\ {} & {} & {12} & { - 6L} \\ {\text{Sym}} & {} & {} & {4L_{{}}^{2} } \\ \end{array} } \right].$$
(12)

The loading vector \(\left\{ {F_{\text{e}}^{{}} } \right\}\) depends on the distributed forces along the element side. For a uniform load \(q(x) = q_{0}^{{}}\), the corresponding loading vector can be formulated as:

$$\left\{ {F_{\text{e}}^{{}} } \right\} = \int\limits_{L} {\left[ {N(x)} \right]_{{}}^{t} } q(x){\text{d}}x = q_{0}^{{}} \int\limits_{L} {\left[ {N(x)} \right]_{{}}^{t} } {\text{d}}x.$$
(13)

Numerical study

Validation of the program

First, the beam on Winkler foundation model is used to show the ability of the numerical program conceived for this concern. This beam is studied using mesh-free method (Binesh 2012). The reaction modulus of the soil foundation is of 104 kN/m3 and the concentrated load of 100 kN applied at 2 m from the left side of the beam (Fig. 7). The beam is of length 14 m and flexural rigidity is assumed of 2604.167 kN/m2, which is equivalent to a beam with thickness of 0.25 m and a Young’s modulus of 2 × 103 MPa under plane stress conditions. The Poisson ratio is assumed of 0.25.

Fig. 7
figure 7

Beam on elastic soil foundation

The model consists of 73 nodes and a background mesh with 14 blocks for numerical integration. Obtained results using mesh-free method and exact solution are shown in Fig. 8a. For comparison, there is well agreement between the exact solution, the mesh-free method and this approach, which confirms the accuracy and efficiency of the proposed numerical solution.

Fig. 8
figure 8

Vertical deflection along the beam. a Binesh (2012). b The present approach

Figure 9 depicted horizontal displacements of interface continuum nodes. The figure shows important displacements in the left region of the applied force but a slight variation accompanied by weak displacements at far nodes.

Fig. 9
figure 9

Horizontal displacement along the beam nodes

Figure 10 shows the horizontal interaction model. In this case, springs are introduced in vertical and horizontal directions, simultaneously.

Fig. 10
figure 10

Horizontal and vertical spring model

Horizontal springs can be justified to avoid the handicap raised in the Winkler model. When the lateral interaction is considered in the analysis, horizontal displacements are different from those obtained by Winkler approach. In this case, an improvement is highly observed varying from 45.13% at x = 0 m to 100% for \(x \ge 5\;{\text{m}}\) along the soil length (Fig. 11). The lateral interaction appears with no effect on vertical displacements (Fig. 12).

Fig. 11
figure 11

Horizontal displacement field

Fig. 12
figure 12

Vertical displacement field

Really, the calibration of horizontal stiffness plays a primordial contribution in this research theme. In this study, the same stiffness values of horizontal and vertical springs are considered due to the continuity of the medium.

Analysis of the beam on a rigid base

A beam of length L = 5 m, a width b = 0.60 m, and height h = 0.80 m, supported on a rigid soil foundation is considered. The load on the beam is 500 kN applied at the middle. For the symmetry reason, only one half of the beam is considered in the analysis (Fig. 13).

Fig. 13
figure 13

Beam on a rigid base

The beam reposing on a rigid base shows important horizontal displacements of upper nodes in the vicinity of the applied load area (Fig. 14) but vertical displacements are very significant at the middle of the beam and decrease at far points (Fig. 15).

Fig. 14
figure 14

Horizontal displacement field

Fig. 15
figure 15

Vertical displacement field

Influence of soil properties on the interaction response

Mechanical properties of the beam and the soil foundation are mentioned on Fig. 16. A fine meshing is selected in the horizontal direction but a uniform meshing is assumed in the vertical direction. The soil foundation dept is of H = 30 m, length L = 54 m, Young’s modulus Es variable and Poisson ratio \(\nu_{\text{s}}^{{}} = 0.30\).

Fig. 16
figure 16

Beam, soil foundation geometry

Figures 17 and 18 show the influence of mechanical properties of soils on the continuum interface responses. When the mechanical properties of the soil are feeble, horizontal and vertical displacements are very important. The enhancement of soil properties of 35.33% is observed when Es increases from 15 MPa to 45 MPa. The maximum horizontal displacement decreases from 1.15 × 10−3 to 0.357 × 10−3 m (Fig. 17) and from 11 × 10−3 to 3.89 × 10−3 m for vertical displacement (Fig. 18). In this context, an improvement of soil behavior is highly observed, the horizontal displacement passes from 0.357 × 10−3 to 0.154 × 10−3 m (43.13%) (Fig. 17) and from 3.89 × 10−3 to 1.85 × 10−3 m (47.55%) for vertical displacement (Fig. 18) when soil properties change from Es = 45–100 MPa.

Fig. 17
figure 17

Horizontal displacement field of the interface level

Fig. 18
figure 18

Vertical displacement field of the interface level

In the same concept, the displacements of upper nodes of the beam are described on Figs. 19 and 20 where horizontal displacements appear with considerable soil properties (Fig. 19) but with a partial independency between soil properties and vertical displacements of beam nodes (Fig. 20). Important displacements using Es = 15 MPa are obtained followed with a considerable bound of displacements with Es = 45 MPa value. Indifferently, the interface medium and the beam behavior are very close when Es varies from 45 to 100 MPa. A unique distinction of the behavior of interface medium and beam nodes is lightly observed for Es = 15 MPa.

Fig. 19
figure 19

Horizontal displacement of the beam

Fig. 20
figure 20

Vertical displacement of the beam

No considerable effect of foundation deepness can be pronounced on the horizontal interaction behavior (Fig. 21). Only a small variation with similitude decay is shown on vertical displacements using different deep values of the soil foundation (Fig. 22).

Fig. 21
figure 21

Horizontal displacements of the interface nodes

Fig. 22
figure 22

Vertical displacements of the interface nodes

Influence of a ballast layer

Nowadays, requirements for strength and stability of railway tracks are increased that is due to train speed and axle load (Petriaev et al. 2017; Sayeed and Shahin 2018). The ballast need to be elaborated for many roles that are: (1) the isolation of structure-borne noise on railway lines in populated areas; (2) the minimization of vibration effects; (3) the stability of railroads and load distribution layer; and (4) provides longitudinal and lateral track support to resist imposed loading from vehicles and thermal rail stress.

In this study, a ballast layer is introduced between the soil foundation and the beam (Fig. 23) and its addition does not present effect on the horizontal and the vertical responses of the interface continuum. The ballast layer has a notable effect on the beam behavior (Figs. 24, 25) and its integration plays a dual reciprocal role; horizontal displacements decrease when the thickness increases (Fig. 24) but vertical displacements increase proportionally with increasing of the thickness of the ballast layer (Fig. 25).

Fig. 23
figure 23

Ballast layer

Fig. 24
figure 24

Horizontal displacements of beam nodes

Fig. 25
figure 25

Vertical displacements of beam nodes

Finally, the diffuse of the displacement fields within the beam and the foundation soil is analyzed. Herein, horizontal displacements have an irregular distribution within the medium between the exciting and the compressive regions (Fig. 26). Maximum vertical displacements are localized at the applied load region and a regular digression occurs in the rest (Fig. 27).

Fig. 26
figure 26

Horizontal displacement field

Fig. 27
figure 27

Vertical displacement field

Conclusions

A procedure to quantify prismatic beams with perfect adhesion to a homogeneous, linearly elastic and isotropic two-dimensional half-space is proposed. Based on the strain energy expressions, shear deformations of the beam element under plane stress state and of plane strain soil foundation are formulated and employed in the analysis.

The numerical results show that the derived formulation leads to a fast convergence and allows studying various problems of soil-structure interaction. A parametric study has been carried out and derived conclusions can be drawn:

  • There is a notable influence of the laterally interaction on the beam-soil foundation behavior.

  • The soil properties have a primordial effect on the beam and on the interaction beam-soil foundation.

  • The deepness of the elastic foundation has a regular effect on the beam-soil foundation and on the beam.

  • The introduction of a ballast layer (height layer) engenders an influence on the beam-foundation interaction in the longitudinal direction.

  • The finite element formulation was established independently of the beam boundary conditions. This approach can be easily used for other boundary conditions of beams.

  • The approach can be considered as issue to nonlinear analysis and vibration analysis of soil-structure interaction due to the impact loads.