# Bending Moment Calculations for Piles Based on the Finite Element Method

Authors: Yu-xin Jie; Hui-na Yuan; Hou-de Zhou; Yu-zhen Yu
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

Using the finite element analysis program ABAQUS, a series of calculations on a cantilever beam, pile, and sheet pile wall were made to investigate the bending moment computational methods. The analyses demonstrated that the shear locking is not significant for the passive pile embedded in soil. Therefore, higher-order elements are not always necessary in the computation. The number of grids across the pile section is important for bending moment calculated with stress and less significant for that calculated with displacement. Although computing bending moment with displacement requires fewer grid numbers across the pile section, it sometimes results in variation of the results. For displacement calculation, a pile row can be suitably represented by an equivalent sheet pile wall, whereas the resulting bending moments may be different. Calculated results of bending moment may differ greatly with different grid partitions and computational methods. Therefore, a comparison of results is necessary when performing the analysis.


---

## Body

## 1. Introduction

As the finite element method (FEM) develops, pile foundations are increasingly being analyzed using FEM [1–8]. Solid elements are used to simulate soil or rock in geotechnical engineering. Other structures embedded in soil such as piles, cut-off walls, and concrete panels are also often simulated with solid elements. However, internal force and bending moment are generally used for engineering design. So it is necessary to calculate the bending moment with stress and displacement obtained using FEM.Theoretically, the following two methods are both appropriate.(a) Calculating Bending Moment with StressThe bending moment is directly calculated by summing the total moments of the elements across the specified pile section. When using this method, sufficient grids are necessary to partition the pile section.(b) Calculating Bending Moment with DisplacementThe bending moment is indirectly calculated using the quadratic differential of deflection (lateral displacement) of the pile. This method uses fewer grids, but the differential process will result in reduced accuracy.The bending moment can also be obtained by integrating the area of the shear force diagram [9] which is a complex process and is not considered in this paper.As we know, shear locking occurs in first-order (linear) fully integrated elements that are subjected to bending, while second-order reduced-integration elements can yield more reasonable results in this case and are often used in the analysis of piles subjected to lateral pressure [1–4, 10]. However, calculating second-order elements is time consuming and increases complexity and computational effort, particularly when the problem involves contact conditions. So we consider that the linear element method with appropriate meshing is still useful for the analysis of piles.A row of piles can be simplified as a plane strain wall (sheet pile wall) and modeled using 2D plane strain elements [11–13]. This simplification can greatly reduce computational effort. However, the influence of bending moment on the computational results merits further research.In this paper, a series of calculations on cantilever beam, pile, and sheet pile wall examples were conducted to study the abovementioned problems. The main aim of the work was to investigate the computational methods for bending moment and the influences of element type and mesh partition. Hence, no interface element was introduced, that is, the pile was assumed to be fully attached to the soil, and the soil and pile were both assumed to have linear elastic behavior.


## 2. Cantilever Beam Example

### 2.1. Analytical Solution

The cantilever beam example is shown in Figure1. The width of the square beam is 1 m. The length L is 30 m. A distributed load p=0.5 kPa is applied to the beam. The analytical solution equations are(1a)M=12p(L-x)2,(1b)ω=pL42EI(k22-k33+k412),where M= bending moment, x= the position coordinate, k=x/L, E= the Young’s modulus, I= the moment of inertia, ω=-uy is the deflection of the beam, and uy= displacement in the y direction.Figure 1
Cantilever beam.The beam parameters are taken as Young’s modulusE=2×104 MPa and Poisson’s ratio ν=0.17 in the computation. The element used in the FEM is the 4-node first-order plane stress element (CPS4). The following two methods were used to calculate the bending moment.(a) Calculating Bending Moment with StressThe bending moment was directly computed with the normal stress on the cross-section (see Figure2):
(2)M=∑σiAili,
where σi= normal stress at the centroid of the element, Ai= corresponding area of the element, and li= distance between the centroid and the midline of the beam section.Figure 2
Computing of bending moment with stress.(b) Calculating Bending Moment with DisplacementThe bending moment was calculated using the following quadratic differential of deflection [14]:
(3)M=EId2ωdx2.Equation (3) can be transformed into a difference scheme, and the bending moment was calculated by the difference operation of lateral displacement. Figure 3 is the comparison of the computed bending moment with analytically exact results, where ω is calculated using (1b) with the number of interpolation points being 10, 30, 60, and 120. This shows that the method used here achieved a good result for the bending moment.Figure 3
Calculation versus analytical solution of bending moment.


### 2.2. Errors in Displacement and Bending Moment

