Skip to contentUnited States Department of Transportation - Federal Highway Administration FHWA Home
Research Home
Report
This report is an archived publication and may contain dated technical, contact, and link information
Publication Number: FHWA-HRT-04-095
Date: November 2004

Manual for LS-DYNA Soil Material Model 147

PDF Version (1.96 MB)

PDF files can be viewed with the Acrobat® Reader®

CHAPTER 1. THEORY MANUAL

INTRODUCTION

This document is the final report for the development of the Federal Highway Administration's (FHWA's) soil model implemented into LS-DYNA. This report is in three sections: (1) a description of the justification and detailed theory of the model, (2) a user's manual that was submitted for inclusion in the LS-DYNA user's manual, and (3) examples that show the expected results of the model.

The original research plan was submitted to Martin Hargrave of FHWA and the FHWA organized Centers of Excellence in Finite Element Crash Analysis in September 1999. The material model was developed and changes requested by the centers of excellence were implemented in 2000. During this time, the user-defined material models in LS-DYNA were used. The model was verified and preliminary validation took place in 2001. The model was implemented as a standard material model in version 970 of LS-DYNA in February 2002. The soil model in the production version was checked against analyses done with the user-defined version in spring 2002. A user ran identical analyses and other analyses to validate the model in summer 2002.

A significant difference between the soil material model discussed in this report and the wood material model discussed in another report is that the wood material model was based on extensive experimental data. In the case of the soil model, there was no material property data available for the applications needed (road-base materials). This lack of data caused the material model to be developed based on the one set of data available and the general behavior of cohesionless soils. In addition, some behaviors could not be validated and default properties for the soils used in road-base testing could not be uniquely defined.

The goal of this research was to develop a soil material model in LS-DYNA that will represent the soils used for National Cooperative Highway Research Program (NCHRP) 350 roadside safety hardware testing. In this section of the report, the theory for the soil material model is presented. Investigative work done prior to development of the model is discussed first. The investigative work was for the purpose of determining the critical behaviors of the material that must be included in a model and establishing whether an existing model could be enhanced to meet these requirements or whether a new material model had to be developed. The critical behaviors and how they can be modeled, followed by details of a review of currently used soil material models, are also discussed. Also in this section, model development areas, including details of the algorithms and implementation, are described.

PRELIMINARY WORK

The investigative work included determining the critical behaviors of NCHRP 350 soil. Material models already in LS-DYNA were investigated to see if one of them, if enhanced, would be suitable, and what enhancements would be necessary for the modeling of soils in roadside safety applications.

Determination of Critical Behaviors

Behaviors that were critical to soil modeling in roadside safety applications were determined through discussions with roadside safety testers and analysts, by performing literature reviews, and by studying road-base laboratory test results. Most soil data that have been used for past soil model development in LS-DYNA are from laboratory tests that have relatively high confinement.(5) For roadside safety applications, the soil will have low or no confinement. Therefore, the laboratory data used to evaluate and develop the model must be at low or no confinement. Maximum pressures will be less than 30 megapascals (MPa).(6)

Most of the basic triaxial shear test data (from the laboratory) used to determine critical behaviors and develop algorithms have been from the U.S. Army Corps of Engineers (USACE) geological database. The specific soil used is a crushed limestone road base that conforms to NCHRP 350 soil, grade B (see table 1). However, the larger aggregate was removed so that uniform stresses/strains could be achieved in standard specimen sizes for triaxial shear tests. Additional data from the roadside safety testing subcontractor were used to verify and validate the new soil model developed for LS-DYNA.

Elastic behavior of soil is isotropic. This requirement is based on the fact that road-base soils are well graded and not stratified. The standard soil is considered to be cohesionless (i.e., it has no tensile strength). This behavior is common to soils that contain little clay. Elastic behavior mainly affects unloading and the isotropic compression (volumetric) behavior because road-base material tends to have very low shear strength at low confinement (< 0.5 MPa). Very little shear stress is needed to initiate nonrecoverable energy dissipation (plasticity and damage).

Soil behavior is greatly affected by void ratio, compaction, and excess pore-water pressure. Experimental data for a standard road base show that undrained conditions (high moisture content) can greatly affect the amount of deformation. Compaction tends to increase the initial yield strength and the ultimate strength. The void ratio is directly related to compaction. Reduction in the void ratio increases the strength and the bulk modulus of the soil.

Figure 1 shows two commonly used models for yield surfaces. The three axes shown are the three principal stress axes. The generator axis for each surface is the pressure axis. The surface on the left is the Mohr-Coulomb model (a typical soil yield surface) and, for comparison, the surface on the right is the standard Von Mises model (a typical metal yield surface). The yield (onset of plasticity) and ultimate (peak) strength of soil are pressure-dependent. This means that the plasticity surface is dependent on both the pressure and the shear stresses. This behavior differs from metals where the plasticity surface is only a function of the shear stresses. Experimental evidence for cohesionless soil shows that, at low pressures, the yield surface is triangular in the deviatoric plane as shown in figure 2.(3) For the Waterways Experiment Station (WES) road-base soil, the ratio of minimum stress (triaxial extension) to maximum stress (triaxial compression) in the deviatoric plane is 0.70.

Table 1. Gradation data for USACE road-base soil tests.

Sieve Size (mm) Percent Passing
9.5200 100.0
6.3500 90.0
4.7600 79.0
3.3600 66.0
2.3800 54.0
2.0000 49.0
1.1900 36.0
0.8400 32.0
0.5900 28.0
0.4200 24.0
0.2970 21.0
0.2100 18.0
0.1490 15.0
0.1050 13.0
0.0740 11.0

 

