ANALYSIS OF BEAMS ON ELASTIC FOUNDATION: THE FINITE DEFFERENCES APPROACH Teodoru I. Bogdan
1
Abstract In the solution of beams on elastic foundation problem, it is usual to use Winkler’s assumption or the concept of modulus of subgrade reaction. The general solution for the Winkler foundation is of limited use, since practical problems (may) involve beams of finite length, or changes in moment of inertia. Moreover, the solution cannot easily adapt to a change of soil’s modulus of subgrade reaction. Because of these shortcomings, the discrete elements methods, among Finite Differences Method (FDM) and Finite Element Method (FEM), are most powerful and popular, are preferred for use, since all kind of contingencies may be accounted for. This paper examines the use of the Finite Differences Method (FDM) for the analysis of elastic beams resting on Winkler medium. To estimate the deformation, soil reaction distribution and internal forces of a continuous footing subjected to external loads, a computer program, based on Matlab code, has been developed. A comparison between FDM, FEM and analytical solutions is also presented.
Keywords beams on elastic foundation, finite differences method, numerical analysis, Winkler’s assumption, continuous beam
1
INTRODUCTION
The acceptance of nu merical analysis in engineering problems is growing. In p articular, the development of numerical analysis and its application to geotechnical problems over the last time have provided geotechnical engineers with an extremely powerful analysis tool. Moreover, the new codes of practice (e.g. Eurocode 7), are not as prescriptive as the older codes and allow the designer to choose an appropriate method of analysis, [3]. Nevertheless, beams on elastic foundation are most usually analysed based on Winkler’s concept in which the soil is treated as a bed of springs. To obtain the theoretical solution of this approach is laborious and classical solutions are not general in their application. Several distinct disadvantages of the classical solution are presented below [1]:
•
Difficult to remove soil effect when footing tends to separate from soil
•
Difficult to apply multiple types of loads to a footing
•
Difficult to change cross section of the footing
•
Difficult to allow for change in subgrade reaction along footing
For these reasons numerical analyses of a beam resting on an elastic foundation are shown in full detail and modelling aspects will be discussed (e.g. discretization dependency). Finally, the results of the numerical analyses are compared with the result of the general solution.
2
BASICS OF FINITE DIFFERENCES FORMULATION
Finite difference schemes provide an alternative route to the conversion of continuum field equations into relationships between discrete numerical values. This method d iscretizes the domain into a regular grid defined by a certain number of nodes which are separated in the coordinates direction by a certain spatial interval. When applying the FDM over the domain, we will be able to approximate the objective function (e.g. displacements) at each one of the nodes. The FDM consists of transforming the partial derivatives in difference equations over a small interval, using the development in Taylor's series of the objective function at each node, leading to an equation that relates the value of the function at a particular node with its value at the neighbouring nodes. This procedure is repeated at each node composing the grid and the assembly of the obtained equations yields a system of equations which may be numerically solved. The method is best illustrated by a physical example, from structural mechanics. Taking the case of a beam in bending, the average slope of the elastic cu rve at distance xi is given by:
1
Teodoru I. Bogdan, Ph.D. student, “Gh. Asachi” Technical Univesity of Ia şi, Faculty of Civil Engineering, Department of Transportation Infrastructure and Foundations, Bd. Dimitrie Mangeron nr. 67, 700050 Iasi IS, Romania, e-mail:
[email protected]
( tan θ )i =
dy
dx i
(1)
where y is deflection. In terms of finite differences, assumed ∆ x = λ as finite values, Equation (1) can be evaluated by three app roximations: 1)
2)
3)
Central difference
dy yi +1 − y i−1 ≈ 2λ dx i
(2)
dy yi+1 − yi ≈ λ dx i
(3)
dy yi − yi−1 ≈ λ dx i
(4)
Forward difference
Backward difference
The last equations are called a finite-difference equation . Solving these equations gives an approximate solution to the differential equation. In an analogous way, we can obtain finite difference approximations to higher order derivatives and differential x + λ /2) and ( x x − λ /2) and then operators. For example, by using the above central difference formula with step λ /2 for ( x applying a central difference formula for the derivative at xi, we obtain the central difference approximation of the second derivative:
d 2 y yi+1 − 2 yi + yi−1 2 ≈ λ2 dx i
(5)
By differentiating Equation (5) one, obtains the central difference expression for the third derivative as
d3 y yi+2 − 2 yi+1 + 2 y i−1 + yi− 2 3 ≈ dx i 2λ 3
(6)
These last two expressions are sufficient to solve the beam resting on an elastic foundations problem by FDM.
3
SOLUTION OF A BEAM ON AN ELASTIC FOUNDATION USING FDM The differential equation of the deflection curve for a bending beam is given b y: EI
d2y dx 2
= −M
(7)
where: -
E is Young’s modulus;
-
I is the moment of inertia;
-
M is bending moment.
By differentiating the above equation, with respect to x and taking into account relationship between bending moment and shearing force, we obtain expression of the sh earing force: EI
d3 y dx 3
= −V
(8)
By substituting Equations (5) and (6) into Equations (7) and (8) we have: EI
( yi+1 − 2 yi + yi−1 ) = − M i
(9)
( yi+2 − 2 yi+1 + 2 yi−1 + y i− 2 ) = − Vi
(10)
λ2 EI 2λ 3
For a foundation beam, by considering Winkler’s concept, one can replace the foundation with a series of concentrated springs on the base of the footing, as shown in Fig.1. The soil pressure at any point on the beam is directly
proportional with beam deflection ( y) and modulus of subgrade reaction ( k ). ). For instance at the point i, soil pressure is related by:
( p ) i = k ⋅ ( y )i
Fig. 1
(11)
Winkler’s concept
By discretizing the beam domain into n elements of ∆ x = λ = constant, and considering pressure distribution of soil given by known function p = f(x) , one can compute the soil reaction against the beam, at each node i = 1 to n + 1. One may use any type of pressure distribution of soil to footing (for details see [4]), but for computational simplicity a stepped pressure distribution is use. Thus, the reactions against the beam of Fig. 2 become:
( R ) k = ( p ) ⋅ λ ⋅ B = k ⋅ ( y ) ⋅ λ ⋅ B for k = 2 to i-1 k k
(12)
1 1 ( R )1 = ( p ) ⋅ λ ⋅ B = k ⋅ ( y ) ⋅ λ ⋅ B 1 1 2 2
(13)
1 1 ( R )i = ( p ) ⋅ λ ⋅ B = k ⋅ ( y ) ⋅ λ ⋅ B i i 2 2
(14)
where B is width of footing.
Fig. 2
Mathematical model for the FD solution for a beam on elastic foundation
By transposing Equation (7) in finite differences, for each node composing the grid, except extremity points (to avoid limit condition), the bending moment at the point k = 2 to i-1 is given by: EI
j=i− 1
λ
( y k−1 − 2 y k + y k+1 ) = − ∑ R j ( k − j) λ − ∑ M P + (− )∑ M 2
where: -
M P is the bending moment of axial loads;
-
M is concentrated moment.
j=1
(15)
Making notation EI notation EI / λ 2 = C and C and arranging, we have: j=i− 1
Cy k −1 − 2Cy k + Cy k+1 +
∑R
j
( k − j) λ = ∑ M P + (− )∑ M
(16)
j=1
Thus we have obtained i-2 equations in i unknown values of displacements y. y. The number of equations with finite differences being less than unknown’s number, we must complete with two equilibrium equations. The sum of the moments with respect to right end, for instance: j=i− 1
∑ j=1
m=1
R j ⋅ ∏ m ⋅ λ =
∑M
(17)
P
m=i− 1
and the sum of forces in the vertical direction: i
∑R = ∑P
(18)
i
1
We have now obtained i simultaneous linear equations in i unknown values of displacements y displacements y.. Cy1 − 2Cy 2 + Cy 3 + R 1 ⋅ λ =
∑(M
) + (− )∑ ( M )1
P 1
Cy 2 − 2Cy 3 + Cy 4 + R 1 ⋅ 2λ + R 2 ⋅ λ =
∑(M
) + (− )∑ ( M )2
P 2
Cy 3 − 2Cy 4 + Cy 5 + R 1 ⋅ 3λ + R 2 ⋅ 2λ + R 3 ⋅ λ =
∑(M
) + (− )∑ ( M )3
P 3
M
(19)
Cy i −2 − 2Cy i−1 + Cy i + R 1 ⋅ ( i − 1) λ + R 2 ⋅ ( i − 2 )λ + ... + R i−1 ⋅ λ = R 1 ⋅ ( i − 1) λ + R 2 ⋅ ( i − 2 )λ + R 3 ⋅ ( i − 3) λ + ... + R i−1 ⋅ λ = R 1 + R 2 + R 3 + ... + R i =
∑ (M
∑(M
) + (− )∑ ( M )i−1
P i −1
) + ( − )∑ ( M )i
P i
∑P
A numerical scheme for solving a system of equations is matrix method. By arranging system of equations (19) for a computer solution by the matrix method, we have:
•
Coefficient matrix:
C + 0, 5CR λ CR λ 1, 5CR λ M 0, 5 ⋅ m ⋅ CR λ 0, 5 ⋅ m ⋅ CR λ 0, 5 ⋅ m ⋅ CR λ 0, 5 ⋅ m ⋅ CR λ 0, 5C R
−2C C + CR λ 2C R λ ( m −1) CR λ ( m −1) CR λ ( m −1) CR λ ( m −1) CR λ CR
0 L 2C Rλ C + C R λ − 2C C 0 0 L L ( m − 2 ) C Rλ 2CR λ C + CR λ − 2C C 0 L L ( m − 2 ) C R λ ( m − 3) C Rλ 2C R λ C + C R λ − 2 C C L L ( m − 2 ) C R λ ( m − 3) C R λ ( m − 4 ) CR λ 2C R λ C R λ 0 L L CR CR CR CR CR C R 0, 5C R L C −2C C + CR λ
0 C −2C
0 0 C
L
L
L
L
0 0 0
0 0 0
0 0 0
(20)
where: -
C R =k B ⋅ ⋅ ⋅λ
-
m is the row’s indices.
The foregoing matrix can be easily carried out by hand, for a low number of division points (in literature [4] it is recommended to use 10 divisions); increasing number of division points leading to a more accurate solution, but in the same time yields difficulty in writing coefficient matrix.
• Free terms vector: is formed by right hand side terms from Equations (19) (moments of all external loads with respect to the each station point i =2 to i-1) i-1) Once the system of equations (19) is solved, one can compute bending moment and shearing force, by back substitution of the beam deflection ( y yi) into Equations (9) and (10) at each division point form the beam domain. Soil reaction against the beam at a point i is given by Equation (11).
4
COMPUTATIONAL EXAMPLE AND COMPARATIVE ANALYSIS
In order to estimate the capabilities and the advantages of the solution described foregoing, let us consider in some detail a continuous footing as a beam on elastic foundation with loading condition shown in Fig. 4. The solution can be adapted to continuous footings with any number of column loads which may include both axial loads as well as moments.
Fig. 4
Assumed loading for Winkler foundation
For the analysis of the longitudinal bending behaviour of the beam, three different methods (FEM, FDM and general method) were performed. In the analytical solution and FE analysis, the beam is divided into 13 elements by 0,4 m (left span), 0,5 m (middle spans) and 0,35 m for the right span. In FD method the beam is also divided into the same number of elements, but of constant length. For the assumed discretization of the given beam, FDM requires the solution of 14 simultaneous equations, thus necessitating the use of a computer. In order to generate the matrix (20) for FD analysis, a computer code listing in Appendix, based on Matlab language, has b een developed. This can be used for any number of division points along beam. Considering for reference the solution obtaining by general method, comparative results obtained by FD and FE analysis are plotted in Fig. 5. From the plot it is evident that either of the FEM or FDM solutions tends, with a good accuracy, to analytical result. The FDM solution yields more appropriate results, regarding negative bending moments, whereas FEM solution gives a good approach for positive moments. For this discretization step, FDM doesn’t give a very accurate solution regarding the positive bending moment, but its accuracy can be easily controlled by changing the grid size. For instance, by increasing five times the number of divisions yields the bending behaviour of the beam plot in Fig. 6.
-80
F DM F EM Ana lytical lytical So lution
] -60 m N k [ -40 t n-20 e m o 0 M g n 20 i d n 40 e B 60
0
Fig. 5
1
2 3 4 Distance along beam [m]
5
6
Graphical comparison of bending behaviour of the beam
FDM computes displacements strictly at predetermined grid points only (unlike FEM, it doesn’t compute displacement functions that can be used to interpolate displacements at the points that are not located at the grid); to obtain an appropriate solution in singularity points (e.g. points of application of external forces), those must be in coincidence with a point from the grid domain.
Referring to beam with loading condition given in Fig. 4, filling its domain with a discretization with step chosen according with foregoing, so with λ = 10, 5, 2, 1, 0,5 (cm) and so on, the values of bending moments are in good agreement with the resuls obtained by FEM for 65…80 division points (see Table 1 ).
Tab. 1
Comparison of bending moments values (in kNm) obtained by FD and FE analyses M+left
M-left
M+middle
M-right
M+right
FDM
52,5
72,6
45,6
39,0
45,1
FEM
52,5
72,7
45,6
39,1
45,0
-80 ] -60 m N k [ -40 t n -20 e m o 0 M g n 20 i d n 40 e B 60
FDM - 65 eleme eleme nts Ana lyti lytica ca l So lutio lutio n
0
Fig. 6
1
2 3 4 Distance along beam [m]
5
6
Discretization step dependency for bending behaviour of the beam
For another discretization steps, different then λ = 10, 5, 2, 1, 0,5 (cm)..., the solution regarding positive bending moments is plotted in Fig. 7. From this, one can see that for a coarse grid the solution is very unstable and the convergence is obtained for a nodes number by 3000 ( λ = 0,1 cm). The discretization with such a huge number of division points (which can increase computing time) can be avoid if one chose different discretization steps, but not in the same time. For instance, if we want to find the value of the left positive bending moment we can perform an analysis with λ = 40 cm (so 15 elements); a discretization step by λ = 30 cm (so 20 elements) will give the value of the middle positive bending moment.
55
] m N50 k [ s t n45 e m o 40 M g n35 i d n e B30
left middle right
0
2
4
6
8 1 0 12 12 14 14 16 16 18 18 20 20 22 22 24 24 26 26 28 28 30 30
Number of divisi divisi on points [x 102 ]
Fig. 7
Discretization step dependency for positive bending moments
The soil pressure distribution evaluated from the numerically computed using FD and FE and pressure distribution obtained by analytical method are compared in Fig. 8. Good agreement between the foregoing is seen, the diagrams of pressure being almost identically for any discretization step.
0
1
Distance Distance al ong beam [m] 2 3 4
5
6
0 FD M FE M Analyti Analytica ca l So lution
50
] a 10 0 P 100 k [ 15 0 e 150 r u s 200 s 20 0 e r p 250 25 0 l i o 30 0 S 300
350 35 0
Fig. 8
5
Soil pressure distribution along beam
CONCLUSION
The results from FD analysis are generally in good agreement with both analytical and FE results. One limitation of the FDM, however, is that it computes displacements at predetermined grid points only and thus accuracy of the solution obtained is affected somehow by grid size. To obtain an app ropriate solution in singularity points (e.g. points of application of external forces), those must be in coincidence with a point from the grid domain; if this can’t be reached one can chose different discretization steps, but not in the same time. Anyway, the FD approach can provide satisfactory prediction of the structural behaviour of beams resting on elastic foundation and thus of use in professional engineering work.
ACKNOWLEDGEMENTS The author is grateful to Professor Mu şat Vasile for his helpful comments, guidance and support in the research work presented in this paper.
Appendix - Computer code to generate coefficient matrix for FD analysis %F.D.M. FOR BEAMS ON ELASTIC FOUNDATION %Coeffs matrix %======================================= %=================== ===================================== ================= % INPUT DATA L=input ('Beam length -- L = '); n=input('Number of elements -- n = '); l=L/n; % << -------- Discretization step <
> k=input('Modulus of subgrade reaction -- k = '); B=input('The width of the footing -- B = '); I=input('Beam''s moment of inertia -- I = '); E=input('Young''s modulus for footing -- E = '); C=E*I/(l^2); CR=k*B*l; %======================================= %=================== ===================================== ================= % Nodes i=n+1; %Domain grid Vl=1:(i-1); VL=Vl*l; %Eq 1 to i-2 for li=1:i-2, for c=1:i, for m=1:i-3, if li-c == -2, MCRl(li,c)=C; elseif li-c == -1, MCRl(li,c)=-2*C; elseif li == c, MCRl(li,c)=C+CR*l; elseif li-c == m, MCRl(li,c)=(m+1).*CR*l; end
end end end Mc=MCRl; %M_end=0 for j=i-1:-1:1 VM0(j)=j*CR*l; VMc=[fliplr(VM0),0]; end %sumFv=0 for j=i:-1:1 VMf(j)=CR; end %Assembly M=[Mc;VMc;VMf]; M(:,1)=M(:,1)*0.5; M(1,1)=C+0.5*CR*l; M(i,i)=0.5*CR; M %< --- Coefficient matrix
Literature [1]
Bowles J. E., - Foundation Analysis and Design, 3rd Ed., New York: MCGRAW-HILL MCGRAW-HILL EDUCATION, 1982, 816 pages, ISBN 19820070661928.
[2]
Muir Wood D., - Geotechnical Modelling, Abingdon: TAYLOR & FRANCIS, 2004, 480 pages, ISBN 0415343046.
[3]
Potts D. M., - Numerical analysis: a virtual dream or practical reality? , GEOTECHNIQUE, Vol. Vol. 53, No. 6, pp. 535572, 2003, ISSN 0016-8505.
[4]
Winterkorn H.F., Fang H.Y. - Foundation Engineering Handbook , New York: VAN NOSTRAND REINHOLD CO., 1975, 751 pages, ISBN 0442295642.
Reviewer Muşat Vasile, Professor, Ph.D., “Gh. Asachi” Technical Univesity of Ia şi, Faculty of Civil Engineering, Department of Transportation Infrastructure and Foundations, Bd. Dimitrie Mangeron nr. 67, 700050 Iasi IS, Romania, Phone: +40 745 574 061, e-mail: [email protected].