We used several meshes to do the computation. Grids along the length of the beam were 32, 64, and 128 (in thex direction, see Figure 1). Sectional partitions were 1, 2, 4, 8, 16, and 32. Two element types, CPS4 and the 8-node reduced-integration element (CPS8R), were employed in the analysis. The results are shown in Tables 1, 2, and 3 and in Figures 4 and 5, where the term m×n denotes that there are m grids along the y direction and n grids along the x direction; B, C, D, and E are positions for error calculation (see Figure 1).Table 1
Deflection error with element CPS4.
Mesh
Deflection error (%)
B
C
D
E
1
×
32
−28.57
−28.69
−28.74
−28.76
4
×
32
−27.20
−27.31
−27.34
−27.36
8
×
32
−27.12
−27.23
−27.27
−27.29
1
6
×
32
−27.10
−27.22
−27.25
−27.27
32
×
32
−27.10
−27.21
−27.25
−27.26
1
×
64
−10.67
−10.83
−10.89
−10.92
4
×
64
−8.45
−8.61
−8.66
−8.68
8
×
64
−8.32
−8.48
−8.54
−8.57
16
×
64
−8.29
−8.45
−8.54
−8.54
32
×
64
−8.28
−8.45
−8.50
−8.53
64
×
64
−8.28
−8.44
−8.50
−8.53
1
×
128
−4.70
−4.87
−4.93
−4.97
4
×
128
−2.12
−2.31
−2.37
−2.40
8
×
128
−1.97
−2.17
−2.23
−2.27
16
×
128
−1.94
−2.13
−2.20
−2.23
32
×
128
−1.93
−2.13
−2.19
−2.23
64
×
128
−1.93
−2.13
−2.19
−2.22
128
×
128
−1.92
−2.13
−2.19
−2.22
2
×
4
−95.94
−95.94
−95.93
−95.92
2
×
8
−85.66
−85.68
−85.68
−85.67
2
×
16
−60.03
−60.08
−60.09
−60.09
2
×
32
−27.50
−27.60
−27.63
−27.65
2
×
64
−8.93
−9.08
−9.12
−9.15
2
×
128
−2.69
−2.86
−2.91
−2.94
2
×
256
−0.99
−1.16
−1.22
−1.25Table 2
Deflection error with element CPS8R.
Mesh
Deflection error (%)
B
C
D
E
1
×
32
0.15
0.05
0.01
−0.01
2
×
32
0.33
0.15
0.09
0.05
4
×
32
0.38
0.18
0.11
0.07
8
×
32
0.40
0.19
0.12
0.08
1
×
64
0.27
0.12
0.07
0.04
2
×
64
0.36
0.17
0.10
0.07
4
×
64
0.39
0.18
0.12
0.08
8
×
64
0.40
0.19
0.12
0.08Table 3
Bending moment error (%).
Element type
Mesh
Computed with stress
Computed with displacement
B
C
D
B
C
D
CPS4
2
×
8
−89.19
−88.99
−88.02
−89.99
−91.72
−95.28
2
×
16
−70.02
−69.89
−69.22
−66.44
−65.26
−57.81
2
×
32
−45.73
−45.67
−45.37
−32.89
−36.96
−59.74
2
×
64
−31.89
−31.87
−31.78
−13.19
−3.40
−620.48
2
×
128
−27.25
−27.24
−27.22
−24.98
−620.48
−3143.38
4
×
64
−14.42
−14.40
−14.27
8
×
64
−10.03
−10.01
−9.87
16
×
64
−8.93
−8.91
−8.77
4
×
128
−8.56
−8.55
−8.52
8
×
128
−3.86
−3.85
−3.81
16
×
128
−2.68
−2.67
−2.64
CPS8R
2
×
32
−33.25
−33.15
−32.59
2
×
64
−33.31
−33.28
−33.11
2
×
128
−33.32
−33.31
−33.24
4
×
64
−6.65
−6.63
−6.52
8
×
64
−1.57
0.02
0.09Bending moment calculated with displacement (with CPS4).
(a)
2
×
8
(b)
2
×
16
(c)
2
×
32
(d)
2
×
64
(e)
2
×
128Figure 5
Bending moment calculated with displacement (2×64, with CPS8R).The errors in the tables are defined as(4a)εω=ωc-ωω×100(%),(4b)εM=Mc-MM×100(%),where εω= relative error of deflection, εM= relative error of bending moment, ωc= calculated deflection using FEM, ω= analytically exact deflection from (1b), Mc= calculated bending moment using FEM, and M= analytically exact bending moment from (1a).From Table1, we can see that the deflection using CPS4 is smaller than the analytical solution, which implies that the shear locking occurs with the first-order element. Shear locking can be easily overcome with second-order reduced-integration elements, that is, using fewer grids can produce approximately the same results (see Table 2).Although using the rectangular first-order element induces shear locking in the beam, we can still obtain a good result with sufficient partitions along the length of the beam; if we partition the beam lengthways in 128 elements, the relative error can be smaller than 3% (see Table1).Table3 and Figure 4 show the computational results of bending moment for the first-order element. They indicate that if the bending moment is to be calculated with stress, sufficient partitions across the cross-section are necessary. Fluctuation can lead to unreliable results when bending moment is calculated with displacement. The effect of variation is even stronger when increasing the grid density. The reason may be that the loss of accuracy occurs for each difference operation and the initial small error will be greatly magnified after two operations.A similar calculation of bending moment was done using CPS8R. Figure5 shows even more notable variation occurring for the bending moments calculated with displacement.


## 2.1. Analytical Solution

The cantilever beam example is shown in Figure1. The width of the square beam is 1 m. The length L is 30 m. A distributed load p=0.5 kPa is applied to the beam. The analytical solution equations are(1a)M=12p(L-x)2,(1b)ω=pL42EI(k22-k33+k412),where M= bending moment, x= the position coordinate, k=x/L, E= the Young’s modulus, I= the moment of inertia, ω=-uy is the deflection of the beam, and uy= displacement in the y direction.Figure 1
Cantilever beam.The beam parameters are taken as Young’s modulusE=2×104 MPa and Poisson’s ratio ν=0.17 in the computation. The element used in the FEM is the 4-node first-order plane stress element (CPS4). The following two methods were used to calculate the bending moment.(a) Calculating Bending Moment with StressThe bending moment was directly computed with the normal stress on the cross-section (see Figure2):
(2)M=∑σiAili,
where σi= normal stress at the centroid of the element, Ai= corresponding area of the element, and li= distance between the centroid and the midline of the beam section.Figure 2
Computing of bending moment with stress.(b) Calculating Bending Moment with DisplacementThe bending moment was calculated using the following quadratic differential of deflection [14]:
(3)M=EId2ωdx2.Equation (3) can be transformed into a difference scheme, and the bending moment was calculated by the difference operation of lateral displacement. Figure 3 is the comparison of the computed bending moment with analytically exact results, where ω is calculated using (1b) with the number of interpolation points being 10, 30, 60, and 120. This shows that the method used here achieved a good result for the bending moment.Figure 3
Calculation versus analytical solution of bending moment.