Figure 1. (A) Pressure-dependent (Mohr-Coulomb) and (B) pressure-independent (Von Mises) yield surfaces. Diagram. This diagram shows two figures, the first, (A), representing pressure-dependent (a typical soil yield surface) yield surface, is cone shaped, with its narrowest point beginning at the conjunction of the three principal stress axes, and projecting outward from there, the cone widening. The second figure, (B), representing pressure-independent yield surface (a typical metal yield surface) is cylindrical shaped, one end simply centered at the point of conjunction of the three principal stress axes. Figure b pressure-independent (Von Mises) yield surfaces

Figure 1 (a) Pressure-dependent (Mohr-Coulomb) and (b) pressure-independent (Von Mises) yield surfaces.

Figure 2. Yield surface in deviatoric plane for cohesionless soils. Diagram. This diagram shows one triangular figure set inside of a larger of the same. The point of conjunction of the axes is centered in the smaller of the triangles. The three starting points of the axes are numbered 1 through 3 and protrude through the triangle at the its corners with number 1 at the top. All three corners of each triangle are rounded off. There are 10-minute circles displaced along each of the triangles three walls, representing data from Monterey Number 0 Sand.

Figure 2. Yield surface in deviatoric plane for cohesionless soils.

Figure 3 shows the peak principal stress difference versus the average normal stress for the crushed limestone road-base material. The peak principal stress difference is also the maximum shear strength at the given pressure. This plot presents the shear strength as a function of pressure for the road-base material. Notice that for the pressure range indicated in the figure, the shear strength of this soil varies linearly with pressure.

Figure 3. Principal stress difference (peak shear strength) versus pressure (average normal stress) for road-base material. Graph. This graph shows a rising red diagonal line. The vertical axis of this graph ranges from 0 to 35 and represents the stress differential in megapascals, while the horizontal axis of this graph ranges from 0 to 20 and represents pressure in megapascals. The diagonal line begins at the 0 points of both the vertical and horizontal axes. The red line stretches diagonally for most the length of the graph, stopping at the points of the 30 mark on the vertical line and the 18 mark on the horizontal line. There are small red circles placed along the length of the red line. The first exists at the points of 2 on the vertical axis and 1 on the horizontal axis. The second is placed at 10 on the vertical axis and 6 on the horizontal axis. The third finds itself at the 28 mark on the vertical axis and just before the 17 mark on the horizontal axis. The final circle sits at the end of the diagonal red line on the points of 30 on the vertical axis and 18 on the horizontal axis.

Figure 3. Principal stress difference (peak shear strength) versus pressure (average normal stress) for road-base material.

At high confinement, the road-base soil can have significant peak shear strength (> 120 MPa).

At low confining pressures, the standard soil dilates (expands) near peak shear stresses. At high confining pressures (> 100 MPa), the standard soil will stop expanding. This change in volumetric behavior is one of the reasons for employing a cap model. For roadside safety analysis applications, there are no pressures greater than 30 MPa; thus, it is believed that a cap model is not needed. A fully associative plasticity model predicts dilation of the material after the yield strength is reached.

The standard soil at low confining pressures typically exhibits strain softening. That is, at pressures below 30 MPa, the soil strain hardens from the yield stress to the ultimate stress, then the strain softens. Strain softening can be a major source of energy dissipation; thus, it is believed that it must be included in the material model.

The strength of the soil increases at high strain rates. No experimental strain rate data were found for road-base materials, but data for cohesionless (sandy) soils show significant strain-rate effects. However, we are not convinced that strain-rate behavior is a critical behavior, because applicable rates are relatively low for roadside safety applications.

The moisture content of the soil can affect the elastic moduli, the shear strength, and the softening behavior of the soil.

Figure 4 shows the force-deflection curves measured in two steel posts in NCHRP 350 soil bogie tests.(8) The tests were identical except for the moisture content (5 percent versus 26 percent). The peak force for the relatively dry soil was much higher, and the stiffness was much greater. Figure 5 shows the energy absorbed during the two tests. The amount of energy absorbed at a given deflection was much larger for the low-moisture-content NCHRP 350 soil.

The effects of moisture are complicated and are different for different soil types. For instance, granular soils with low relative densities show little effect on bulk modulus, while clayey soils show significant effects. Fortunately, the NCHRP 350 soil is granular and the test facilities typically run their tests at a low moisture content (3 percent to 7 percent). However, we believe that the ability to simulate actual field conditions is important, so techniques to simulate moisture effects were implemented. The degree of saturation and the void ratio are critical parameters in the determination of moisture effects. Moisture effects on shear strength can be introduced by including an excess pore-water pressure algorithm.

Compaction is the decrease in the void ratio and the increase in the relative density. Compaction in granular soils tends to increase the shear strength slightly and lessen the amount of volumetric strain that causes the onset of excess pore-water pressure effects.

Figure 4. Force deflection for two steel posts in soil tests with different moisture contents (5 percent and 26 percent). Graph. This graph shows two lines, one red that depicts 5 percent moisture content and one black that depicts 26 percent moisture content. The vertical axis of this graph ranges from a negative 10 to a positive 80 and represents force in kilonewtons, while the horizontal axis ranges from 0 to 700 and represents deflection in millimeters. Both red and black lines have the same starting point, 0 on the vertical axis, though the red line spikes immediately to 70 while the black line takes slightly longer before spiking only to 30, it highest peak. From its highest peak at 70, the red line then steadily decreases, and the line finally stops near the 250 mark on the horizontal axis. The black line continues to decrease and increase, though it remains mostly steady along the length of the graph, where it finally drops off at the 700 mark on the horizontal axis.

Figure 4. Force deflection for two steel posts in soil test with different moisture contents (5 percent and 26 percent).

