• Federal Highway Administration >
• Publications >
• Research Publications >
• 04094 >
• Alttext.Cfm
• Evaluation of LS-DYNA Soil Material Model 147

Evaluation of LS-DYNA Soil Material Model 147

PDF files can be viewed with the Acrobat® Reader®

Alternative Text for Images

Figure 1. Large-scale direct shear testing device. Photo.

The picture shows a large-scale device for direct shear testing in an indoor laboratory setting. The 500-millimeter red tester is bolted to the floor with four large bolts and nuts. Force is applied pneumatically from a pressure source in the foreground. The cylinder holds a square block with four collars around it that represent the 18.5-kilopascal overburden. This device replicates at a scale 10 times larger than American Society for Testing Materials (ASTM) specifications for direct shear testing, which is needed to accurately capture shear behavior of National Cooperative Highway Research Program Report 350 strong soil, the standard for roadside safety crash testing.

Return to Figure 1

Figure 2. General direct shear testing performance. Graph.

The graph shows a steeply rising curve that depicts direct shear performance during testing. The horizontal axis is displacement or strain, while the vertical axis is resisting shear stress with residual shear stress at the midpoint and peak shear stress at the top of the graph. The curve rises sharply to the peak shear stress until leveling off at the midpoint residual shear stress, which indicates a break from soil continuity when shear force causes fracture. The level line represents sliding over the newly formed soil interface caused by frictional force.

Return to Figure 2

Figure 3. Direct shear testing results (performed at the user’s facility) Graph.

The graph shows the results of two cylinder halves with continuous soil volume moved diametrically in a direct shear test. At first, the force increases steadily until the peak, then begins to fracture and sustain damage that accounts for the sliding with friction and the drop in the graph. The two shear tests plot approximately the same. The horizontal axis is displacement in millimeters ranging from 0 to 250. The vertical axis is force in kilonewtons ranging from 0 to 10. Each line starts at zero, rises dramatically to 8.5 kilonewtons at 40 millimeters of displacement and then falls gradually to 0.5 kilonewtons at 250 millimeters. According to the developer, the soil model is only designed to handle the simulation up to the point of material separation, not the majority of the force-displacement history shown in this figure.

Return to Figure 3

Figure 4. Finite element model of large-scale direct shear test. Drawing.

A three-dimensional computer model grid of the direct shear tests with 10,600 elements was created to measure the effects of the various soil parameters. The model included solid rigid elements to simulate both the steel and the type 147 soil. The drawing shows the two-part steel cylinder with the bottom half yellow (where the motion condition of 1 millimeter per millisecond was applied) and the top half green, filled with a blue cylinder of soil. The cylinder sits on a brown input deck and is surrounded by the mechanical components of the raspberry-colored testing device.

Figure 5. Soil displacement: Baseline model. Drawing.

The two-part computer model shows the test cylinders unchanged and plumb (green stacked exactly on blue) at (A), which represents 0 milliseconds. At (B), or 49.992 milliseconds, the top green cylinder is deflected to the left about 50 millimeters. In the left drawing, each surface cylinder of soil begins in contact with the other and at the end of the time analysis, they are contacting metal and air. The area between the two cylinder halves becomes a noncontinuous slide surface.

Return to Figure 5

Figure 6. Direct shear force: Baseline model. Graph.

The graph predicts the initial increase of driving force to peak and the decrease as soil sustains damage. The horizontal axis is displacement in millimeters ranging from 0 to 50. The vertical axis is X-force (kilonewtons times internal energy) ranging from 0 to 1.4. Section IDI or A1 starts at zero and rises steadily, except for one dip at 25 millimeters and 0.4 X-force, to end at 1.3 X-force at 50 millimeters. According to the developer, the soil modeling appears to accurately correlate with the actual tests in figure 3.

Return to Figure 6

Figure 7. Direct shear soil internal energy: Baseline model. Graph.