## 2.2. Errors in Displacement and Bending Moment

We used several meshes to do the computation. Grids along the length of the beam were 32, 64, and 128 (in thex direction, see Figure 1). Sectional partitions were 1, 2, 4, 8, 16, and 32. Two element types, CPS4 and the 8-node reduced-integration element (CPS8R), were employed in the analysis. The results are shown in Tables 1, 2, and 3 and in Figures 4 and 5, where the term m×n denotes that there are m grids along the y direction and n grids along the x direction; B, C, D, and E are positions for error calculation (see Figure 1).Table 1
Deflection error with element CPS4.
Mesh
Deflection error (%)
B
C
D
E
1
×
32
−28.57
−28.69
−28.74
−28.76
4
×
32
−27.20
−27.31
−27.34
−27.36
8
×
32
−27.12
−27.23
−27.27
−27.29
1
6
×
32
−27.10
−27.22
−27.25
−27.27
32
×
32
−27.10
−27.21
−27.25
−27.26
1
×
64
−10.67
−10.83
−10.89
−10.92
4
×
64
−8.45
−8.61
−8.66
−8.68
8
×
64
−8.32
−8.48
−8.54
−8.57
16
×
64
−8.29
−8.45
−8.54
−8.54
32
×
64
−8.28
−8.45
−8.50
−8.53
64
×
64
−8.28
−8.44
−8.50
−8.53
1
×
128
−4.70
−4.87
−4.93
−4.97
4
×
128
−2.12
−2.31
−2.37
−2.40
8
×
128
−1.97
−2.17
−2.23
−2.27
16
×
128
−1.94
−2.13
−2.20
−2.23
32
×
128
−1.93
−2.13
−2.19
−2.23
64
×
128
−1.93
−2.13
−2.19
−2.22
128
×
128
−1.92
−2.13
−2.19
−2.22
2
×
4
−95.94
−95.94
−95.93
−95.92
2
×
8
−85.66
−85.68
−85.68
−85.67
2
×
16
−60.03
−60.08
−60.09
−60.09
2
×
32
−27.50
−27.60
−27.63
−27.65
2
×
64
−8.93
−9.08
−9.12
−9.15
2
×
128
−2.69
−2.86
−2.91
−2.94
2
×
256
−0.99
−1.16
−1.22
−1.25Table 2
Deflection error with element CPS8R.
Mesh
Deflection error (%)
B
C
D
E
1
×
32
0.15
0.05
0.01
−0.01
2
×
32
0.33
0.15
0.09
0.05
4
×
32
0.38
0.18
0.11
0.07
8
×
32
0.40
0.19
0.12
0.08
1
×
64
0.27
0.12
0.07
0.04
2
×
64
0.36
0.17
0.10
0.07
4
×
64
0.39
0.18
0.12
0.08
8
×
64
0.40
0.19
0.12
0.08Table 3
Bending moment error (%).
Element type
Mesh
Computed with stress
Computed with displacement
B
C
D
B
C
D
CPS4
2
×
8
−89.19
−88.99
−88.02
−89.99
−91.72
−95.28
2
×
16
−70.02
−69.89
−69.22
−66.44
−65.26
−57.81
2
×
32
−45.73
−45.67
−45.37
−32.89
−36.96
−59.74
2
×
64
−31.89
−31.87
−31.78
−13.19
−3.40
−620.48
2
×
128
−27.25
−27.24
−27.22
−24.98
−620.48
−3143.38
4
×
64
−14.42
−14.40
−14.27
8
×
64
−10.03
−10.01
−9.87
16
×
64
−8.93
−8.91
−8.77
4
×
128
−8.56
−8.55
−8.52
8
×
128
−3.86
−3.85
−3.81
16
×
128
−2.68
−2.67
−2.64
CPS8R
2
×
32
−33.25
−33.15
−32.59
2
×
64
−33.31
−33.28
−33.11
2
×
128
−33.32
−33.31
−33.24
4
×
64
−6.65
−6.63
−6.52
8
×
64
−1.57
0.02
0.09Bending moment calculated with displacement (with CPS4).
(a)
2
×
8
(b)
2
×
16
(c)
2
×
32
(d)
2
×
64
(e)
2
×
128Figure 5
Bending moment calculated with displacement (2×64, with CPS8R).The errors in the tables are defined as(4a)εω=ωc-ωω×100(%),(4b)εM=Mc-MM×100(%),where εω= relative error of deflection, εM= relative error of bending moment, ωc= calculated deflection using FEM, ω= analytically exact deflection from (1b), Mc= calculated bending moment using FEM, and M= analytically exact bending moment from (1a).From Table1, we can see that the deflection using CPS4 is smaller than the analytical solution, which implies that the shear locking occurs with the first-order element. Shear locking can be easily overcome with second-order reduced-integration elements, that is, using fewer grids can produce approximately the same results (see Table 2).Although using the rectangular first-order element induces shear locking in the beam, we can still obtain a good result with sufficient partitions along the length of the beam; if we partition the beam lengthways in 128 elements, the relative error can be smaller than 3% (see Table1).Table3 and Figure 4 show the computational results of bending moment for the first-order element. They indicate that if the bending moment is to be calculated with stress, sufficient partitions across the cross-section are necessary. Fluctuation can lead to unreliable results when bending moment is calculated with displacement. The effect of variation is even stronger when increasing the grid density. The reason may be that the loss of accuracy occurs for each difference operation and the initial small error will be greatly magnified after two operations.A similar calculation of bending moment was done using CPS8R. Figure5 shows even more notable variation occurring for the bending moments calculated with displacement.


## 3. 2D Analysis of a Sheet Pile Wall

