ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION VEDRAN JAGODNIK, GORDON JELENIC and ZÊLJKO ARBANAS about the authors Vedran Jagodnik University of Rijeka, Faculty of Civil Engineering R. Matejcic 3, 51 000 Rijeka, Republic of Croatia Gordan Jelenic University of Rijeka, Faculty of Civil Engineering R. Matejcic 3, 51 000 Rijeka, Republic of Croatia Zeljko Arbanas University of Rijeka, Faculty of Civil Engineering R. Matejcic 3, 51 000 Rijeka, Republic of Croatia Abstract In this paper the deformation of a Bernoulli beam resting on Winkler's soil is reviewed in terms of the mixed finite-element methodology. While the stiffness matrix of the Bernoulli beam problem utilizing the standard displacement-based approach, in which only the displacement field is interpolated, may be alternatively obtained using a mixed-type approach to the absolutely shear-stiff second-order Timoshenko beam (in which the rotation and shear-stress resultant fields are additionally interpolated), the two approaches lead to different Winkler-type soil-stiffness contributions. Furthermore, extending the mixed-type formalism to both of these elements by additionally interpolating the distributed soil-reaction field, the soil-stiffness contributions also differ. In this way four different elements are obtained, with one, two, three or four independently interpolated fields, in which the beam-stifness matrix is equal, but the soil-stiffness matrices are different. It is demonstrated that the displacement-based one-field element is the least convergent, while the mixed-type element with four interpolated fields is the most convergent. Keywords Bernoulli beam, Winkler soil, mixed finite-element method 1 INTRODUCTION Among many other engineering fields, the well-known Bernoulli beam theory [1, 2] finds widespread application in the numerical analyses of slender engineering structures resting on deformable soil. In contrast to beam structures that are subject to known point or distributed external loading, in this case the actual deformation of the structure results in a soil-induced reaction that depends on the actual constitution of the soil. To avoid the very complex (nonlinear, anisotropic, heterogeneous and stress-dependent) behavior of the soil, the subsoil is often modelled by a simpler system called a subgrade reaction model [3]. In [3] Winkler proposed a model that assumes a constant ratio between the contact pressure and the associated deflection of the soil (settlement) defined by the modulus of subgrade reaction ks. Many researchers [4-11] have investigated the modulus of subgrade reaction, and it was found that the geometry, the foundation dimensions and the soil layering below the foundation structure are the most important parameters needed to define this modulus. Terzaghi made recommendations for how to obtain the modulus of the subgrade reaction from a 1 ft x 1 ft rigid plate test placed on diferent soil layers [5]. Biot solved the problem with an infinite beam model on a 3D elastic soil continuum loaded with a concentrated force [4]. The modulus of the subgrade reaction can be measured using diferent experiments such as the plate-load test, the oedometer test, the triaxial compression test and the California Bearing Ratio (CBR) test. The ranges of values for the modulus of the subgrade reaction for typical soil types are given in Table 1. The desired value of the modulus of the subgrade reaction to be used in the present uniaxial beam model, of course, is measured in force-per- length-squared, rather than force-per-length-cubed. To distinguish between the two, the former modulus will be denoted simply as k, and is related to ks via k = Bks with B as the foundation-strip width. ACTA GEOTeCHNICA SLOVENICA, 2013/2 15. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION Table 1. Range of values for the modulus of the subgrade reaction ks [12], Soil ks [kN/m3] Loose sand 4800 - 16000 Medium dense sand 9600 - 80000 Dense sand 64000 - 128000 Clayey medium dense sand 32000 - 80000 Silty medium dense sand 24000 - 48000 Clayey soil: qu < 200 kPa 12000 - 24000 200 < qu < 400 kPa 24000 - 48000 qu > 400 kPa > 48000 (qu = uniaxial compressive strength) In [6] Vesic gave an expression for k as a direct function of the material properties of the soil (Young's modulus Es, Poisson's ratio v) as well as the material and geometric properties of the foundation strip (Young's modulus E, foundation width B, and its second moment of area I) as: k = 0.65 12 E ■ B4 E ■ I 1-v (1) The solution to the problem of a beam resting on Winkler's soil changes considerably with respect to the problem of a beam with no such soil contribution and standard, point-wise, supporting conditions. The resulting deformation ceases to be polynomial and becomes a combination of trigonometric and hyperbolic functions [13]. Using the finite-element method for the problems of beam-soil interaction is as popular as elsewhere in structural analyses, and standard beam finite elements [1, 2] are regularly used in academic and commercial software. In addition, a number of special-type finite elements using a non-polynomial interpolation of the deformation field have been proposed and successfully tested against the exact results, which they are required to reproduce by design [14, 15]. Even though such 'exact' beam-soil finite elements exist, they by no means obviate the need to assess the performance of the standard finite elements based on a polynomial interpolation. For one reason, the former elements are designed for a situation of limited applicability (e.g., Bernoulli beam on Winkler's soil) and there is no guarantee that they will perform better than the standard beam elements when applied to a different problem, e.g., a non-linear beam or soil model. Additionally, the shape functions used in these elements are non-standard and sometimes contain singularities [15] for extreme values of beam-to-soil stifness. With this motivation in mind, after defining the model problem in Section 2, in Section 3 we recall that the solution to the Bernoulli beam problem utilizing the standard displacement-based approach may be completely recovered using a mixed-type approach to the second-order Timoshenko beam, leading to results with a well-defined shear-rigid limit [16]. In Section 4, we investigate the application of the two approaches to the Bernoulli beam resting on Winkler's soil and show that the resulting soil-stiffness contributions become different, even though the beam-stiffness contribution is the same. Additionally, if the distributed soil-reaction field is interpolated independently from the interpolation of the displacement of the beam reference line using a mixed-type technique, we show that both of the above soil-stiffness contributions experience further changes, effectively leading to four different solutions for the soil-stiffness contribution associated with the same beam stiffness. In Section 5 we numerically analyze the results and assess their performance on two simple test examples. 2 THEORETICAL PROBLEM SET-UP Let us consider a straight beam of length L and uniform cross-sectional bending stiffness EI, loaded by a distributed static loading q and an arbitrary point loading resting on a reacting soil with a distributed soil-reaction field f (x), which is proportional to the amount of displacement, i.e., f (x) = kw(x), with k as the modulus of the subgrade reaction. As a result, the beam deects by an amount w(x) and its cross-sections rotate by an amount Q(x) (here assumed as positive in the counter-clockwise direction), from where the stress-couple resultant may be obtained as M(x) =EI&(x), while the shear-stress resultant follows from T(x) = M(x), where the dash (') indicates a dieren-tiation with respect to the longitudinal coordinate x. Figure 1. Straight beam on elastic soil. l6. ACTA G£OT£CHNICfl SLOVENICA, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION 3 SUMMARY OF THE BERNOULLI BEAM THEORY 3.1 displacement-based approach to the bernoulli beam problem In the Bernoulli beam theory (see, e.g., [17]) the cross-sectional rotations are the same as the rotations of the beam centroidal line (d(x) =-w'(x)) giving the stress-couple resultant in a cross-section as M(x) = -EIw"(x). Equilibrium is achieved when the total potential energy of the problem V = Vdef- U is station- pL ary, where Vief = 2 J EIw"2dx is the strain energy and U = JL qdx + U0L is the work of the applied loading with U0L = w(0)F0 + w '(0)M0 + w(L)FL + w' (L)Ml as the work of the boundary point forces F0, Fl and moments M0, ML. By dividing the beam length into Nei finite elements of length li = xi+1 - xi > 0 each, and assuming a distribution of the displacement field within each element using the standard Hermitean polynomials, the above becomes Nel Edpt Kp> - R ) = 0 (2) i=1 with pi as the vector of the standard nodal degrees of freedom (vertical displacement and rotation at both ends), Spi as its variation, K - EI K - T 12 —61 -12 -61 -61 4/2 61 212 — 12 61 12 61 —61 212 61 412 (3) as the element stiffness matrix and Ri as the corresponding element load vector. tional V* = Vd + VdT - VT - U with Vd - - f EI9'2dx, T L ^ 2 J L 2 VdT =J(w' + q)Tdx and VT =1J~dx. 0 0 0 Assuming a quadratic Lagrangian interpolation for the displacements and the rotations (I1 (X) = -X , 12(X) = 1 - x2and 13(x) = x ) and a linear shear-stress resultant field (T = NtT1 with Nt = (1 X)) thus leads to N" E[öPt (KdiPi + KmT — R.) + ST! (pt + KtT)] - 0 (4) with K d, i - - 2 r dN! f h —1 dX 0 0 0 EI K' (N fdX (5) dX NTdX (6) kt , i - 2GA f NT NTdX (7) where N (X)- 11 (X) 0 12 (X) 0 13 (X) 0 0 11 (X) 0 12 (X) 0 13 (X) Since (4) must hold true for any variations Spi and STi it turns out that for any discontinuous interpolation for the shear stress resultant field (which is admissible owing to the fact that no derivatives of T with respect to x appear in V* ) the term within the second parentheses must vanish. This results in 3.2 bernoulli beam as a shear-rigid second-order timoshenko beam The above result may be reproduced using the Timosh-enko beam theory (see, e.g., [17]) with an infinite shear stiffness. In the Timoshenko beam theory, the cross-sectional rotation is in general assumed to differ from the amount of rotation of the beam centroidal line by a shear angle y(x) = 6(x) + w'(x), resulting in a shear-stress resultant T(x) = GAy(x) where G is the shear modulus of the beam material and A is the shear area of the cross-section. We consider the shear-stress resultant T as an independent field and, instead of V = Vdey- U, we start from the condition of stationarity of a mixed func- * • "i £ s p! Kd,i — KdT,iK— jKdT,i )pi — R, - 0 (8) with Kb. as the stiffness matrix of a three-node element in which the displacement and the rotation at the internal node may be condensed out to eventually give Kb.i - EI (1 + f)) 12 ~6lt —12 ~6lt —61, (4 + f)) 6/i (2- —f)2 — 12 6/i 12 6lt —61, (2- —f)/2 6/i (4 +f)i , (9) ACTA GeOTeCHNICA SLOVeNICA, 2013/2 17. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION with f: 12EI GAlJ , which coincides with the result obtained in [16] using the stiffness-based approach. This result can also be obtained by consistently deriving the appropriate shape functions needed to obtain the exact solution, as shown by Reddy [18]. This stiffness matrix has a well-defined shear-rigid limit that coincides with (3). K Si = kl, 13 35 210 ! 9 70 13 420 -i! l, 210 i ^ l2 105i --1, 420i l2 140 i ^ l, 420i 140i 11 9 70 13 420 13 35 11 210 210 i ^ l2 105i (13) 4 BERNOULLI BEAM ON 1 I WINKLER S SOIL The differential equation for a Bernoulli beam resting on Winkler's soil (see, e.g., [19]) is given as EIwIV (x) + kw(x) = q(x) (10) defined over a domain 0 < x < L with known values for w or its derivatives at x = 0,L as the boundary conditions. For known boundary conditions, the above dierential equation is easily solved [13]. 4.1 displacement-based approach to the original bernoulli-winkler problem (one-field interpolation) The problem may be variationally approached in a manner completely analogous to that presented in Section 3.1 with the only dierence being that the strain energy of the problem is now 1 L 1 L vdef = - J (EIw"2 + kw2) dx rather than - J EIw " 2dx. 2 0 0 From SV = 0, therefore, it now follows (Kipi - R ) = 0 (11) ,=1 with K, = Kb, + Ks,i and 4.2 mixed approach to the original bernoulli-winkler problem (two-field interpolation) Instead of minimizing the total potential energy of the problem, it is possible to approach the task from the standpoint of the Hellinger-Reissner complementary energy principle [1]. In particular, let us introduce a new function f(x) (the distributed soil-reaction field), which, although obviously uniquely related to the unknown displacement function w(x), serves to define a new two-field potential V* = V + Vf - Vf - U, where 1 l l 1 L Vb = - J EIw"2dx, Vbf = J fwdx and Vf =— J f 2dx in 00 0 which w(x) and f(x) may be treated (and interpolated) as independent fields. Furthermore, owing to the absence of any derivatives on f(x), this function only requires C"1 continuity, i.e., it does not have to be continuous between the elements. From SV* = 0 we therefore obtain L L L 6Vb + J Swfdx + J 6fwdx - J 6ffdx - SU = 0 . (14) Dividing the beam into Nej finite elements as before, and now additionally assuming the interpolation of the distributed soil-reaction field within each element as ft (X) = Nf i (X) ft with ft = f... fM ) where M-1 is the order of the polynomial used to describe the distributed soil-reaction field, the above equation becomes /kl t -fN NdX , (12) E|SPÎ (( P + K bf,, f - R ) + 6f> (( p + K f,, f = 0 (15) where Ks,j denotes the soil part of the stiffness matrix. Using the standard Hermitean polynomials N. = ,2(1-x)-x(1-x2) - l,(1 -x2)(1-x) v 4 8 (1+x)+x(1-x2 ) l (1-x2)(1+x). this matrix is computed as with Kbf ,i = ~ JNtNf,, dx (16) l Kf,=-iJNf,Nf,dx . (17) -1 Since (15) must hold true for any variations Spj and Sf it turns out that for any discontinuous interpolation for the N N =1 18. ACTA GeOTeCHNICfl SLOVeNICfl, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION distributed soil-reaction field the term within the second parentheses must vanish. As a result f =-K-, X/ i , pt (18) and (15) becomes Ni E « P i=i and N/ i (X) = (1 X X2 X3 ) all make the possible choices. For a constant approximation of the distributed soil-reaction field, the soil-stiffness contribution to the element stiffness matrix thus follows as K = klt 144 36 ~6li 36 6lf —6li l,2 — 6li —2 36 —6lf 36 6lf 61 — l2 61 l2 (21) while for a linear approximation of the distributed soil-reaction field, it becomes K= kl 1800 666 —93lt 234 57li —93lt 14/2 — 57li — 11l2 234 — 57li 666 93li 57li — 11l2 931, 14l2 (22) using the displacement-based approach (13), in which the displacement field, from which the distributed soil-reaction is computed, is cubic. 4.3 bernoulli-winkler beam as a second-order shear-rigid timosh-enko beam on winkler's soil (three-field interpolation) The problem can also be variationally approached in a manner completely analogous to that presented in Section 3.2, with the only difference being that the original three-field mixed functional V* is now substituted 1 1 L with V* + - J kw2dx . From SV* + J Swkwdx = 0, 2 0 0 therefore, it now follows E «Pt [(Kb, +KsJ ) — R ] = 0 (24) with kl, K Si = — Si' 2 / — 1 N1 1 0 0 0 NdX . (25) The soil contribution to the element's stiffness matrix for a shear-rigid element (GA upon condensation of the mid-node degrees of freedom reads K. = ki s'' 6 l l 2 —L 1 -j- l 4 l2 l 4 l2 4 20 4 20 1 I 2 l l 4 l2 l 4 l2 4 20 4 20 (26) For a quadratic approximation of the distributed soil-reaction fieldi the soil-stiffness contribution to the element stiffness matrix becomes Ks,i = kl 1800 666 —93li 234 57li —931, 16.5l2 —57li — 13.5l? 234 — 57li 666 93 li 57li — 13.5l2 931, 16.5l? i (23) while for a cubic and all further higher-order approximations of the distributed soil-reaction field, the soil-stiffness contribution to the element stiffness matrix is, as expected, identical to the one obtained It should be noted by comparing (13) and (26) that even though the shear-rigid second-order Timoshenko beam reproduces the stiffness matrix of the Bernoulli beam, the same is not true for the ensuing soil part. 4.4 dual mixed approach to the shear-rigid timoshenko-winkler problem (four-field interpolation) Obviously, the mixed approaches described in Sections 3.2 (where displacement, rotation and shear-stress resultant have been interpolated as independent fields) and 4.2 (where displacement and distributed soil-reaction field have been interpolated as independent fields) can be easily ACTA GeOTeCHNICA SLOVeNICA, 2013/2 19. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION combined to form a four-field functional V** =Vd + VdT - Vt+Vbf - Vf - U . From SV** = 0 we now obtain N" r / \ EdPt ((P, + KdTJ, + Kbf jfi - R)+ (27) 1=1 + ST< (,i p,+ Kt T) + f (K f,, p, + Kf,, f,,)] = 0 which, after condensing out the shear-stress resultant and distributed soil-reaction parameters, turns into EN=15Pt [(Kb,, + Ks,)p, -R,] = 0 with (28) Ks,i — Kbf ,iKf1 Kbf ,i . It should be noted that the soil part of the stiffness matrix (28) is not necessarily the same as in (20) for the corresponding choice of interpolation for the distributed soil-reaction field, since the interpolation for the displacement field is different. For a three-node shear-rigid beam element with a constant and a linear approximation of the distributed soil-reaction field, the soil part of the stiffness matrix upon condensation of the mid-node degrees of freedom becomes, respectively, Ksi — kli 144 36 36 9l. —9l, 36 9l 912 —9l; — 912 4 -9li —912 4i 36 91 4 911 912 4i and (29) 48 -6li 24 6li —6lt l2 —6lt —l2 24 -6l 48 6l 6lt f2 6lt l2 kli 144 while for the quadratic and higher-order approximation of the distributed soil-reaction field we obtain the same result as in (26). 5 NUMERICAL EXAMPLES In this section, two problems will be analyzed: a rectangular beam on an elastic foundation with the concentrated load in the middle of the span and a beam on an elastic foundation with the concentrated force at a boundary, thus simulating a laterally loaded pile. For each example a displacement convergence analysis is undertaken for different values of the relative soil-kL4 stiffness parameter b = . 4EI Six values of the relative soil-stiffness parameter ft have been considered: 1, 5, 10, 50, 100 and 500. Using the expression for ft and back-calculating the desired value for the modulus of the subgrade reaction k, and then plugging it into the equation for the modulus of the subgrade reaction proposed by Vesic [6] (1), approximate values of the Young's modulus can be obtained, which makes it possible to classify the soils in a range between soft and hard. Which soil the particular value of the relative soil-stiffness parameter ft is describing is given in Table 2. Table 2. Soil description based on the relative soil-stiffness parameter ji. ß Soil Description 1 Very loose sand 5 Stiff clay 10 Medium dense sand 50 Very dense gravel 100 Weak rock mass 500 Hard rock mass 5.1 thin germ on winkler s foundation subject to a concentrated force in the middle The problem is sketched in Figure 2. The material and geometric properties of the beam and loading are summarized in Table 3. Table 3. Beam parameters. Length 3.0 m Width 0.3 m Young's modulus 3.15 • 107kPa Second moment of area 0.000675m4 Force 1.0 kN Due to its symmetry, only one half of the problem is analyzed with the rotation at the symmetry line being set to zero. In the absence of the soil stiffness, all the approaches would return exact nodal results for the displacements and rotations, provided a suitable supporting (e.g., a simple support at the free end) were defined. The results for the displacement in the middle of the beam will be normalized with respect to the exact solution provided by Hetenyi [13], expressed in terms of the parameter ft: l6. ACTA G£OT£CHNICfl SLOVENICA, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION Figure 2. Mesh of n finite elements for a beam on Winkler's foundation with a point force in the middle. P ■ L (b) + cosh((b) + 2 8 ■ b3/4 ■ EI sin ( (b) + sinh ( (b (30) 5.1.1 Displacement convergence analysis We first analyze the convergence properties of the different formulations and in Table 4 show the normalized displacements in the middle of the beam for meshes with two, four, eight, sixteen, thirty-two and sixty-four elements for a soft, moderately hard and hard soil (3 = 5, 3 = 50 and 3 = 500, respectively). It is immediately obvious that a constant approximation of the distributed soil-reaction field (equations (21) and (29)j) is inferior to the linear (equations (22) and (29)2) or quadratic interpolation (23). The use of the quadratic approximation of the distributed soil-reaction field in Table 4. Values of normalized displacements at mid-span for 3 = 5, 50 and 500. Formulation Elements 3 = 5 3 = 50 3 = 500 2 0.998348089607226 0.987206820721338 0.919540084654748 4 0.999897498365916 0.999168178332181 0.991938846421126 One-field (13) 8 16 0.999994057876285 1.000000000000000 0.999948078751554 0.999996754921972 0.999469328447786 0.999962094889128 32 1.000000000000000 1.000000000000000 0.999993682481521 64 1.000000000000000 1.000000000000000 0.999993682481521 2 1.000506566046710 1.104105348213100 2.838429464906180 4 1.000401093350760 1.020562977770130 1.100707562069620 Two-field constant (21) 8 16 1.000120328005230 1.000031196149500 1.005037442791980 1.001255845196800 1.022281887674520 1.005489923558030 32 1.000008913185570 1.000313690876030 1.001370901509890 64 1.000002971061860 1.000077881872670 1.000341145997850 2 0.998404539782518 0.989357225761160 0.925636489986733 4 0.999903440489631 0.999378026711319 0.993031777117948 Two-field (linear) (22) 8 0.999994057876285 0.999962140756341 0.999583043780403 16 1.000000000000000 0.999997836614648 0.999974729926085 32,64 1.000000000000000 1.000000000000000 0.999993682481521 2 0.998351060669083 0.987360421081325 0.921144734348348 4 0.999897498365916 0.999172505102885 0.992039926716786 Two-field quadratic (23) 8 0.999994057876285 0.999948078751554 0.999469328447786 16 1.000000000000000 0.999996754921972 0.999962094889128 32,64 1.000000000000000 1.000000000000000 0.999993682481521 2 0.998599144334185 0.998335274971687 1.010076441973590 4 0.999910868144275 0.999731740216360 0.997757280940047 Three-field (26) 8 16 0.999995543407214 1.000000000000000 0.999979447839157 0.999998918307324 0.999772569334765 0.999981047444564 32 1.000000000000000 1.000000000000000 0.999993682481521 64 1.000000000000000 1.000000000000000 1.000000000000000 2 0.986178620238873 0.999778253001427 2.229098490113080 4 0.996439182363777 0.993338936501395 1.030696822288210 Four-field constant (29)1 8 16 0.999108681442747 0.999777170360687 0.998246576172258 0.999561914466233 1.004378040305770 1.001029755512030 32 0.999945035355636 0.999890749039727 1.000252700739150 64 0.999986630221641 0.999971875990425 1.000063175184790 2 0.998651137916691 1.000393736134050 1.013683745024950 4 0.999915324737061 0.999937261824794 0.998742813822731 Four-field (linear) (29)2 8 16 0.999995543407214 1.000000000000000 0.999994591536620 0.999998918307324 0.999892602185862 0.999987364963043 32 1.000000000000000 1.000000000000000 0.999993682481521 64 1.000000000000000 1.000000000000000 1.000000000000000 ACTA GeOTeCHNICA SLOVeNICA, 2013/2 21. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION the two-field approach gives marginally less-convergent results than the use of the linear approximation of this field and in what follows only such a linear approximation of the distributed soil-reaction field will be analyzed. The corresponding formulations (22) and (29)2 will be denominated as the two-field formulation and the four-field formulation, respectively. The results from Table 4 for the one-field, two-field, three-field and four-field formulations are shown graphically in Figures 3, 4 and 5. From these figures it can be easily concluded that the three-field and four-field element formulations behave visibly better than the one-field and two-field element formulations. £ yv jy yv* /T-y /// J' /> // "s -NiijiiIht of 1 Uments Figure 3. Normalized displacements in the middle of the beam depending on the soil stiffness for various elements (fi = 5). Additionally, it can be observed that the four-field element converges better than the three-field element and, apart from the two-element case with fi = 500 (hard soil), is also the most accurate one, regardless of the mesh size. Likewise, the two-field element is more accurate than the one-field element. It is interesting to note that the mixed character of the formulations becomes increasingly pronounced as the soil gets stiffer, whereby, for the three-field and four-field formulations, the convergence ceases to be monotonous from the stiff side. The superiority of the three-field and four-field formulations compared to the one-field and two-field formulations increases as the soil gets stiffer. 5.1.2 parametric analysis for different soil-to-beam stiffness ratios We further compare the performance of the four formulations for a fixed eight-element mesh and various soil stiffnesses, including the results for fi = 1, fi = 10 and fi = 100. As an illustration, we also give the results for some other typical soils, given by the values fi = 3, fi = 25, fi = 35 and fi = 60. Figure 6 shows that between themselves the formulations rank in performance in direct correspondence to the number of fields interpolated, with the four-field formulation consistently the most effective and one-field formulation the least effective. As already noted, the difference in performance between the formulations increases as the soil gets stiffer, but now it can be seen that this goes hand-in-hand with the Figure 4. Normalized displacements in the middle of the beam depending on the soil stiffness for various elements (fi = 50). Figure 5. Normalized displacements in the middle of the beam depending on the soil stiffness for various elements (fi = 500). 22. acta GeOTeCHNICA SLOVeNICA, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION Figure 6. Normalized displacement in the middle of the beam for 8 elements varying with soil stiffness. general deterioration of the performance. To understand why this is so, it should be noted that any finite-element formulation based on a polynomial interpolation is bound to lose the absolute accuracy as the soil gets stiffer, as in this case the exact solution is dominated by the trigonometric and hyperbolic functions. 5.2 thin pile in winkler s soil subject to a horizontal force at the head In this section, a beam will be subjected to an end force. Such a set-up can be thought of as a pile (which is essentially a beam embedded in soil) subjected to a horizontal load at the head. This is shown in Figure 7. The material and geometric properties as well as the value of the loading are the same as in the previous model problem (Table 3). The results for the displacement at the head of the beam will be normalized with respect to the exact solution provided by Hetenyi [13], expressed in terms of the parameter ft: Figure 7. Horizontally loaded pile in Winkler's soil - geometric and material parameters. ACTA GeOTeCHNICA SLOVeNICA, 2013/2 23. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION _ P. L3 sinh(2^b )- sin(2^b ) w = 2. b3,4. EI ' cos(2Vb ) + cosh(2Vb )- 2 ' (31) 5.2.1 Displacement convergence analysis In this example, the constant approximation of the distributed soil-reaction field is also considerably less effective than the linear approximation, while the quadratic interpolation of that field is marginally less effective than the linear approximation, hence these results are omitted from any further discussion. The results for the one-field, two-field, three-field and four-field formulations (equations (13), (22), (26) and (29)2) for the finite-element meshes of two, four, eight, sixteen, thirty-two and sixty-four elements for a soft, moderately hard and hard soil (fi = 5, fi = 50 and fi = 500, respectively) are given in Table 5. Table 5. Values of normalized displacements at the head of pile for ft = 5, 50 and 500. Formulation Elements ß = 5 ß = 50 ß = 500 2 0.998641145803926 0.990463886169159 0.967540844437416 4 0.999898368245723 0.999206289432062 0.993935963912077 One-field (13) 8 16 0.999992471721906 0.999996235860953 0.999947858429844 0.999997103246102 0.999523932385082 0.999970037702558 32 1.000000000000000 1.000000000000000 0.999998335427920 64 1.000000000000000 1.000000000000000 1.000000000000000 2 0.998701372028683 0.993577896609060 1.010959542575590 4 0.999902132384770 0.999409062204893 0.996364574576991 Two-field (22) 8 0.999992471721906 0.999962342199332 0.999685395876855 16 0.999996235860953 1.000000000000000 0.999980025135038 32,64 1.000000000000000 1.000000000000000 1.000000000000000 2 0.998678787194399 0.993928403830667 0.961158875082188 4 0.999905896523818 0.999594454454339 0.997574718479247 Three-field (26) 8 0.999992471721906 0.999973929214922 0.999838536508227 16 0.999996235860953 1.000000000000000 0.999990012567519 32,64 1.000000000000000 1.000000000000000 1.000000000000000 2 0.998742777558203 0.997030827254978 0.999930087972634 4 0.999909660662865 0.999797227227169 0.999978360562958 Four-field (29)2 8 0.999992471721906 0.999988412984410 0.999998335427920 16 0.999996235860953 1.000000000000000 1.000000000000000 32,64 1.0000000000000000 1.000000000000000 1.000000000000000 g 0.9996 / S /// J' Ay ß <&/ V / / / / Number of Elements Number of Elements -Four - field (2S)i ---— Three - field (26) FE convergence plot — JS = 5 ........Two - field (22) ---- One - field /13) Figure 8. Normalized displacements at the head of the pile depending on the soil stiffness for various elementss (ft = 5). Figure 9. Normalized displacements at the head of the pile depending on the soil stiffness for various elementss (ft = 50). l6. ACTA G£OT£CHNICfl SLOVENICA, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION 0 ! a a. 0,99 1 — I ffl E c £ , / - - // " / / 4 VmuhiTof Klc-mcnrs fl cokiwwahca plor fi jt)lt -Four - field (23b ----Thin: fold (26J ........Two field (22) ---One - field ¡131 Figure 10. Normalized displacements at the head of the pile depending on the soil stiffness for various elementss (fi = 500). These results are shown graphically in Figures 8, 9 and 10. It can be concluded that, as in the previous example, the three-field and four-field element formulations behave better than the one-field and two-field element formulations and this is especially true for the values of fi of 50 and 500. Additionally, it can be observed that the four-field element converges better than three-field element, especially for the case of fi = 500 where it gives results very close to the exact displacement. As in the previous example, the mixed character of the formulations becomes increasingly pronounced as the soil gets stiffer. As in the previous example, the superiority of the three-field and four-field formulations increases as the soil gets stiffer. 5.2.2 parametric analysis for different soil-to-pile stiffness ratios Again, we also compare the performance of the four formulations for a fixed eight-element mesh and various soil stiffnesses, including the results for fi = 1, fi = 10 and fi = 100. As before, we now also give the results for some other typical soils, given by the values fi = 3, fi = 25, fi = 35 and fi = 60. Figure 11 shows that the formulations again rank in performance as in the previous example, with the four-field formulation being by far the most effective. Figure 11. Normalized displacement at the head of the pile for 8 elements varying with soil stiffness. ACTA GeOTeCHNICA SLOVeNICA, 2013/2 25. V. JAGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION The difference in performance between the formulations again increases as the soil gets stiffer, but now it can be seen that, in contrast to the previous example, the four-field formulation gives surprisingly good results, even for very hard soils. 6 CONCLUSIONS The classic problem of a Bernoulli beam resting on Winkler's soil has been reviewed in the light of a mixed finite-element design methodology. It has been recalled that the same standard stiffness matrix of the Bernoulli beam may be obtained either from the classic interpolation of the displacement field using Hermitean polynomials (one-field approach) or as the shear-rigid limit of the second order mixed-type Timoshenko beam in which displacements, rotations and shear-stress resultants are interpolated independently (three-field approach). Still, applying these interpolations to the Bernoulli-Winkler problem results in different respective soil-stiffness contributions. It has been shown on the two model problems analyzed (a beam on an elastic foundation loaded by a central force and a pile horizontally loaded at the head) that the three-field approach is superior to the one-field approach. Additionally, in both of these finite-element design principles we have investigated the possibility to interpolate the distributed soil-reaction field as an independent field, thus effectively obtaining what we have named the two-field and the four-field approaches. We have analyzed the performance of the finite elements in which this field is approximated by polynomial functions of various orders (constant, linear and, in the two-field approach, quadratic) and concluded that a linear approximation gives the best results. In both cases, interpolating the distributed soil-reaction field independently improves the results of the underlying formulation, i.e., the two-field approach turns out to be more effective than the one-field approach, while the four-field approach is more effective than the three-field approach. In both the numerical problems analyzed, the four-field approach gives the best results for the displacements regardless of the soil stiffness. The results are surprisingly accurate for the horizontally loaded pile embedded in hard soil. acknowlcdgcmcnts This work has been financially supported by the Ministry of Science, Education and Sports of the Republic of Croatia through the Research Project 114-0000000-3025 (Improved accuracy in non-linear beam elements with finite 3D rotations). REFERENCES [1] Bathe K. J. (1985). Finite Element Procedures, Prentice-Hall, New Jersey. [2] Zienkiewicz O. C., Taylor R. L. (2005). The Finite Element Method for Solid and Structural Mechanics, Elsevier Butterworth-Heinemann, Oxford. [3] Winkler E. (1867). Die Lehre Von Der Elastizizat Und Festigkeit, Dominicus, Prague. [4] Biot M. A. (1937). Bending of infinite beams on an elastic foundation, Journal of Applied Mechanics, 57, A, pp. 1-7. [5] Terzaghi K. (1955). Evalution of coefficients of subgrade reaction, Geotechnique, 5, pp. 297-326. [6] Vesic A. (1961). Beams on elastic subgrade and the Winkler's hypothesis, Proceedings of 5th International Conference on Soil Mechanics and Foundation Engineering, vol. 1, Michigan, pp. 845-850. [7] Horvath J. S. (1983). Modulus of subgrade reaction: new perspective, Journal of Geotechnical Engineering, 109, pp.1591-1596. [8] Horvath J. S. (1989). Subgrade models for soil-structure interaction analysis, Foundation Engineering: Current Principles and Practices, ASCE, Illinois, pp. 599-612. [9] Daloglu A. T., Vallabhan C. V. G. (2000). Values of k for slab on Winkler foundation, Journal of Geotechnical and Geoenvironmental Engineering, 126, pp.463-471. [10] Vallabhan C. V. G., Das Y. C. (1988). Parametric study of beams on elastic foundations, Journal of Engineering Mechanics, 114, pp. 2072-2082. [11] Vallabhan C. V. G., Das Y. C. (1991). Modified Vlasov model for beams on elastic foundations, Journal of Geotechnical Engineering, 117, pp. 956-966. [12] Bowles J. E. (2001). Foundation Analysis and Design, McGraw-Hill, New Jersey. [13] Hetenyi M. (1964). Beams on Elastic Foundation, University of Michigan Press: Ann Arbor, Michigan. [14] Ting B. Y., Mockry E. F. (1984). Beam on elastic foundation finite element, Journal of Structural Engineering, 110, pp. 2324-2339. [15] Eisenberger M., Yankelevsky D. Z. (1985). Exact stiffness matrix for beams on elastic foundation, Computers and Structures, 21, pp. 1355-1359. [16] Przemieniecki J. (1968). Theory of Matrix Structural Analysis, McGraw-Hill, New York. l6. ACTA G£OT£CHNICfl SLOVENICA, 2013/2 V. lflGODNIK ET AL.: ON THE APPLICATION OF A MIXED FINITE-ELEMENT APPROACH TO BEAM-SOIL INTERACTION [17] Timoshenko S. P. (1958). Strength of Materials, Part 1: Elementary Theory and Problems, Van Nostrand Reinhold, New York. [18] Reddy J. N. (1997). On locking-free shear deform-able beam finite elements, Computer Methods in Applied Mechanics and Engineering, 149, pp. 113-132. [19] Timoshenko S. P. (1958). Strength of Materials, Part 2: Advanced Theory and Problems, Van Nostrand Reinhold, New York. ACTA GeOTeCHNICA SLOVeNICA, 2013/2 2.