In this graph, the horizontal axis is displacement in millimeters ranging from 0 to 50. The vertical axis is internal energy (E plus 3) ranging from 0 to 60. Component A for internal energy shows a concave curve that rises gradually from zero to 20 internal energy at 30 millimeters displacement to 68 internal energy at 50 millimeters displacement.

Return to Figure 7

Figure 8. Graphical representation of the Mohr-Coulomb failure criteria. Graph.

The Mohr theory states that materials rupture from a critical combination of normal and shearing stress, not from a maximum of either alone. The graphic depicts the normal stress on a horizontal line. Two half circles sit on this line, and the points where they touch are labeled normal stress subscripts 3, 1, 3, and 1. The shear stress forms a vertical line at 90 degrees from the normal stress. A skewed line starting at the left of the normal stress line rises over the half circles as a hypotenuse to form an acute triangle. To the right, the angle of internal friction is a line off this hypotenuse that is parallel to the normal stress line. On the left, C for cohesion (the bond between soil particles) is the distance from the normal stress line to where the vertical shear stress line intersects the hypotenuse of the triangle.

Return to Figure 8

Figure 9. Cohesion parameter (COH) variations. Graphs.

Both these shear test simulations of displacement (horizontal axis ranging from 0 to 30) versus X-force (vertical axis ranging form 200 to 1000) or internal energy (vertical axis ranging from 0 to 14) show little variation. For the X-force graph, Line A (red) equals 0.0, Line B (green) equals 6.2E minus 8, Line C (dark blue) equals 6.2E minus 6 and Line D (turquoise) equals 6.2 and E minus 5. A and B are nearly congruent, with C and D showing slightly higher peaks at 700 and 725 kilonewtons at between 17 and 20 millimeters of displacement in the top example. Both graphs show similar configurations to the baseline models in figures 6 and 7. Decreasing the cohesion past 6.2E minus 7 did not show significant effects.

Figure 10. Angle of internal friction variations. Graphs.

Both these shear test simulations of angle of internal friction (horizontal axis ranging from 0 to 30 millimeters displacement) versus X-force (vertical axis ranging from 200 to 1000) or internal energy (vertical axis ranging from 0 to 14) show little variation. Shear forces increase for decreasing angles of internal friction. Line A (red) equals 63 degrees, Line B (green) equals 43 degrees, Line C (dark blue) equals 23 degrees, and Line D (turquoise) equals 3 degrees. In the X-force graph A, B, and C are nearly congruent until 17 millimeters of displacement, with D showing a slightly higher peak at 800 kilonewtons. Both graphs show similar configurations to the baseline models in figures 6 and 7.

Return to Figure 10

Figure 11. Hyperbolic approximation of Mohr-Coulomb. Diagram.

The Mohr-Coulomb failure criterion is shown as a pyramid on its side that starts at the horizontal axis and ascends and descends to the right, crossing the vertical axis. The point where the line cuts the yield surface axis corresponds to the hexagonal Mohr-Coulomb pyramid, where the gradient of the yield surface is undefined. To avoid such angularity, Drucker-Prager is an inscribed cone that has smoothed ridge corners. In the diagram, this hyperbola fits inside the Mohr-Coulomb pyramid to give better approximations of real failure conditions. Significant distances are labeled A (the distance between the peak of the pyramid and peak of the cone), B (the distance from the X-axis to where the cone intersects the Y-axis), and D (the distance from the peak of the pyramid on the X-axis to where the line intersects the Y-axis).

Return to Figure 11

Figure 12. AYHP influence on yield surface. Diagram.