Figure6 shows a sheet pile wall subjected to a load p=1 MPa. The bottom of the domain and the pile tip are fully restrained from moving in any direction while both sides of the domain are restrained in the x direction, while free in the y direction. The length of the pile is the same as that of the above beam, 30 m, and the width is 1 m. We used different mesh partitions to do the calculation with the rectangular first-order element (CPE4).Figure 6
A sheet pile wall subjected to surface load.Since there is no analytical solution for this problem, referring to the above analysis of the cantilever beam, the results using a grid partition of8×60, the 8-node plane strain element, reduced integration (CPE8R) were considered as “exact.” The displacement at the top of the wall, Ut, and the maximum displacement in the middle, Um, are shown in Table 4. The bending moment at the tip of the wall, Mb, and the maximum bending moment in the middle, Mm, are shown in Table 5. Figure 7 shows the distribution of bending moment calculated along the wall. Though a small variation occurs in the bending moment distribution calculated with displacement (CPE8R-8-60-W), it is quite close to that calculated with stress (CPE8R-8-60-Y).Table 4
Horizontal displacement of the wall (with CPE8R).
Item
Horizontal displacement
U
m (m)
0.03647
U
t (m)
−0.00528Table 5
Calculated bending moment in the wall (with CPE8R).
Item
Computed with displacement
Computed with stress
M
b  (kN·m)
2939.947
2650.037
M
m (kN·m)
−906.667
−902.366Figure 7
Bending moment calculated with element CPE8R.The calculated results using CPE4 are presented in Tables6 and 7 and Figures 8, 9, and 10. The relative errors in the tables are relevant to those calculated with CPE8R. We find that the displacements calculated with first-order elements are smaller than those with 8-node elements, reduced-integration, whereas the values are close if they have the same partition (8×60) along the height of the wall, which implies that the shear locking is not distinct as for the cantilever beam. The bending moments calculated with displacement approximate those with stress and have insignificant variation. However, if the wall is suspended in the soil, obvious variation occurs for the element CPE8R as shown in Figure 11.Table 6
Calculated horizontal displacement of the wall (with CPE4).
Item
Grid number along the length
Grid number across wall section
2
4
8
Displacement(m)
Error(%)
Displacement (m)
Error(%)
Displacement(m)
Error(%)
U
m
15
0.02875
−21.16
0.02835
−22.25
0.02825
−22.52
30
0.03475
−4.71
0.03411
−6.47
0.03395
−6.89
60
0.03673
0.74
0.03601
−1.24
0.03584
−1.72
U
t
15
−0.00144
−72.71
−0.00122
−76.86
−0.00117
−77.88
30
−0.00450
−14.71
−0.00420
−20.48
−0.00412
−21.90
60
−0.00541
2.49
−0.00509
−3.63
−0.00501
−5.14Table 7
Calculated bending moment in the wall (with CPE4).
Item
Grid number alongthe length
Grid number across wall section
2
4
8
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
M
b
Computed with stress
15
1096.813
−58.61
1335.297
−49.61
1393.136
−47.43
30
1714.913
−35.29
2052.434
−22.55
2133.281
−19.50
60
2073.363
−21.76
2454.809
−7.37
2546.233
−3.92
Computed with displacement
15
1073.083
−63.50
1048.292
−64.34
1042.292
−64.55
30
1966.567
−33.11
1880.900
−36.02
1860.867
−36.70
60
2724.200
−7.34
2521.400
−14.24
2478.267
−15.70
M
m
Computed with stress
15
−535.650
−40.64
−658.054
−27.07
−687.977
−23.76
30
−657.713
−27.11
−805.006
−10.79
−840.821
−6.82
60
−694.513
−23.03
−850.134
−5.79
−887.903
−1.60
Computed with displacement
15
−690.833
−23.81
−679.167
−25.09
−676.667
−25.37
30
−851.667
−6.07
−840.000
−7.35
−836.667
−7.72
60
−900.000
−0.74
−886.667
−2.21
−886.667
−2.21Comparison of horizontal displacements.
(a)
Grid number being 15 along the wall
(b)
Grid number being 30 along the wall
(c)
Grid number being 60 along the wallComparison of bending moments calculated with displacement.
(a)
Grid number being 15 along the wall
(b)
Grid number being 30 along the wall
(c)
Grid number being 60 along the wallComparison of bending moments calculated with stress.
(a)
Grid number being 15 along the wall
(b)
Grid number being 30 along the wall
(c)
Grid number being 60 along the wallBending moments calculated with displacement for suspended sheet pile wall.
(a)
Suspended sheet pile wall
(b)
Bending moments along the wall


## 4. Pile Examples

### 4.1. 3D Analysis of a Pile