Figure 5. Energy versus deflection for two steel posts in soil bogie tests with different moisture contents. Graph. This graph shows two roughly diagonal lines, one red that depicts 5 percent moisture content and one black that depicts 26 percent moisture content. The vertical axis of this graph ranges from 0 to 20 and represents energy absorbed in kilojoules, while the horizontal axis ranges from 0 to 700 and represents deflection in millimeters. The red line begins near the 20 mark on the horizontal axis and spikes to 17 on the vertical axis, where it also stops at the 400 mark on the horizontal axis. The black line begins near the 80 mark on the horizontal axis and peaks at the 15 mark on the vertical axis, where it begins to straighten out and finally ends at the 700 mark on the horizontal axis.

Figure 5. Energy versus deflection for two steel posts in soil bogie test with different moisture contents.

Evaluation of the Utility of Models Already in LS-DYNA

Several material models available in LS-DYNA were reviewed as possible soil model candidates, from the simplest (material model 5) to the most complex (material model 25). Most of the candidate models in LS-DYNA are extensions of material model 5 (soil and foam). Two of the exceptions to this are model 16 (pseudo-tensor model) and model 25 (geological cap model). Model 5 and its extensions, model 14 (soil and foam with failure), and model 79 (hysteretic soil) are basically an analytical pressure-dependent yield surface. They all must have confinement to be stable (see the LS-DYNA user's manual for model 5). For the roadside application, the top surface of the soil is not confined and is at zero pressure during a significant portion of the analysis. However, both model 5 and model 79 were evaluated with single-element runs. Both of these models were indeed found to be unstable in unconfined states. Material model 14 and model 78 (soil and concrete) are for analyses in high-pressure regions. Model 14 was evaluated previously and it was found that it does not simulate low-/zero-pressure behavior well.

Material model 25 (geological cap model) was also evaluated. This model is complex, but does handle low-confinement behavior well. This model was first evaluated with single-element runs. It performed well, it was stable (see figure 6), and it predicted peak stress level accurately. However, model 25 cannot simulate strain softening. Despite this shortcoming, it is the only existing model that could simulate much of the basic behavior needed for roadside safety applications. However, model 25 is very inefficient for this application since it uses a cap surface. Because of the relatively low confinement that soil has in roadside safety applications, a cap surface is not needed

Because of the cap/shear surface intersection (non-smooth) in model 25, there are many corners in the yield surface. These corners cause the algorithm to be very complex and inefficient. In addition, model 25 would need extensive enhancements, including:

  • Three invariant yield surfaces instead of two invariant surfaces to simulate lower strength in triaxial extension than in triaxial compression and triaxial shear.

  • Smooth behavior of the yield surface at very low shear stresses.

  • Strain softening behavior.

  • Isotropic hardening behavior on the yield surface, instead of kinematic hardening.

  • Strain-rate-dependent strength enhancements.

  • Pore-water pressure behavior.

Based on these observations, it was concluded that none of the existing models are adequate, and it was decided to develop a new soil material model. The theoretical basis for this new model is presented in the next section.

Figure 6. Model 25 single-element run, triaxial compression at 3.4 megapascals compared to WES data. Graph. This graph shows two lines, one black that depicts WES 3.4, and the other red, which depicts DYNA. The vertical axis of this graph ranges from 0 to 15 and represents stress differential in megapascals, while the horizontal axis ranges from 0 to 20 and represents the percent of axial strain. The black line, beginning at 0, spikes immediately to the 12 mark on the vertical axis where it remains steady, stretched out until it stops at the 8 mark on the horizontal axis. The red line, also beginning at 0, curves up to the 12 mark on the vertical axis gradually where it remains constant across the length of the graph, finally stopping at the 11 mark and dropping straight down into the horizontal axis.

Figure 6. Model 25 single-element run, triaxial compression at 3.4 MPa compared to WES data.

MODEL DEVELOPMENT

The objectives of the material model development effort in order of priority were: accuracy, robustness, efficiency (speed), and ease of use. The choices made in the model are balanced between these objectives. The following subsections describe the characteristics of the algorithms necessary to model soil for roadside safety applications.

Elastic Constitutive Behavior

We assumed that the elastic properties of the soil are isotropic. Bulk and shear moduli were used as input parameters. Standard soil tests (i.e., triaxial shear tests and uniaxial strain tests) produce these parameters directly. To simulate the effects of voids, the bulk modulus was made to be a function of volumetric strain. As the volumetric strain increases, the modulus increases to simulate the collapse of voids and the stiffening of the material.

The effects of moisture content/excess pore pressure were also simulated with changes to the elastic moduli. As the remaining voids of the soil become filled with moisture, the material becomes more incompressible. To simulate the effects of excess pore-water pressure, a function that involves the nonporous bulk modulus (inverse of soil compressibility), the porosity, and the degree of saturation was used:

 
Equation 1. The bulk modulus, K, equals nonporous bulk modulus, K subscript I, divided by 1 plus the product of K subscript I times the material constant controlling the stiffness before the air voids are collapsed, D subscript 1, times the current porosity of the soil, N subscript cur.
(1)

where:

Ki = nonporous bulk modulus

ncu = current porosity = Max open bracket 0, (w minus epsilon) close bracket

w = volumetric strain corresponding to the volume of air voids =n(1- – S)

epsilon symbol v = total volumetric strain

D1 = material constant controlling the stiffness before the air voids are collapsed

n = porosity of the soil = equation 1b

e = void ratio = gamma subscript sp Rho subscript w (1 plus m subscriptc) over Rho - 1

S = degree of saturation = rho m subscript c over n rho subscript w (1 plus m subscript c)

rho symbol, gamma symbolsp, m c, P w = soil density, specific gravity, moisture content, and water density.

Figure 7 shows the effect of the D1 parameter on the pressure-volumetric strain relationship (bulk modulus). The elastic moduli are used to determine the elastic stresses and the elastic trial stresses. The bulk modulus is always a monotonically increasing value (i.e., j is the time-step index),

 
Equation 2. K subscript time index J plus 1 equals K subscript I divided by 1 plus the product of K subscript I times D subscript 1 times N subscript cur if the total volumetric strain, E subscript V subscript J plus 1 is greater than E subscript VJ, and K subscript time index J plus 1 equals K subscript J if E subscript V subscript J plus 1 is less than or equal to E subscript VJ.
(2)

Note that the standard practice of treating compressive stresses and strains as positive quantities is followed.

Figure 7. Pressure versus volumetric strain showing effects of the D subscript 1 parameter. Graph. This graph shows five separate lines. The first, black, labeled K subscript J equals 0.468 and D subscript 1 equals 0.68; the second, red, labeled K subscript J equals 0.468 and D subscript 1 equals 5.0; the third, orange, labeled K subscript J equals 2.5 and D subscript 1 equals 5.0; the fourth, green, labeled K subscript J equals 11 and D subscript 1 equals 5.0; and the fifth, blue, labeled K subscript J equals 15 and D subscript 1 equals 10.0. The vertical axis of this graph ranges from 0.0 to 2.0 and represents pressure, while the horizontal axis of this graph ranges from 0.00 to 0.20 and represents the volumetric strain. All five colored lines begin near the 0.0 mark and gradually rise. At the 0.11 mark on the horizontal axis, both the green and blue lines bend and drastically spike, the blue rising higher where it leaves the graph at the 1.5 mark on the vertical axis, while the green leaves the graph at the 1.2 mark on the vertical axis. The orange line bends only slightly at the 0.11 mark and rises gradually, where it leaves the graph just above the 0.3 mark on the vertical axis. Both black and red lines remain constant with each other, rising to their highest point at 0.1 on the vertical axis.

Figure 7: Yield Surface Behavior

The initial yield surface is where the soil initially starts to dissipate nonrecoverable strain energy. A well-documented yield surface used for soils is the Mohr-Coulomb surface. However, the standard Mohr-Coulomb surface has two significant deficiencies for our use:

The first deficiency is that the surface comes to a point (singularity) at the intersection with the pressure axis (zero shear strength). This region of the yield surface is critical in roadside safety applications because of the low confinement. This type of singularity can cause both numerical and efficiency problems in the plasticity algorithm. To ensure an accurate, robust, and efficient algorithm, the yield surface needs to be convex and smooth. Second, the Mohr-Coulomb surface is hexagonal or circular in the deviatoric plane (see figure 8). Based on the experimental evidence, the yield surface should be able to become triangular in shape at low confinement pressures. (9)

To correct these deficiencies, a modified Mohr-Coulomb surface was adopted. The yield surface was modified based on the work of Abbo and Sloan.(9) The standard Mohr-Coulomb yield surface, , is represented as:

 
Equation 3. The standard Mohr-Coulomb yield surface, F, equals negative pressure, P, times the sine of the internal friction angle, phi, plus the product of the function of the angle theta in the deviatoric plane, K theta, times the square root of the second invariant of the stress deviator, J subscript 2, minus the amount of cohesion, C, times the cosine of phi, which equals .
(3)

where:

P = pressure

phi symbol = internal friction angle

K (theta symbol) = function of the angle theta symbol in the deviatoria plane

square root of J subscript 2= square root of the second invariant of the stress deviator

c = amount of cohesion

The modified yield surface is a hyperbola fitted to the Mohr-Coulomb surface. At the crossing of the pressure axis (zero shear strength), the modified surface is a smooth surface. At this point, it is perpendicular to the pressure axis. The equation for the modified Mohr-Coulomb surface is:

 
Equation 4. F equals negative P times the sine of phi plus the square root of the following: J subscript 2 times K theta squared, plus the product of a parameter for determining how close the modified surface is fitted to the standard Mohr-Coulomb yield surface, A, squared, times the sine, squared, of theta. Then subtract C times the cosine of theta, and this equals 0.
(4)

Here, is a parameter for determining how close the modified surface is fitted to the standard Mohr-Coulomb yield surface. If is zero, then the standard Mohr-Coulomb surface is recovered. The input parameter should be set close to zero, based on numerical considerations.

Figure 9 shows the modified Mohr-Coulomb surface in shear stress versus pressure space. It is almost identical to the original surface, except at low shear stresses.

Figure 8. Standard Mohr-Coulomb yield surface in principal stress space. Diagram. This diagram is cone shaped, with its narrowest point beginning at the conjunction of the three principal stress axes, and projecting outward from there, the cone widens.

Figure 8. Standard Mohr-Coulomb yield surface in principal stress space.

Figure 9. Comparison of Mohr-Coulomb yield surfaces in shear stress-Pressure space (standard-(A) green, modified-(B) red). Graph. This graph shows two lines, one green that depicts standard yield surface, and one red that depicts modified yield surface. The vertical axis of this graph is represented by shear stress, while the horizontal axis is represented by pressure. There are no marks or specific counting measures on this graph. Though both lines run nearly parallel to each other as they rise diagonally along the graph, almost touching and end at the same spot on the graph. Their starting points differ, with the green line beginning on the vertical axis and the red line beginning on the horizontal axis.

Figure 9. Comparison of Mohr yield surfaces in shear stress — Pressure space (standard — (A)/green, modified —(B)/red.

To remedy the second deficiency, the shape in the deviatoric plane, the standard Mohr-Coulomb K (theta symbol) function was changed to a function used by Klisinski:(10)

 
Equation 5. The standard Mohr-Coulomb K theta function 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 2 times E minus 1, all divided by 2 times the sum of 1 minus E squared, times the cosine of theta, plus the sum of 2 times E minus 1 times the following, raised to the half power: 4 times the sum of 1 minus E squared, times the cosine, squared, of theta, plus 5 times E squared minus 4 times E.
(5)

Here, cos equation 5a = third invariant of the stress deviator, and e= material input parameter describing the ratio of triaxial extension strength to triaxial compression strength. If e is set equal to 1, then a circular cone surface is formed (see figure 8). If is set to 0.55, then a triangular surface is found (see figure 10). K (theta symbol) > is defined for 0.5,< e less than or equal to 1.0.

Figure 10. Yield surface with E equals 0.55. Diagram. This diagram is triangularly cone shaped in nature. Its narrowest point begins at the conjunction of the three principal stress axes and projects outward from there; the triangular cone widens.

Figure 10. Yield surface with e = 0.55.

Excess Pore-Water Pressure Behavior

For some soils, excess pore-water pressure can make a significant difference on the shear strength of the soil, especially for near-saturated conditions. Because of the loading typical for roadside safety applications (i.e., many soils have low permeability), water in the voids will not have time to flow; therefore, the water will cause an excess pressure increase as the impact load collapses the air voids in the soil. To simulate this behavior, a standard (practical) soil mechanics technique(11) is used for reducing the total pressure, , by the excess pore-water pressure, , to get an "effective pressure":

 
Equation 6. Effective pressure, P prime, equals total pressure, P, minus excess pore-water pressure, U.
(6)

 

Figure 11. Effects on pressure because of pore-water pressure. Diagram. This graph shows one straight line with a positive slope equal to the internal angle of friction. This line sets the yield strength of the soil as a function of pressure. The vertical axis of this graph is the yield strength, denoted as the square root of J subscript 2. The horizontal axis is pressure. Two vertical lines extend from the horizontal axis to the sloped line. One line is located at the effective pressure, denoted P superscript prime. The other line is located at the total pressure P. The total pressure is larger than the effective pressure. The pressure difference between these lines (pressures) is designated as the excess pore water pressure U. The graph demonstrated the yield stress at the effective pressure is lower than the yield stress at the total pressure.

Figure 11. Effects on pressure because of pore-water pressure.

Figure 11 shows how pore-water pressure affects the algorithm for the plasticity surface. The excess pore-water pressure reduces the total pressure, which will lower the shear strength, square root of J subscript 2 . A large excess pore-water pressure can cause the shear strength to become zero.

The water in the voids of the soil causes the excess pore-water pressure. As the air void volume is reduced to zero during loading, the pore-water pressure increases. The water in the remaining voids causes the effective load on the soil particles to be reduced.

To calculate the pore-water pressure, , an expression similar to the equation used for the moisture effects on the bulk modulus was used:

 
equation 7
(7)

where:

K sk = bulk modulus for soil without air voids (skeletal bulk modulus)

n cur = current porosity = Max[0, (w - Ev) ]

w    = volumetric strain corresponding to the volume of air voids

      = n(1-S)

epsilon symbolv = total volumetric strain

D2 = material constant controlling the pore-water pressure before the air voids are collapsed

n = porosity of the soil = equation 1b

e = void ratio = gamma subscript sp Rho subscript w (1 plus m subscriptc) over Rho - 1

S = degree of saturation = rho m subscript c over n rho subscript w (1 plus m subscript c)

rho symbol, gamma symbolsp,m = soil density, specific gravity, and moisture content, respectively

Pore-water pressure is not allowed to become negative ( u greater than or equal to 0).

Figure 12 is a plot of the pore pressure versus volumetric strain for different parameter values. With the D2 parameter set relatively high compared to Ksk, there is no pore pressure until the volumetric strain is greater than the strains associated with the air voids. However, as D2 is lowered, the pore pressure starts to increase before the air voids are totally collapsed. The Ksk parameter affects the slope of the post-void collapse pressure-volumetric behavior

Parameter D2 can be found from Skempton pore-water pressure parameter B, where B is defined as:(11)

 
Equation 8. Skempton pore-water pressure parameter B equals 1 divided by the sum of 1 plus N times the quotient of K subscript SK divided by K.
(8)
 
Equation 9. Therefore parameter D subscript 2 equals the quotient of 1 minus B divided by B times K subscript SK times the product of N times 1 minus S.
(9)

where: Ksk=bulk modulus of the soil without air voids

This method does not include the dissipation of excess pore-water pressure as a function of time. The rate of dissipation can be a function of the loading rate and soil parameters, such as permeability. However, at this time, the lack of experimental data on road-base material on typical roadside tests would make this type of dissipative model useless. However, if experimental data and the required parameters for roadside tests become available, then a dissipative model could be easily inserted.

Figure 12. Effects of parameters on pore-water pressure. Graph. This graph shows three lines, the first is black and labeled K subscript SK equals 15 and D subscript 2 equals 1500; the second is red and labeled K subscript SK equals 15 and D subscript 2 equals 15; and the third is green and labeled K subscript SK equals 30 and D subscript 2 equals 15. The vertical axis of this graph ranges from 0.0 to 3.0 and represents pressure, while the horizontal axis of this graph ranges from 0.00 to 0.20 and represents volumetric strain. The green and red lines overlap from the beginning at 0.0/0.00, with a gradual increase of 0.1 on the vertical axis to the point of 0.11 on the horizontal axis. From that point, both lines spike, (the green more drastically than the red). The green line leaves the graph at 2.8 on the vertical axis, while the red leaves the graph at 1.4 on the vertical axis. The black line appears to begin at the 0.11 mark on the horizontal axis and diagonally increases until it leaves the graph at 1.3 on the vertical axis.

Figure 12. Effects of parameters on pore-water pressure.

Strain Hardening Behavior

Figure 13 shows the principal stress difference versus principal strain difference results of a triaxial compression test at a lateral stress of 6.9 MPa for a standard soil. The unloading portion of the curve shows that there was very little elastic (recoverable) strain in this test. The nonlinear part of the loading portion of the curve is pre-peak (plastic) hardening. The amount of hardening increases as the lateral confinement increases (i.e., there is less hardening at lower confining pressures). To simulate this nonlinear hardening behavior, the friction angle was increased as a function of the effective plastic strain:

 
Equation 10. The change in friction angle, phi, equals the input parameter, H, times the sum of 1 minus phi minus phi subscript init divided by N times phi subscript max, times the change in effective plastic strain, E subscript eff plas.
(10)

where:

epsiloneff plas = effective plastic strain

N = fraction of the peak strength internal friction angle where nonlinear behavior begins (0 < N less than or equal to 1 )

H (the input parameter) determines the stiffness of the nonlinear hardening

Figure 14 shows the effect on the yield surface of an increase in phi symbol.

Figure 13. Principal stress difference versus principal strain difference for triaxial compression test at sigma subscript 2 equals 6.9 megapascals of WES road-base material. Graph. This graph shows one black line. The vertical axis of this graph ranges from 0 to 30 and represents principal stress differential in megapascals, while the horizontal axis ranges from 0 to 40 and represents the percent of principal strain differential. The line rises along the vertical axis to the 5 mark, where it then curves upward and peaks at the points of 25 on the vertical axis and 15 on the horizontal axis. The line then begins its decent, where it stops at the points of 20 on the vertical axis and 35 on the horizontal axis. From there, the line drops straight down, leaning backward, where it leaves the graph at the point of 33 on the horizontal axis.

Figure 13. Principal stress difference versus principal strain difference for triaxial compression test at sigma symbol 2=6.9 MPa of WES road-base material.

Figure 14. Hardening of yield surface. Diagram. This diagram shows two images, both triangular and cone-shaped in fashion, appearing to be identical in shape, with one residing inside the other. Both narrowest points of these images begins at the conjunction of the three principal stress axes and project outward from there, where both images widen.

Figure 14. Hardening of yield surface.

Strain Softening Behavior

Figure 15 is a plot of the principal stress difference versus axial strain result for the same experiment as shown in figure 13. The principal stress difference softens (i.e., decreases) after it has reached its peak. The area under the curve in figure 15 after the peak stress is reached is the strain energy dissipated by the material because of strain softening. The strain energy dissipated in this post-peak region is almost as great as the strain energy dissipated in the pre-peak region.

To simulate this behavior, a continuum damage algorithm was implemented. The strain softening (damage) algorithm is based on the work of J.W. Ju and J.C. Simo. They proposed a strain-energy-based damage criteria.(12-13) The major advantage of their method is that the strain softening is uncoupled from the plasticity algorithm. The plasticity algorithm uses undamaged stresses. This means that the plasticity algorithm can be implemented and verified independently of the damage algorithm.

Figure 15. Principal stress difference versus axial strain for triaxial compression test at sigma subscript 2 equals 6.9 megapascals of WES road-base material. Graph. This graph shows one black line. The vertical axis of this graph ranges from 0 to 30 and represents principal stress differential in megapascals, while the horizontal axis ranges from 0 to 25 and represents the percent of axial strain. The line rises along the vertical axis to the point 5, where it curves away and upward until it peaks at the points of 26 on the vertical axis and 10 on the horizontal axis. From there, it gradually descends to the points of 23 on the vertical axis and 20 on the horizontal axis. From there, it drops vertically until it leaves the graph, bending backward slightly at the point of 19 on the horizontal axis.

Figure 15. Principal stress difference versus axial strain for triaxial compression test at sigma symbol2 6.9 MPa of WES road-base material.

For the damage criterion, xi symbol, we used:

 
Equation 11. The damage criterion equals negative 1 divided by K subscript I times the integral of the pressure, P bar, with respect to the plastic volumetric strain, E subscript PV.
(11)

where: P = pressure and epsilon symolpv = plastic volumetric strain

When epsilon symolpv < 0, the soil dilating. The damaged stress is found from the undamaged stresses, namely

 
Equation 12. Sigma equals sigma bar times the sum of 1 minus the isotropic damage parameter, D.
(12)

where d = isotropic damage parameter and the damage parameter is found at step j+1 as:

 
Equation 13. D subscript J plus 1 equals D subscript J if the damage criterion subscript J plus I is less than or equal to R subscript J. However, D subscript J plus 1 equals the quotient of the damage criterion subscript J plus I minus the damage criterion subscript 0 divided by the sum of alpha minus the damage criterion subscript 0 if the damage criterion subscript J plus 1 is greater than R subscript J.
(13)

Here, r j+1 is a damage threshold surface, which is updated in this manner:

 
Equation 14. A damage threshold surface, R subscript J plus 1, equals the maximum of the set R subscript J, the damage criterion subscript J plus 1, and the damage criterion subscript 0 equals R subscript 0.
(14)

Typically, the damage, d, varies from 0 to a maximum of 1. However, some soils can have residual strength that is pressure-dependent. The residual strength is represented by , an internal friction angle.

The maximum damage allowed is related to the internal friction angle of residual strength by:

 
Equation 15. The maximum damage allowed, D subscript max, equals the sine of theta minus the sine of theta subscript res, all divided by the sine of theta.
(15)

If Q res > 0, then d max, the maximum damage, will not reach 1, and the soil will have residual strength.

The damage parameter is used to reduce the effective internal stress equation 15a.

If damage parameter d is allowed to become 1, then the internal stress is zero, which for a finite element code such as LS-DYNA (explicit) causes the internal forces (element stiffness) to become zero. By not allowing the damage parameter d to become 1, this keeps a residual stiffness in the element. Therefore, by setting phi symbolres > 0, then dmax , the maximum damage, will not reach 1, and the soil will have some residual strength. If the strains continue with approximately the same behavior, the effective internal stresses will be almost constant. However, if the strains drastically increase or decrease, then the effective internal stresses can change, because , the undamaged stresses, are changing drastically

When material models include strain softening, special techniques must be used to prevent mesh sensitivity. Mesh sensitivity is the tendency of the finite element model/analysis to produce significantly different results as the element size is reduced. Mesh sensitivity occurs because softening in the model concentrates in one element. As the element size is reduced, the failure becomes localized in smaller volumes; this causes less energy to be dissipated by the softening, leading to instabilities or at least mesh-sensitive behavior.

To eliminate or reduce the effects of strain-softening mesh sensitivity, the softening parameter, alpha symbol(the strain at full damage), must be modified as the element size changes. Normally, a material property that is independent of the test specimen size is used. For many materials (e.g., metal, concrete, wood, composites), the material property used is the fracture energy. However, for soil, there seems to be no corresponding softening property that is independent of the test specimen size. Therefore, we assume a property—void formation—and use it as an input parameter. Figure 16 shows graphically the definition of the void formation parameter, Gf. The void formation parameter is the area under the softening region of the pressure-volumetric strain curve times the cube root of the element volume, V 1/3. For the linear softening model, is just the area of the shaded triangle:

Figure 16. Definition of the void formation parameter. Diagram. This diagram shows a single triangle split in half. The right half is shaded in a light blue and labeled G subscript F divided by V. The left half of the triangle is white with no shading or labeling, while the image itself rests on a plane labeled E subscript VP times alpha. A vertical plane rises from the beginning of the horizontal plane where the left corner of the triangle sits. This plane is labeled P.

Figure 16. Definition of the void formation paramater.

 
Equation 16. The void formation parameter, G subscript F, equals the cube root of the element volume, V, times the integral of the pressure, P, with respect to the volumetric strain, E subscript V. The integration is performed from volumetric strain at peak . pressure, E subscript VP, to alpha. The result of this integration equals P subscript peak times the quantity alpha minus E subscript VP times the cube root of V, all divided by 2.
(16)

where: epsilon symbolvp = volumetric strain at peak pressure

Then, alpha symbol can be found as a function of the volume of the element V:

 
Equation 17. Alpha equals 2 times G subscript F divided by K times E subscript VP times the cube root of V, all plus input parameter E subscript VP.
(17)

where: epsilon symbolvp = an input parameter

If G f is made very small relative to K epsilon symbolvpV 1/3, then the softening behavior will be brittle.

Strain-Rate Behavior

There are some experimental data that suggest that soil strength is strain-rate-dependent. (11) Based on the earlier evaluation of a steel post in soil experiments,(14) strength enhancements caused by high strain rates may not be needed for roadside safety applications. However, since the development and implementation of strain-rate-dependent effects are relatively easy and will not affect the overall efficiency if they are not used, strain-rate effects were implemented into the soil model.

The two-parameter Devaut-Lions viscoplastic update algorithm developed by Y. Murray was used. (15) This algorithm interpolates between the elastic trial stress and the inviscid stress. The inviscid stresses are on the plasticity surface eq17b, with equation 17c and equation 17d

Figure 17 shows the behavior of zeta symbol . As zeta symbol becomes 1, then the viscoplastic stress becomes the elastic trial stress.

Figure 17. Zeta versus strain rate for different parameters. Graph. This figure shows three different lines. The first is black and labeled gamma equals 1.0E minus 3N equals 2.0; the second is red and labeled gamma equals 1.0E minus 3N equals negative 0.5; and the third is green and labeled gamma equals 1.0N equals negative 2.0. The vertical axis of this graph ranges from 0.0 to 1.2 and represents zeta, while the horizontal axis of this graph ranges from 10 to the negative tenth to 10 to the sixth and represents strain rate. The green line in this graph begins its ascent at the point of 10 to the negative third on the horizontal axis, where it climbs to 0.1 on the vertical axis and spikes nearly straight upward to the point of 0.9 on the vertical axis. It then bends up to the points of 1.0 on the vertical axis and 10 to the negative one on the horizontal axis. From there, it plateaus until it stops at the point of 10 to the fifth on the horizontal axis. The red line begins its ascent at the point of 10 to the negative fourth, where it climbs to 0.05 on the vertical axis and then to 0.1 on the vertical at the point of 10 to the negative second on the horizontal axis. From there, it rises almost vertically to the points of 0.8 on the vertical axis and 10 to the 0 power on the horizontal axis, at which point the line curves upward until it stops at the points of 1.0 on the vertical axis and 10 to the fifth on the horizontal. The black line begins its ascent at the point of 10 to the negative one, where it spikes upward in a straight line vertically until it plateaus at the points of 1.0 on the vertical axis and 10 to the 0 power on the horizontal axis. The plateau carries on to the point of 10 to the first on the horizontal axis, where it runs parallel with the green line.

Figure 17. Zeta versus strain rate for different parameters.

INCORPORATION INTO LS-DYNA

The above discussion describes the equations and the form of the soil model. In this section, the implementation of the model into LS-DYNA is discussed. The intermediate equations were determined using the symbolic algebra program Mathematica ®.

First, the elastic moduli and undamaged stresses are found (see figure 18). The bulk modulus function needs the current total volumetric strain. The undamaged stresses are recovered from the damaged stresses based on the current value of the damage variable d.

Next, the elastic trial stresses are determined (see figure 19Error! Reference source not found.) and it is determined whether the state of stress is within the current yield surface. If it is, the routine bypasses the plasticity algorithm.

Figure 18. Elastic moduli and undamaged stresses. Instructions. Enter routine with: sigma subscript J, the change in E subscript J plus 1. Get undamaged stresses: sigma bar equals the quotient of sigma subscript J divided by the sum of 1 minus D subscript J. Determine elastic moduli: E subscript V subscript J plus 1 equals E subscript VJ plus the change in E subscript V subscript J plus 1. K subscript J plus 1 equals K subscript I divided by the sum of 1 plus K subscript I times D subscript 1 times N cur.

Figure 18. Elastic moduli and undamaged stresses.

Figure 19. Elastic trail stresses. Instructions. Compute elastic trail stress: change in P bar equals K minus the quotient of K subscript SK, divided by 1 plus the product of K subscript SK times D subscript 2, times N subscript cur, all times the change in E subscript V subscript J plus 1. P bar subscript J plus 1 equals P bar subscript J plus the change in P bar. Sigma bar subscript E equals sigma bar subscript J plus the change in P bar plus 2G times the sum of the change in E subscript J plus 1 minus the change in E subscript V subscript J plus 1. Determine if elastic trial stress is within the yield surface: If F of P bar, J subscript 2, theta, and E subscript EP are less than or equal to 0, go to box 6; if more than 0, go to box 3.

Figure 19. Elastic trial stresses

The total strains are split into elastic and plastic strains:

 
Equation 18. Total strain, E, equals elastic strain, E subscript E, plus plastic strain, E subscript P.
(18)

The stress increment at step j+1 is then:

 
Equation 19. The change in sigma subscript J plus 1 equals C times the sum of the change in E subscript J plus 1 minus the change in E subscript P subscript J plus 1.
(19)

To compute the stress increments, it is first necessary to determine the plastic strain increments. The latter are determined by assuming an associated flow rule (see figure 20). Note that the gradient of the yield function is taken at the stress state of the previous time step; this is consistent with the explicit algorithm used in LS-DYNA. The explicit form of the gradients is presented in appendix A. The increment in the hardening parameter phi symbol is assumed to have a form similar to the flow rule. C is the elasticity tensor. The function h will be determined later (see figure 21). For the determination of the increment of the scale parameter, the gradients of the yield function are again evaluated at step j. Once the plasticity scale parameter is determined, the stresses and history variables can be updated (see figure 21).

Figure 20. Determination of plastic strains. Instructions. The plastic strain increment is found from the associated flow rule: the change in E subscript P subscript J plus 1 equals the change in lambda times the quotient of Partial F divided by Partial sigma, straight line, subscript J. The increment in the hardening parameter phi is: the change in phi equals the change in lambda times H times E subscript EP. Employing standard techniques, the increment of the scalar parameter lambda is: the change in lambda equals Partial F divided by Partial sigma, times C times the change in E, all divided by Partial F divided by Partial sigma, times C, times Partial F divided by Partial sigma, minus Partial F divided by Partial phi, times H, straight bar.

Figure 20. Determination of plastic strains.

Figure 21. Update of stresses and history variables. Instructions. To update the stresses and plasticity variables: sigma bar subscript J plus 1 equals sigma bar subscript E minus K times the change in E subscript PV minus 2G times the change in E subscript P minus the change in E subscript PV. The change in phi equals the change in lambda times H of E subscript EP, where H of E subscript EP equals Partial phi divided by Partial E subscript EP times the square root of two-thirds times the square root of Partial F divided by Partial sigma times Partial F divided by Partial sigma. Determine if stress is on the yield surface: If F of P bar, J subscript 2, theta, and E subscript EP are less than or equal to 0, go to box 5; if more than 0, go to box 3.

Figure 21.Update of stresses and history variables.

Figure 22. Viscoplasticity update. Instructions. Determine effective strain rate: E dot subscript E equals the change in E subscript T divided by the change in T. If E dot subscript E is less than or equal to tolerance, then skip viscoplasticity update. N equals the quotient, raised to the N minus 1 divided by N power, of Lambda divided by E dot subscript E. Psi equals 1 divided by the change in T divided by N plus 1. Sigma bar subscript V-P equals 1 minus psi times sigma bar minus psi times sigma bar subscript E.

Figure 22. Viscoplasticity update.

Next, the effect of the strain rates on the strength of the material is determined. The effective strain rate must be found and used to determine the viscoplasticity scale factor. If the effective strain rate is less than 0.00001, the step is assumed to be static and the algorithm is skipped. The viscoplasticity scale factor interpolates between the elastic trial stress and the inviscid stress (see figure 22).

The strain softening algorithm transforms the effective (undamaged) stresses to the damaged stresses (see figure 23). The damage is assumed to be isotropic. The damage scale factor is determined based on total volumetric strain energy. It is desirable for strain to soften based on tensile volumetric (dilation) strain energy. In this routine, compressive strain values are considered to be positive and, therefore, the increment of volumetric strain energy is subtracted. When the element starts to dilate, the volumetric strain energy starts to increase. The initial damage threshold (a user input parameter) is input with a value close to zero.

Finally, the history variable to be plotted is updated, erosion (failure) is determined, and the routine is exited.

Once the model was verified and preliminarily validated, the material model source code was implemented into the production version of LS-DYNA in February 2002. We next checked that the soil model was correctly implemented into LS-DYNA. An independent evaluator then began the investigation of the validity of the soil model.

Figure 23. Damage update. Instructions. Update total volumetric strain energy: Psi subscript J plus 1 equals psi subscript J minus one-half times the sum of P bar subscript J plus 1 minus P bar subscript J, times the change in E subscript V subscript J plus 1. If psi subscript J plus 1 is less than R subscript J, then skip to end. D subscript J plus 1 equals psi subscript J plus 1 minus US subscript 0 divided by the sum of alpha minus psi subscript 0. Sigma equals sigma bar times the sum of 1 minus D subscript J plus 1. R subscript J plus 1 equals the maximum of R subscript J, psi subscript J plus 1, end.

Figure 23. Damage update

Table of Contents | Next

ResearchFHWA
FHWA
United States Department of Transportation - Federal Highway Administration