The diagram shows the effect of choosing different values for the hyperbolic coefficient, AHYP. Values greater than zero place the curve inside the Mohr-Coulomb pyramid, and values less than zero place it outside the pyramid. Significant distances are labeled AHYP subscript 1 (the distance between the peak of the pyramid and peak of the inside cone), AHYP subscript 2 (the distance between the peak of the pyramid and peak of the outside cone), B subscript 1 (the distance from the X-axis to where the inside cone intersects the Y-axis), B subscript 2 (the distance from the X-axis to where the outside cone intersects the Y-axis), and D (the distance from the peak of the pyramid on the X-axis to where the line intersects the Y-axis). According to the developer, this diagram is misleading, because the modified yield surface (hyperbolic approximation) cannot be set outside the Mohr-Coulomb pyramid. Since parameter AHYP is squared in the yield function, it is always the absolute value. The developer modifies this diagram with this caution, “The input parameter AHYP should not be set so large that the initial stress lies outside of the modified Drucker-Prager surface (hyperbolic approximation surface).”

Return to Figure 12

Figure 13. AYHP parameter variations. Graphs.

Both these shear test simulations of AHYP parameter variations (horizontal axis ranging from 0 to 30 millimeters displacement) versus X-force (vertical axis ranging from 200 to 1000) or internal energy (vertical axis ranging from 0 to 20) show little variation. Line A (red) equals negative 1.0E minus 7, Line B (green) equals 0.0, Line C (dark blue) equals 1.0E minus 7, and Line D (turquoise) equals 1.0E minus 3. In the both graphs A, B, and C are congruent. The X-force graph is a curve that peaks at 18 millimeters at 700 kilonewtons for all AHYP values except Line D, which rises steeper to 825 kilonewtons at 16 millimeters displacement. For internal energy, the curves for A through C rise to 11 internal energy at 25 millimeters, while Line D is steeper and ends at 20 internal energy at 28 millimeters displacement. Significant increases in forces were found when AHYPwas increased. This parameter does not have significant variations from Mohr-Coulomb when set to small values (1.0E minus 7).

Return to Figure 13

Figure 14. Itermax parameter variations. Graphs.

Both of these shear test simulations of itermax (plasticity iterations) parameter variations (horizontal axis ranging from 0 to 25 millimeters displacement) versus X-force (vertical axes ranging from 200 to 800 and from 0 to 600) show little variation. In the first graph where cohesion equals 6.2E minus 6, Line A (red) equals 1, Line B (green) equals 5, Line C (dark blue) equals 10, and Line D (turquoise) equals 100. All lines are practically congruent, peaking at 700 kilonewtons at 18 millimeters displacement and ending at 400 kilonewtons at 27 millimeters displacement. Line D drops off at 675 kilonewtons at 17 millimeters displacement. In the second graph, cohesion equals 0.0, Line A equals 10, Line B equals 50, and Line C equals 100. All the lines plot practically congruent with peak values of 700 kilonewtons at 14 millimeters displacement, and falling to 400 kilonewtons at 24 millimeters displacement.

Return to Figure 14

Figure 15. DS-4 simulation results from the developer. Graph.

The developer used LS-DYNA to prepare this graph comparing the physical test results with the modeled shear test DS-4. The horizontal axis is deflection in millimeters ranging from negative 1.00E plus 01 to 8.00E plus 01. The vertical axis is shear force, expressed in kilopascals, ranging from 0 to 60. However, a note cautions that shear force is really shear stress. It is assumed that the shear stress was made by dividing the shear force (the X-direction cross-sectional force defined in the model) by the original cross-sectional area of the sample (0.2 square meters). The test data is a red line that peaks between 1.00E plus 01 and 2.00E plus 01 at 45 kilopascals and then levels off to 35 kilopascals. The blue LS-DYNA line rises higher more gradually to a peak of 52 kilopascals between 4.00E plus 01 and 5.00E plus 01, but ends abruptly because the model was stopped at 46 millimeters deflection. To the user, this situation suggests the developer’s graph represents instability in the soil material model, but the latter disputes this assertion.

Return to Figure 15

Figure 16. Instability of developer’s soil model. Drawing.