Figure12 shows a pile and its surrounding soil in 3D view. The length of the pile is 30 m, and the width of the square pile is 1 m, which is also the same as that of the above-mentioned beam. The width of the computational domain is 9 m. The applied load, material properties, and boundary conditions are the same as those of the above-mentioned sheet pile wall.Figure 12
3-D finite element mesh.3-D analyses were made with different mesh partitions of the pile shaft. The results with a grid partition of8×60, a 20-node brick element and reduced integration (C3D20R) were considered as “exact”. Displacement at the top of the pile, Ut, and maximum displacement in the middle, Um, are shown in Table 8, while the bending moment at the tip of the pile, Mb, and maximum bending moment in the middle, Mm, are shown in Table 9. Figure 13 is the distribution of moment calculated along the pile. Again, the moment calculated with displacement (C3D20R-8-60-W) is close to that calculated with stress (C3D20R-8-60-Y) and insignificant variation is found.Table 8
Horizontal displacement of the pile (with C3D20R).
Item
Horizontal displacement (m)
U
m
0.045113
U
t
−0.007570Table 9
Calculated bending moment in the pile (with C3D20R).
Item
Computed with displacement
Computed with stress
M
b (kN·m)
3967.200
4706.848
M
m (kN·m)
−1146.667
−1101.314Figure 13
Bending moment calculated with element C3D20R.The calculation results with first-order element (C3D8) are presented in Tables10 and 11 and Figures 14, 15, and 16. The relative errors in the tables are relevant to those calculated with C3D20R. As before, displacements calculated with first-order elements are smaller than those with C3D20R and the values are close if they have the same partition (8×60) along the length of the pile. Bending moments calculated with displacement approximate those with stress and have insignificant variation. This indicates that the response of the pile is similar to that of the above-mentioned wall for displacement and bending moment computations.Table 10
Calculated horizontal displacement of the pile (with C3D8).
Item
Grid number along the length
Grid number across pile section
2
4
8
Displacement (m)
Error(%)
Displacement(m)
Error(%)
Displacement(m)
Error(%)
U
m
15
0.042865
−4.98
0.042764
−5.21
0.042740
−5.26
30
0.044506
−1.35
0.044408
−1.56
0.044384
−1.62
60
0.044896
−0.48
0.044796
−0.70
0.044774
−0.75
U
t
15
−0.006532
−13.71
−0.006488
−14.29
−0.006478
−14.43
30
−0.007278
−3.85
−0.007238
−4.38
−0.007229
−4.51
60
−0.007452
−1.56
−0.007414
−2.06
−0.007405
−2.17Table 11
Calculated bending moment in the pile (with C3D8).
Item
Grid number along the length
Grid number across pile section
2
4
8
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
M
b
Computed with stress
15
2184.875
−53.58
2668.929
−43.30
2787.809
−40.77
30
3172.975
−32.59
3831.406
−18.60
3992.286
−15.18
60
3793.750
−19.40
4545.541
−3.43
4729.851
0.49
Computed with displacement
15
1881.417
−52.58
1863.667
−53.02
1859.417
−53.13
30
3248.333
−18.12
3179.000
−19.87
3163.000
−20.27
60
4082.267
2.90
3933.733
−0.84
3895.867
−1.80
M
m
Computed with stress
15
−816.513
−25.86
−1004.993
−8.75
−1051.735
−4.50
30
−847.748
−23.02
−1041.848
−5.40
−1089.858
−1.04
60
−858.825
−22.02
−1053.395
−4.35
−1102.111
0.07
Computed with displacement
15
−1063.333
−7.27
−1060.000
−7.56
−1059.583
−7.59
30
−1108.333
−3.34
−1105.000
−3.63
−1103.333
−3.78
60
−1126.667
−1.74
−1126.667
−1.74
−1120.000
−2.33Comparison of horizontal displacements.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaftComparison of bending moments calculated with displacement.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaftComparison of bending moments calculated with stress.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaft


### 4.2. 3-D Analysis of Piles for a Bridge Abutment

Figure17 shows the cross-section of the piles of a bridge and a polder dike. Figure 18 shows one of the meshes for the computation. The length of the computational domain is 700 m and the width is 60 m, which equals the abutment span. The elevation of the ground surface is −5.3 m, the elevation of the bottom of the domain is −120 m, and the elevation of the pile tip is −90 m. Each round pile is represented by an equivalent square pile with a width of 1.33 m, and the four piles are connected by a pile cap (see Figure 18(b)). The flyash is filled to 4.63 m (see Figure 17).Figure 17
Cross-section of the piles of a bridge and a polder dike (in m).3-D finite element mesh for the piles.
(a)
Computation domain
(b)
Meshes for pilesThe weight of the dike and the flyash was simulated with a distributed load acting on the ground surface, respectively. The effective unit weight of the dike is 18 kN/m3 above the water level and 11 kN/m3 below the water level. The effective unit weight of the flyash is 13.5 kN/m3 and 5.9 kN/m3, respectively. Again, linear elastic model was used to simulate the soil and the pile. The parameters of the soil strata are presented in Table 12, while Young’s modulus and Poisson’s ratio for the pile are the same as those of the above-mentioned sheet pile wall.Table 12
Soil layers and parameters.
No.
Soil layer
Elevation of layer top (m)
Elevation of layer bottom (m)
Poisson’s ratio\nu
Young’s modulus (MPa)
1
Muck
−5.3
−7.2
0.35
1.0
2
Mucky soil
−7.2
−24.2
0.40
1.0
3
Silty fine sand
−24.2
−36.3
0.25
10
4
Clay
−36.3
−54.1
0.33
10
5
Clay
−54.1
−56.6
0.33
25
6
Medium sand
−56.6
−62.2
0.25
30
7
Clay
−62.2
−71.4
0.3
30
8
Clay
−71.4
−77.2
0.3
30
9
Clay
−77.2
−86.1
0.3
30
10
Clay
−86.1
−120.0
0.3
50Five meshes (M1 to M5) were used to make the computation. The partitions of the piles and the total number of elements and nodes in each mesh are listed in Table13. Mesh M3 has the most number of elements and nodes and the most grids (4×42) for the pile. Again, the results for M3, a 20-node brick element with reduced integration (C3D20R-M3), were considered as “exact” to calculate the relative errors.Table 13
Meshes for the analysis.
Item
C3D20R-M3
C3D8-M1
C3D8-M2
C3D8-M3
C3D8-M4
C3D8-M5
Element number
22620
9541
15855
22620
13497
14283
Node number
99996
11208
18264
25824
15600
16488
Grids for pile
4 × 42
4
×
16
4
×
26
4 × 42
1
×
26
2
×
26
CPU-time (s)
1451.8
67.1
132.6
217.5
124.0
121.0The calculated displacements with different meshes and element types are presented in Table14 and the distributions are shown in Figure 19. The results agree closely with each other, that is, with fewer linear elements, one can achieve satisfactory displacement results.Table 14
Calculated horizontal displacement of pileA.
Item
At elevation of −5.3 m
At elevation of −90 m
C3D20R-M3
Displacement (m)
0.02024
0.00642
C3D8-M1
Displacement (m)Error (%)
0.020093−0.71
0.006339−1.22
C3D8-M2
Displacement (m)Error (%)
0.019967−1.33
0.006355−0.96
C3D8-M3
Displacement (m)Error (%)
0.019897−1.68
0.006368−0.77
C3D8-M4
Displacement (m)Error (%)
0.019971−1.31
0.006355−0.96
C3D8-M5
Displacement (m)Error (%)
0.019967−1.33
0.006355−0.97Figure 19
Comparison of horizontal displacements.Figure20 shows the distribution of bending moment with C3D20R-M3. As a whole, the moment calculated with displacement (C3D20R-M3-W) is similar to that calculated with stress (C3D20R-M3-Y). However, notable variation occurs for the moment calculated with displacement, particularly at the upper part of the pile.Figure 20
Bending moment calculated with element C3D20R.Figure21 shows the bending moment calculated with stress obtained from element type C3D8 and meshes M1, M2, M3, and M5. Mesh M4 has no partition across the pile section so we cannot calculate the bending moment. It is obvious that the results for M3 and M2 are closer to that of C3D20R-M3 than for M1 and M5. Mesh M5 has only two grids across the pile section, and the moment calculated for M5 is unreliable and quite different from other results.Figure 21
Bending moments calculated with stress.Figure22 shows the bending moment calculated with displacement. It is found that unlike C3D20R, the distribution has insignificant variation for the linear element (C3D8) with different meshes and agrees closely with each other, which implies that calculations with linear elements may produce fewer and smaller fluctuations than high-order elements in this case.Bending moments calculated with displacement.
(a)
With mesh M1, M2, and M3
(b)
With mesh M4 and M5


