1966). The third again used a cubic fit over stress-strain data, but this time the data were obtained from an unconfined test. The most realistic field situation was obtained from the third constitutive relationship which allowed the most deformation for a given load. The other two constitutive relationships restricted the values of strain to less than 20% which is low for agricultural purposes. These two relationships also predicted that the elastic strain would be much greater than the

more »

... strain which is unreasonable for compactable agricultural soil. Clearly more work needs to be done before the finite element method can be used to accurately predict soil compaction. With recent advances at the National Soil Dynamics Laboratory and Auburn University in development of a constitutive relationship for agricultural soils, portions of the necessary technology have now been developed to make the finite element method more attractive. Therefore, the objectives of this research were to: 1. Construct a finite element computer program to take into effect axisymmetric geometry and linear elastic constitutive relationships. Verify the linear elastic finite element program by comparing results against the linear elastic theory developed by Bousinnesq. A slightly different method of storing the stiffness matrix was also employed by this program. Most stiffness matrices for finite element problems are banded, symmetric, and positive definite. Utilizing these characteristics allowed all the necessary variables to be stored in a matrix with the same length but with the width equal to the half-bandwidth. For our example (which is discussed below) the storage space was decreased from 1229K to only llOK for double precision arithmetic using this storage technique. Program verification The program was verified using a static flexible circular footing. This is a common technique utilized in several civil engineering applications mentioned earlier. The parameters necessary to make this program work were the finite element mesh (Figure 5), loading type and amount, and soil properties. This mesh was intended to be large enough so that the total radial distance was 6 times greater than the radius of the pressure load. The factor '6' was thought to be sufficient to introduce a semi-infinite soil medium that allowed large stresses and strains near the central axis without being affected by fixed boundary conditions. The mesh was composed of 169 elements and 196 nodes and extended radially from the centerline 1.44 m and vertically into the soil 1.44 m. The loading was distributed evenly over the centermost four 24 elements at the top of the mesh. The load corresponded to a similar problem solved in Girijavallabhan and Reese (1968). Values input into their program were: Pressure Load = 4.788 kPa Young's modulus = 4788 kPa Poisson's ratio = 0.3 Radius of pressure load = 1.524 m Their results are shown in Figure 6. They evaluated their program by plotting the surface settlement against the theoretical linear elastic equations developed by Boussinesq and against other equations assuming a rigid boundary developed by Steinbrenner (Terzaghi, 1943). Note how the finite element method underpredicted the surface settlement compared to the Boussinesq solution. This was undoubtedly due to the assumption of an infinite medium that the theoretical solution made. The presence of a boundary that was needed by the finite element method probably would result in less deformation of the surface. The technique we used to compare the current finite element program and the Boussinesq solution was similar to the one just described. The same soil strength parameters were used as in the above example. The loading and radius of loading were changed to be more similar to a tire (Pollack et al., 1985). These included loading 25 to 125 kPa over a radius of 24 cm. The soil surface displacements obtained from the finite element method compared very favorably against the theoretical values obtained from linear elastic theory (Figure 7). Results are shown both for a model that starts at the central axis and for a model that is offset from this axis by the width of one element (0.06 m). The accurate prediction of surface deformation by both of these models indicated that the program was working correctly and accurately predicted surface deformation. The prediction of vertical stress and volumetric strain (Figures 8 and 9) was also very similar to the Bousinnesq solutions (Figures 10 and 11). Note the location and value of the maximum values and the general shape of the pressure and strain bulbs. Upon further investigation, however, a problem was noted. When the model was subjected to only gravitational loading, the stress and strain contours were not straight across as they should have been, but tended to deviate toward the center of the mesh both for vertical stress and volumetric strain (Figures 12 and 13). These errors indicated that perhaps the model was not working correctly. Perhaps the averaging method was giving faulty results near the center of the finite element mesh. Because this area was of primary interest, it was decided that a new type of element would be incorporated into the model that would allow numerical integration much easier. [K] = 2 a, . G. . 2 n. X. (45) ij J J k * * where G^j is the matrix EG] = CBJ^CCICB] det CJ] evaluated at the Gauss-Legendre sampling points r^ and sj, and are the given Gauss-Legendre weighting factors. Two subroutines which used the above numerical procedure for calculating the stiffness matrix were obtained from Bathe (1982). Node 4 Node 1 Node Node 3 Node 2 FIGURE 4. Method of using 4 triangular elements to form 1 quadrilateral element 36 FIGURE 5, Axiayrametric finite-element mesh that is used to model a flexible circular footing 37 RADIUS = 3 «#. COMPUTER SOLUTION _ (RIOID BOUNDARY, DA : 9) z H Z 111 z bJ STEINBRENNER'S SOLUTION (RIGID BOUNDARY, D/rs 3) 005 ui in BOUSSINESO'S SOLUTION ( SEMI-INFINITE , D/r » « ) ui u < I u. K 3 M 010 FIGURE 6. Settlement of the surface when a circular flexible footing is applied (After Girijavallabhan and Reese, 1968) Volumetric Strain -1.20 -0.90 -0.60 -0.30 -0.00 0.30 0.60 0.90 1.20 -0.30 0.30 -0.60 -0.60 0.90 -0.90 ^'^-1.20 -0.90 -1.20 0.60 -0.30 -0.00 0.30 0.60 0.90 1.20 FIGURE 9. Contours of volumetric strain obtained from finite element solution of flexible circular footing problem using triangular elements % 1 1 1 1 N. 1 / % \ / 1 X-/ 1 \ » 1 \ / / / RESULTS The entire data set obtained from the laboratory experiment is included in Appendix F. A portion of the results from the laboratory experiment carried out in the soil bin at the NSDL are shown in Table 4. This table shows the displacements of the soil surface when the load was applied and the displacements and final locations of the stress state transducers. Six values were averaged for the surface displacements and three values were averaged for the stress state transducers. The displacements of the flat disc were uniform over the surface of the disc, but the displacements given for the spherical disc were maximum in the center and decreased near the edges due to its geometry (Figure 38). The values given for the surface displacements and the shape of the loading device were used to load the finite element mesh. They were split into multiple load steps. Values of stresses and strains were obtained at the centroid of each element for each respective load step. The stresses obtained at the element centroids were compared to the values obtained from the soil bin experiment. The vertical stresses obtained from the soil bin experiment were the largest and probably the most reliable. The resolution of the stress state transducers was approximately ± 5 kPa. The largest values that were recorded for the horizontal stresses were only about 50 kPa. The largest values for the vertical stresses were about 300 kPa. The resolution of the SSTs was much less important in the