The three-dimensional drawing shows horizontal and vertical shooting nodes in the model that represent failure modes for some of the material parameters. The user’s theory is that damage routines caused the instability, and thus parameters were altered to find the instability source. Suspect parameters include PHIRES (residual shear strength), DINT (volumetric strain), VDFM (void formation energy), DAMLEV (deletion damage), and EPSMAX (failure strain).

Return to Figure 16

Figure 17. Direct shear results with the developer’s soil parameters. Graph.

This figure shows the shear stress displacement for the run in figure 16. The horizontal axis is displacement in millimeters ranging from 0 to 80. The vertical axis is shear stress in kilopascals ranging from 0 to 60. The shear stress rises from zero to 41 at 30 millimeters displacement. The shear stress calculated on the user’s computer was significantly less than that calculated on the developer’s computer.

Return to Figure 17

Figure 18. Hourglass control type 1, qm = 0.1. Graph.

The graph compares Line A (red) internal energy to Line B (blue) hourglass control coefficient. The horizontal axis is time (in milliseconds) ranging from 0 to 40. The vertical axis is matsum data ranging from 0 to 300. With the hourglass default value of 0.1, both lines plot in a similar pattern. The internal energy rises to 250 at 47 seconds, while the hourglass energy is a flatter curve, only reaching 150 at 47 seconds.

Return to Figure 18

Figure 19. Hourglass control type 4, qm = 0.005. Drawing.

This three-dimensional graphic representation of attempts at hourglass control shows significant distortion to the cylinder, with many of the mesh components skewed and out of alignment. A mesh refinement study using finer mesh compromised the contact between the soil and the direct shear testing device, so no firm results were obtained.

Return to Figure 19

Figure 20. Lagrangian approach: Stable large deformations, modified honeycomb material. Drawings.

Soils in roadside safety applications require large deformation capabilities. The user suggests that, although this requirement is difficult to meet for Lagrangian-based codes, type-126 material (modified honeycomb) allows large deformations that maintain stability without presenting analysis difficulties. The first model shows the steady state with a pink half-cylinder on top of the blue oblong, then simulates extreme crushing in two stages with deformation of the base material and the mesh compromised to where the half-cylinder is imbedded in the oblong. The developer counters that this logic is irrelevant, because the material crushing example involves only deformation of a continuous material without the complication of material separation that must occur in a direct shear test. The developer contends that this example would have to shear a hole in the crushing material rather than just crushing it without developing new material surfaces.

Return to Figure 20

Figure 21. Taylor problem. Drawings.

The Taylor problem is a classic example of large distortions. The first model shows a square with red grids as the initial condition. The Lagrangian approach shows highly distorted elements in the mesh that lead to calculation inaccuracies. In the third example, the ALE formulation maintains uniform mesh that prevents local instabilities caused by distorted elements. This theory was not tested with the FHWA soil material model.

Return to Figure 21

Figure 22. Multimaterial Eulerian formulation with Lagrangian coupling. Drawings.

The multimaterial Eulerian formulation allows materials to exist within each solid element and lets the material flow from element to element. The solid element mesh is fixed in space in this approach, because Langrangian coupling allows structural elements to be placed in the Eulerian mesh. Contact algorithms handle the interaction between the fluid Eulerian and Lagrangian structural elements. The contrived example is a two-dimensional pink post in red soil with a blue half-circle acting as the shear source. In the three-dimensional example, neither the grid in the light blue posts nor the pink half-cylinder mesh is distorted. This theory was not tested with the Federal Highway Administration soil material model.

Return to Figure 22

Figure 23. Direct shear test. Drawings.

The top half of these models duplicate figure 5, which is the modeling for the shear test. The two-part computer model shows the test cylinders unchanged and plumb (green stacked exactly on gray at the outset of the test, which represents 0 milliseconds). In the second example, at 49.992 milliseconds, the top green cylinder is deflected to the left about 50 millimeters. In the left drawing, each surface cylinder of soil begins in contact with the other, and at the end of the time analysis, they are contacting metal and air. The area between the two cylinder halves becomes a noncontinuous slide surface. In the second model, which represents physical testing, at the beginning of the test the red cylinder is normal, but becomes skewed to the left with its grids distorted. The user’s conclusion is that the soil deforms uncharacteristically because the model does not replicate the shear seen in physical testing.