### 4.3. 2-D Analysis of a Pile Row

Generally, a pile row can be replaced by a sheet pile wall with stiffness chosen as the average of the pile stiffness and that of the soil between the piles [11–13],
(5)EI=EpIp+EsIs,
where E= equivalent modulus of the sheet pile wall, I= moment of inertia of the sheet pile wall, Ep, Es= Young’s moduli of the pile and the soil, respectively, and Ip, Is= moments of inertia of the pile and the soil, respectively.If the piles are at a spacing ofu and each pile is squared with a width of d, the equivalent modulus can be
(6)E=Epd+Es(u-d)u.In the analysis, the spacing of the pile row was assumed to beu=2 m, 9 m, 100 m, and 1000 m. Other parameters are the same as those of the above analyses. The configurations are shown in Figures 6 and 12. As in the preceding studies, the grid partition for the pile shaft was 8×60. The computational results are shown in Figure 23 and Tables 15 and 16, in which “-D” denotes that the 2-D analysis was conducted with the equivalent modulus of the piles and (×Ep/E) denotes that the bending moment was modified by multiplying by Ep/E. It is found that the calculated displacements of the “equivalent sheet pile wall” are in close agreement with that of the pile row. However, the results for the bending moment show some difference, particularly at the fixed tip of the pile.Table 15
Comparison of horizontal displacements calculated by 2-D and 3-D analyses.
Pile spacing
2 m
9 m
100 m
1000 m
Element type
CPE4-D
C3D8
CPE4-D
C3D8
CPE4-D
C3D8
CPE4-D
C3D8
U
m (m)
0.04008
0.04016
0.04391
0.04477
0.04543
0.04575
0.04608
0.04575
U
t (m)
−0.00675
−0.00672
−0.00776
−0.007405
−0.00782
−0.00732
−0.00766
−0.00732Table 16
Comparison of bending moments by 2-D and 3-D analyses (u=9 m).
Element type
CPE4-D-W(EI)
CPE4-D-W(EpIp)
C3D8-W
CPE4-D-Y
CPE4-D-Y(×Ep/E)
C3D8-Y
M
b (kN·m)
4422.84
4387.73
3895.87
636.48
5682.86
4729.85
M
m (kN·m)
−1061.76
−1053.33
−1120.00
−117.02
−1044.83
−1102.11Comparison of computational results with 2-D and 3-D analyses.
(a)
Horizontal displacement
(b)
Bending moment


## 4.1. 3D Analysis of a Pile