Return to Figure 23

Figure 24. Soil modulus failure test. Drawings.

The model shows a green post in red soil with the test mechanism applying pressure to test soil modulus failure. In this comparison between the model and physical testing, the post bends in the simulation while it rotates in the physical testing. This anomaly makes the user question the validity of the model.

Return to Figure 24

Figure 25. Soil shear failure test. Drawings.

The model shows a green post in red soil with the test mechanism applying pressure to test shear failure. In this comparison between the model and physical testing, the post bends in the simulation while it rotates in the physical testing. This anomaly makes the user question the validity of the model.

Return to Figure 25

Figure 26. Hydrostatic tension: Internal energy. Graph.

In this graph, the horizontal axis is time ranging from 0 to 50 milliseconds and the vertical axis is internal energy (E minus 08). The results are nearly identical until 35 milliseconds, when the SGI computer shows a steeper curve reaching 0.9 at 50 milliseconds and PC results diverge to 0.7 at 50 milliseconds. After the elements begin to fail and plasticity becomes significant, the conjecture is that roundoff errors between the computer platforms begin to influence the hydrostatic tension results.

Return to Figure 26

Figure 27. Triaxial compression results: Internal energy. Graph.

In this graph, the horizontal axis is time in milliseconds ranging from 0 to 100. The vertical axis is internal energy (E minus 03) ranging from 0 to 1.4. For the internal energy triaxial results, the runs are nearly identical, rising in a curve to 1.3 at 85 milliseconds, then drop vertically to zero, except that the single element in the user’s SGI run fails a few cycles before the PC runs.

Figure 28. Triaxial compression results: Effective stress. Graph.

In this graph, the horizontal axis is time in milliseconds ranging from 0 to 100. The vertical axis is effective stress (V-M) ranging from 0 to 0.012. The runs are nearly identical, rising in a curve to 0.011 at 40 milliseconds, remaining level until 85 milliseconds, then dropping vertically to zero. The single element in the user’s SGI run fails a few cycles before the PC runs.

Return to Figure 28

Figure 29. Internal energy. Graph.

This graph shows how the results from all computers used for testing internal energy match well up to 260 milliseconds of simulation. The horizontal axis is time in milliseconds from 0 to 300, and the vertical axis is internal energy ranging from 0 to 10. When the elements begin to fail, the results start to diverge. The user’s SGI results plot vertically from 6 to 8 starting at 20 milliseconds, while the PCs show more congruent plots. Overall agreement between results is acceptable.

Return to Figure 29

Figure 30. Cross section force through cylinder. Graph.

This graph shows how the results from all computers used for testing internal energy match well up to 260 milliseconds of simulation. The horizontal axis is time in milliseconds from 0 to 300, and the vertical axis is resultant force ranging from 0 to 400. When the elements begin to fail, the results start to diverge. The user’s SGI results begin to drop sooner at 250 milliseconds, while the PCs show more congruent plots. Overall agreement between results is acceptable.

Return to Figure 30

Figure 31. Deformed geometry. Drawing.

These models show how the results from all computers used for testing match well up to 260 milliseconds of simulation. When the elements begin to fail, the results start to diverge. The developer’s PC shows deformation in the upper right corner, the user’s SGI results show a deformed left side and buckling in the middle, and the user’s PC also shows deformation in the right corner. However, overall agreement between results is acceptable.

Return to Figure 31

Equations

Equation 1. Mass equals PV, which equals the sum of 2.35E minus the quotient of 6 kilograms divided millimeters, cubed, times 1.0 millimeters, cubed, which equals 2.35E minus 6 kilograms.

Equation 2. Sigma bar subscript VP equals the sum of 1 minus zeta times sigma bar plus sigma bar subscript trial.