Figure12 shows a pile and its surrounding soil in 3D view. The length of the pile is 30 m, and the width of the square pile is 1 m, which is also the same as that of the above-mentioned beam. The width of the computational domain is 9 m. The applied load, material properties, and boundary conditions are the same as those of the above-mentioned sheet pile wall.Figure 12
3-D finite element mesh.3-D analyses were made with different mesh partitions of the pile shaft. The results with a grid partition of8×60, a 20-node brick element and reduced integration (C3D20R) were considered as “exact”. Displacement at the top of the pile, Ut, and maximum displacement in the middle, Um, are shown in Table 8, while the bending moment at the tip of the pile, Mb, and maximum bending moment in the middle, Mm, are shown in Table 9. Figure 13 is the distribution of moment calculated along the pile. Again, the moment calculated with displacement (C3D20R-8-60-W) is close to that calculated with stress (C3D20R-8-60-Y) and insignificant variation is found.Table 8
Horizontal displacement of the pile (with C3D20R).
Item
Horizontal displacement (m)
U
m
0.045113
U
t
−0.007570Table 9
Calculated bending moment in the pile (with C3D20R).
Item
Computed with displacement
Computed with stress
M
b (kN·m)
3967.200
4706.848
M
m (kN·m)
−1146.667
−1101.314Figure 13
Bending moment calculated with element C3D20R.The calculation results with first-order element (C3D8) are presented in Tables10 and 11 and Figures 14, 15, and 16. The relative errors in the tables are relevant to those calculated with C3D20R. As before, displacements calculated with first-order elements are smaller than those with C3D20R and the values are close if they have the same partition (8×60) along the length of the pile. Bending moments calculated with displacement approximate those with stress and have insignificant variation. This indicates that the response of the pile is similar to that of the above-mentioned wall for displacement and bending moment computations.Table 10
Calculated horizontal displacement of the pile (with C3D8).
Item
Grid number along the length
Grid number across pile section
2
4
8
Displacement (m)
Error(%)
Displacement(m)
Error(%)
Displacement(m)
Error(%)
U
m
15
0.042865
−4.98
0.042764
−5.21
0.042740
−5.26
30
0.044506
−1.35
0.044408
−1.56
0.044384
−1.62
60
0.044896
−0.48
0.044796
−0.70
0.044774
−0.75
U
t
15
−0.006532
−13.71
−0.006488
−14.29
−0.006478
−14.43
30
−0.007278
−3.85
−0.007238
−4.38
−0.007229
−4.51
60
−0.007452
−1.56
−0.007414
−2.06
−0.007405
−2.17Table 11
Calculated bending moment in the pile (with C3D8).
Item
Grid number along the length
Grid number across pile section
2
4
8
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
Bending moment(kN·m)
Error(%)
M
b
Computed with stress
15
2184.875
−53.58
2668.929
−43.30
2787.809
−40.77
30
3172.975
−32.59
3831.406
−18.60
3992.286
−15.18
60
3793.750
−19.40
4545.541
−3.43
4729.851
0.49
Computed with displacement
15
1881.417
−52.58
1863.667
−53.02
1859.417
−53.13
30
3248.333
−18.12
3179.000
−19.87
3163.000
−20.27
60
4082.267
2.90
3933.733
−0.84
3895.867
−1.80
M
m
Computed with stress
15
−816.513
−25.86
−1004.993
−8.75
−1051.735
−4.50
30
−847.748
−23.02
−1041.848
−5.40
−1089.858
−1.04
60
−858.825
−22.02
−1053.395
−4.35
−1102.111
0.07
Computed with displacement
15
−1063.333
−7.27
−1060.000
−7.56
−1059.583
−7.59
30
−1108.333
−3.34
−1105.000
−3.63
−1103.333
−3.78
60
−1126.667
−1.74
−1126.667
−1.74
−1120.000
−2.33Comparison of horizontal displacements.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaftComparison of bending moments calculated with displacement.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaftComparison of bending moments calculated with stress.
(a)
Grid number being 15 along the pile shaft
(b)
Grid number being 30 along the pile shaft
(c)
Grid number being 60 along the pile shaft


## 4.2. 3-D Analysis of Piles for a Bridge Abutment

Figure17 shows the cross-section of the piles of a bridge and a polder dike. Figure 18 shows one of the meshes for the computation. The length of the computational domain is 700 m and the width is 60 m, which equals the abutment span. The elevation of the ground surface is −5.3 m, the elevation of the bottom of the domain is −120 m, and the elevation of the pile tip is −90 m. Each round pile is represented by an equivalent square pile with a width of 1.33 m, and the four piles are connected by a pile cap (see Figure 18(b)). The flyash is filled to 4.63 m (see Figure 17).Figure 17
Cross-section of the piles of a bridge and a polder dike (in m).3-D finite element mesh for the piles.
(a)
Computation domain
(b)
Meshes for pilesThe weight of the dike and the flyash was simulated with a distributed load acting on the ground surface, respectively. The effective unit weight of the dike is 18 kN/m3 above the water level and 11 kN/m3 below the water level. The effective unit weight of the flyash is 13.5 kN/m3 and 5.9 kN/m3, respectively. Again, linear elastic model was used to simulate the soil and the pile. The parameters of the soil strata are presented in Table 12, while Young’s modulus and Poisson’s ratio for the pile are the same as those of the above-mentioned sheet pile wall.Table 12
Soil layers and parameters.
No.
Soil layer
Elevation of layer top (m)
Elevation of layer bottom (m)
Poisson’s ratio\nu
Young’s modulus (MPa)
1
Muck
−5.3
−7.2
0.35
1.0
2
Mucky soil
−7.2
−24.2
0.40
1.0
3
Silty fine sand
−24.2
−36.3
0.25
10
4
Clay
−36.3
−54.1
0.33
10
5
Clay
−54.1
−56.6
0.33
25
6
Medium sand
−56.6
−62.2
0.25
30
7
Clay
−62.2
−71.4
0.3
30
8
Clay
−71.4
−77.2
0.3
30
9
Clay
−77.2
−86.1
0.3
30
10
Clay
−86.1
−120.0
0.3
50Five meshes (M1 to M5) were used to make the computation. The partitions of the piles and the total number of elements and nodes in each mesh are listed in Table13. Mesh M3 has the most number of elements and nodes and the most grids (4×42) for the pile. Again, the results for M3, a 20-node brick element with reduced integration (C3D20R-M3), were considered as “exact” to calculate the relative errors.Table 13
Meshes for the analysis.
Item
C3D20R-M3
C3D8-M1
C3D8-M2
C3D8-M3
C3D8-M4
C3D8-M5
Element number
22620
9541
15855
22620
13497
14283
Node number
99996
11208
18264
25824
15600
16488
Grids for pile
4 × 42
4
×
16
4
×
26
4 × 42
1
×
26
2
×
26
CPU-time (s)
1451.8
67.1
132.6
217.5
124.0
121.0The calculated displacements with different meshes and element types are presented in Table14 and the distributions are shown in Figure 19. The results agree closely with each other, that is, with fewer linear elements, one can achieve satisfactory displacement results.Table 14
Calculated horizontal displacement of pileA.
Item
At elevation of −5.3 m
At elevation of −90 m
C3D20R-M3
Displacement (m)
0.02024
0.00642
C3D8-M1
Displacement (m)Error (%)
0.020093−0.71
0.006339−1.22
C3D8-M2
Displacement (m)Error (%)
0.019967−1.33
0.006355−0.96
C3D8-M3
Displacement (m)Error (%)
0.019897−1.68
0.006368−0.77
C3D8-M4
Displacement (m)Error (%)
0.019971−1.31
0.006355−0.96
C3D8-M5
Displacement (m)Error (%)
0.019967−1.33
0.006355−0.97Figure 19
Comparison of horizontal displacements.Figure20 shows the distribution of bending moment with C3D20R-M3. As a whole, the moment calculated with displacement (C3D20R-M3-W) is similar to that calculated with stress (C3D20R-M3-Y). However, notable variation occurs for the moment calculated with displacement, particularly at the upper part of the pile.Figure 20
Bending moment calculated with element C3D20R.Figure21 shows the bending moment calculated with stress obtained from element type C3D8 and meshes M1, M2, M3, and M5. Mesh M4 has no partition across the pile section so we cannot calculate the bending moment. It is obvious that the results for M3 and M2 are closer to that of C3D20R-M3 than for M1 and M5. Mesh M5 has only two grids across the pile section, and the moment calculated for M5 is unreliable and quite different from other results.Figure 21
Bending moments calculated with stress.Figure22 shows the bending moment calculated with displacement. It is found that unlike C3D20R, the distribution has insignificant variation for the linear element (C3D8) with different meshes and agrees closely with each other, which implies that calculations with linear elements may produce fewer and smaller fluctuations than high-order elements in this case.Bending moments calculated with displacement.
(a)
With mesh M1, M2, and M3
(b)
With mesh M4 and M5