Equation 3. Zeta equals nu divided by the sum of the change in T plus nu.

Equation 4. Nu equals the quotient of gamma times R divided by epsilon dot, all raised to the power of the quotient of V subscript N minus 1 divided by V subscript N.

Equation 5. E equals 3K times the sum of 1 minus 2V.

Equation 6.E equals 1.5K.

Equation 7. K equals the quotient of 2E divided by 3.

Equation 8. G equals E divided by the product of 2 times the sum of 1 plus V.

Equation 9. 2G times the sum of 1 plus V equals 3K times the sum of 1 minus V.

Equation 10. G equals K times the quotient of 3 times the sum of 1 minus 2V divided by the product of 2 times the sum of 1 plus V.

Equation 11. G equals K times the quotient of 3 times the sum of 1 minus 2 times 0.25 divided by the product of 2 times the sum of 1 plus 0.25.

Equation 12. G equals K times the quotient of 3 times the sum of 1 minus 0.50 divided by the product of 2 times 1.25. This then means that G equals K times the quotient of 1.50 divided by 2.50, which means that G equals 0.6K.

Equation 13. Tau subscript phi (the angle of friction) equals cohesion, C, plus normal stress, sigma, times the tangent of phi.

Equation 14. The yield surface, sigma subscript Y, equals negative pressure, P, times the sine of the angle of friction, phi, plus the square root of the sum of the second invariant of the stress deviator, J subscript 2, times the function of the angle in the deviatoric plane, K of theta, squared, plus the product of the Drucker-Prager coefficient, AHYP, squared, times the sine, squared, of phi. Then subtract cohesion, C, cosine, which equals 0.

Equation 15. AHYP equals C divided by 20 times the cotangent of phi.

Equation 16. K of theta equals the quotient of 4 times the sum of 1 minus E squared, times the cosine, squared, of theta, plus the square of the sum of 2E minus 1, all divided by the following: 2 times the sum of 1 minus E squared, times the cosine of theta, plus the sum of 2E plus 1, times the square root of the following: 4 times the sum of 1 minus E squared, times the cosine, squared, of theta, plus 5E, squared, minus 4E.

Equation 17. The cosine of 3 times theta equals 3 times the square root of 3 times J subscript 3, all divided by 2 times the square root of J subscript 2, cubed.

Equation 18. Change in phi equals E subscript T times the sum of 1 minus the quotient of phi minus phi subscript init divided by A subscript N times phi subscript max, all times the change in epsilon subscript eff plastic.

Equation 19. K equals nonporous bulk modulus, K subscript I, divided by the sum of 1 plus the product of K subscript I times the parameter controlling the stiffness before the air voids are collapsed, D subscript 1, times the current porosity, N subscript cur.

Equation 20. U equals nonporous bulk modulus, K subscript SK, divided by the sum of 1 plus the product of K subscript SK times the parameter for pore-water pressure before the air voids are collapsed, D subscript 2, times the current porosity, N subscript cur, all times the total volumetric strain, epsilon subscript V.

Equation 21. B equals 1 divided by the sum of 1 plus the product of N times the quotient of K subscript SK divided by K subscript V.

Equation 22. PWD2 equals D subscript 2, which equals K subscript SK divided by the sum of 1 plus the product of K subscript SK times the parameter for pore-water pressure before the air voids are collapsed, D subscript 2, times the current porosity, N subscript cur, all times the total volumetric strain, epsilon subscript V.

Equation 23. Residual shear strength, S subscript residual, equals effective stress, sigma prime, times the tangent of the residual angle of friction, phi subscript ult.

Equation 24. Alpha equals the quotient of 2 times G subscript F divided by K times initial damage threshold strain, Xi subscript O, times V raised to the one-third power, all added to Xi subscript O.

Equation 25. D subscript max equals the quotient of the sum of the sine of phi minus the sine of phi subscript res, all divided by the sine of phi.

Table of Contents