## 4.3. 2-D Analysis of a Pile Row

Generally, a pile row can be replaced by a sheet pile wall with stiffness chosen as the average of the pile stiffness and that of the soil between the piles [11–13],
(5)EI=EpIp+EsIs,
where E= equivalent modulus of the sheet pile wall, I= moment of inertia of the sheet pile wall, Ep, Es= Young’s moduli of the pile and the soil, respectively, and Ip, Is= moments of inertia of the pile and the soil, respectively.If the piles are at a spacing ofu and each pile is squared with a width of d, the equivalent modulus can be
(6)E=Epd+Es(u-d)u.In the analysis, the spacing of the pile row was assumed to beu=2 m, 9 m, 100 m, and 1000 m. Other parameters are the same as those of the above analyses. The configurations are shown in Figures 6 and 12. As in the preceding studies, the grid partition for the pile shaft was 8×60. The computational results are shown in Figure 23 and Tables 15 and 16, in which “-D” denotes that the 2-D analysis was conducted with the equivalent modulus of the piles and (×Ep/E) denotes that the bending moment was modified by multiplying by Ep/E. It is found that the calculated displacements of the “equivalent sheet pile wall” are in close agreement with that of the pile row. However, the results for the bending moment show some difference, particularly at the fixed tip of the pile.Table 15
Comparison of horizontal displacements calculated by 2-D and 3-D analyses.
Pile spacing
2 m
9 m
100 m
1000 m
Element type
CPE4-D
C3D8
CPE4-D
C3D8
CPE4-D
C3D8
CPE4-D
C3D8
U
m (m)
0.04008
0.04016
0.04391
0.04477
0.04543
0.04575
0.04608
0.04575
U
t (m)
−0.00675
−0.00672
−0.00776
−0.007405
−0.00782
−0.00732
−0.00766
−0.00732Table 16
Comparison of bending moments by 2-D and 3-D analyses (u=9 m).
Element type
CPE4-D-W(EI)
CPE4-D-W(EpIp)
C3D8-W
CPE4-D-Y
CPE4-D-Y(×Ep/E)
C3D8-Y
M
b (kN·m)
4422.84
4387.73
3895.87
636.48
5682.86
4729.85
M
m (kN·m)
−1061.76
−1053.33
−1120.00
−117.02
−1044.83
−1102.11Comparison of computational results with 2-D and 3-D analyses.
(a)
Horizontal displacement
(b)
Bending moment


## 5. Conclusions

The bending moment computational methods for piles were investigated using a series of calculation examples in this study, and the following conclusions were reached.(1)
Compared to a cantilever beam, shear locking is not significant for the passive pile embedded in soil, so higher-order elements are not always necessary for the computation. Computation with first-order (linear) elements and appropriate grid partition can produce similar good results as for higher-order elements.(2)
The number of the grids along the length of the pile plays an important role in the analysis. With an increase in grid number, the calculated displacement and bending moment are closer to theoretical results. Increasing the grid number across the pile section is helpful for increasing the accuracy of the bending moment calculated with stress, while it has insignificant influence on displacement and the related bending moment calculation.(3)
Calculating bending moment with stress can produce good results, but many grids are needed to partition the pile section. Calculating bending moment with displacement needs fewer grids across the pile section, but it may result in fluctuations of the results, especially for the cantilever beam presented in Section2. The reason may be that the bending moment calculated with stress corresponds to the “integration” operation of stress, while the bending moment calculated with displacement corresponds to the “difference” operation of displacement. The difference operation may amplify the error, and the initial small error will be greatly magnified after two operations. Consequently, if the fluctuations of bending moment calculated with displacement are evident, it is suggested that the bending moment should be calculated with stress.(4)
When calculating the displacements of the piles, a pile row can be suitably represented by an equivalent sheet pile wall which has the same flexural stiffness per unit width as the piles and the soil it replaces. The displacements of the wall can agree closely with that of the pile row, while bending moments may differ from each other.(5)
A special attention should be given to meshing and the computational method for bending moment. Calculated results may differ greatly with different grid partitions and computational methods. Comparison of results using different meshes is necessary when performing the analysis.It should be noted that, in order to clearly reveal the influences of element type and mesh partition, only linear elastic model was used in this study to simulate the soil and the pile. Obviously, introduction of constitutive models for the analysis of actual soil and pile will further complicate the problem, and thus, more attention should be paid to the calculation methods of bending moment.



---

*Source: 784583-2013-06-17.xml*