U.S. Department of Transportation
Federal Highway Administration
1200 New Jersey Avenue, SE
Washington, DC 20590
202-366-4000


Skip to content
Facebook iconYouTube iconTwitter iconFlickr iconLinkedInInstagram

Federal Highway Administration Research and Technology
Coordinating, Developing, and Delivering Highway Transportation Innovations

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 2. USER'S MANUAL

The user's manual was written as the model was being implemented, verified, and validated. The user's manual consists of a user input guide (much like material model sections in the LS-DYNA user's manual); a brief theory manual (LSTC theory manual), which is a condensed version of the first section of this report; and a discussion of the use of the model. Both manuals will be added to an updated LS-DYNA manual. The user's manual addresses the basics of the model, input parameters, and basic equations.

Table 2 contains a brief description of the user input variables for the soil model, along with the corresponding symbols used in the LSTC theory manual. The bold text is the LSTC theory manual symbol, which is typically followed by a brief description and then the user input value symbol. The parameters that need to be specified are dependent on the soil and the specific application.

Table 2. Input parameters for soil model.

Elastic and Soil Characteristics

K (bulk modulus or nonporous bulk modulus if pore-water effects are used, K)

G (shear modulus, G)

gamma symbol sp(specific gravity, Spgrav)

mc (moisture content, 0.0-1.00, Mcont)

rho symbol (density of soil, RO)

Plasticity

phi symbol (friction angle, radians, Phimax)

c (cohesion, units of stress, Coh)

ahyp (coefficient for modified Drucker-Prager surface, units of stress, Ahyp)

e (eccentricity parameter for third invariant effects, Eccen)

Pore-Water Effects

D1 (parameter for pore-water effects onbulk modulus, Pwd1)

Ksk (skeleton bulk modulus pore-water parameter, PwKsk)

D2 (parameter for pore-water effects on effective pressure, Pwd2)

Strain Hardening

An (strain hardening, percent of phimax where nonlinear effects start, An)

Et (strain hardening, amount of nonlinear effects, Et)

Strain Softening

xi symbol0(volumetric strain at initial damage threshold, Dint)

Gf (void formation energy, Vdfm)

phi symbol(minimum internal friction angle used for residual strength, radians, Phires)

Strength Enhancement Caused by Strain-rate Effects

gamma symbol (viscoplasticity parameter, strain-rate-enhanced strength, Gammar)

n (viscoplasticity parameter, strain-rate-enhanced strength, V n)

Element Deletion

Damlev: Level of damage that will cause element deletion (0.0-1.0)

Epsmax: Maximum principal failure strain

Miscellaneous

Nplot: Element plotting variable to put into effective plastic strain variable

Rhowat: Density of water in model units, used to determine air void strain (saturation)

Itermax: Maximum number of iterations used in plasticity iterations

USER INPUT GUIDE

FHWA Soil Material Model Input

*MAT_FHWA_SOIL_OPTION

  1. Available options include:
    1. NEBRASKA
    2. < BLANK >
  2. such that the keyword cards appear as
    1. *MAT_FHWA_SOIL
    2. *MAT_FHWA_SOIL_NEBRASKA

This is material type 147. This is an isotropic material with damage and is available for solid elements in LS-DYNA. The model has a modified Mohr-Coulomb surface to determine the pressure-dependent peak shear strength. It was developed for applications involving road-base soils.

*MAT_FHWA_SOIL_NEBRASKA

It is an option to use the default properties determined for soils used at the University of Nebraska at Lincoln. The default units used for this material are millimeter (mm), millisecond (ms), and kilogram (kg). If different units are desired, the conversion factors must be input.

Card Format

Card 1 1 2 3 4 5 6 7 8
Variable MID FCTIM FCTMAS FCTLEN        
Type I F F F        
Default None 1.0 1.0 1.0        

 

Variable Description
MID Material identification (a unique number has to be chosen)
FCTIM Factor by which to multiply milliseconds to get desired time units
FCTMAS Factor by which to multiply kilograms to get desired mass units
FCTLEN Factor by which to multiply millimeters to get desired length units

*MAT_NCHRP_SOIL_blank

Define the following cards:

Card Format

Card 1 1 2 3 4 5 6 7 8
Variable MID RO Nplot Spgrav Rhowat Vn Gammar Itermax
Type I F I F F F F I
Default None None 1 None 1 0 0 1
Card 2 1 2 3 4 5 6 7 8
Variable K G Phimax Ahyp Coh Eccen An Et
Type F F F F F F F F
Default None None None None None None None None
Card 3 1 2 3 4 5 6 7 8
Variable Mcont Pwd1 PwKsk Pwd2 Phires Dint Vdfm Damlev
Type F F F F F F F F
Default None None None None 0 None None None
Card 4 1 2 3 4 5 6 7 8
Variable Epsmax              
Type F              
Default None              

 

Variable Description
MID Material identification (a unique number has to be chosen)
RO Mass density
Nplot

Plotting options:

1 Effective strain
2 Damage criterion threshold
3 Damage (diso)
4 Current damage criterion
5 Not used
6 Current friction angle (phi)

Spgrav Specific gravity of soil used to get porosity
Rhowat Density of water in model units, used to determine air void strain (saturation)
Vn Viscoplasticity parameter (strain-rate-enhanced strength)
Gammar Viscoplasticity parameter (strain-rate-enhanced strength)
Itermax Maximum number of plasticity iterations (default 1)
K Initial bulk modulus or nonporous bulk modulus if pore-water effects are used (non-zero)
G Shear modulus (non-zero)
Phimax Peak shear strength angle (friction angle) (radians)
Ahyp Coefficient for modified Drucker-Prager surface
Coh Cohesion, shear strength at zero confinement (overburden)
Eccen Eccentricity parameter for third invariant effects
An Strain hardening percent of phimax where nonlinear effects start
Et Strain hardening amount of nonlinear effects
Mcont Moisture content of soil (determines amount of air voids) (0-1.00)
Pwd1 Parameter for pore-water effects on bulk modulus
PwKsk Skeleton bulk modulus, pore-water parameter, set to zero to eliminate effects
Pwd2 Parameter for pore-water effects on effective pressure (confinement)
Phires Minimum internal friction angle (radians) (residual shear strength)
Dint Volumetric strain at initial damage threshold (xi symbol0 )
Vdfm Void formation energy (like fracture energy)
Damlev Level of damage that will cause element deletion (0.0-1.0)
Epsmax Maximum principal failure strain

THEORY MANUAL

MAT_FHWA_SOIL

A brief discussion of the FHWA soil model is given. The elastic properties of the soil are isotropic. The implementation of the modified Mohr-Coulomb plasticity surface is based on the work of Abbo and Sloan. (9) The model is extended to include excess pore-water effects, strain softening, kinematic hardening, strain-rate effects, and element deletion.

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 and it is perpendicular to the pressure axis. The yield surface is given as

 
Equation 20. F equals negative P times the sine of phi plus the square root of the following: J subscript 2 times K theta squared, plus parameter AHYP, squared, times the sine, squared, of theta. Then subtract C times the cosine of theta, and this equals 0.
(20)

where:

P=pressure

phi symbol = internal friction angle

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

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

c = amount of cohesion, eq20a

J3 = third invariant of the stress deviator

ahyp = parameter for determining how close to the standard Mohr-Coulomb yield surface the modified surface is fitted

If ahyp is input as zero, the standard Mohr-Coulomb surface is recovered. The input parameter ahyp should be set close to zero, based on numerical considerations, but always less than . It is best not to set the cohesion, c, cot phi symbol to very small values since this causes excessive iterations in the plasticity routines.

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

 
Equation 21. 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.
(21)

where:

J3 = third invariant of the stress deviator

e = material parameter describing the ratio of triaxial extension strength to triaxial compression strength

If e is set to 1, then a circular cone surface is formed. If e is set to 0.55, then a triangular surface is formed. K (theta symbol) is defined for 0.5. < e less than or equal to1.0.

To simulate nonlinear strain hardening behavior, the friction angle is increased as a function of the effective plastic strain:

 
Equation 22. The change in friction angle, phi, equals the input parameter, E subscript T, times the sum of 1 minus phi minus phi subscript init divided by a fraction of the peak strength internal friction angle where nonlinear behavior begins, A subscript N, times phi subscript max, times the change in effective plastic strain, E subscript eff plas.
(22)

where:

xi symbol eff plas = effective plastic strain

An = fraction of the peak strength internal friction angle where nonlinear behavior begins, 0 < An less than or equal to - 1

The input parameter E1 determines the rate of the nonlinear hardening. If there is no strain hardening, then phi symbol max = phi symbol init= phi symbol .

To simulate the effects of moisture and air voids, including excess pore-water pressure, both the elastic and plastic behaviors can be modified. The bulk modulus is:

 
Equation 23. 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.
(23)

where:

Ki = nonporous bulk modulus

n cur = current porosity =equation 23a

w = volumetric strain corresponding to the volumne of the air voids

= n (1-S)

epsilon symbol v = total volumetric strain

D 1 = 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 24 shows the effect of the D 1 parameter on the pressure-volumetric strain relationship (bulk modulus). The bulk modulus will always be a monotonically increasing value, that is:

 
Equation 24. 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.
(24)

Note that the model is following the standard practice of assuming that compressive stresses and strains are positive. If the input parameter is zero, then the standard linear elastic bulk modulus behavior is used.

Figure 24. Pressure versus volumetric strain showing the effects of the D subscript 1 parameter. Graph. This figure shows five separate lines. The first is a solid black line labeled K subscript J equals 0.468, D subscript 1 equals 0.68. The second is a broken red line labeled K subscript J equals 0.468, D subscript 1 equals 5.0. The third, a dotted orange line, is labeled K subscript J equals 2.5, D subscript 1 equals 5.0. The fourth is a broken and dotted green line labeled K subscript J equals 11, D subscript 1 equals 5.0. The fifth is a broken blue line labeled K subscript J equals 15, 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 volumetric strain. Both the green and blue lines rise gradually from the start to the points of 0.1 on the vertical axis and 0.11 on the horizontal axis, where they both ascend diagonally; the blue line leaves the graph at the point just below 1.5 on the vertical axis, and the green leaves the graph at the point of 1.2 on the vertical axis. The dotted orange line ascends very gradually from the start, curving upward slightly until it leaves the graph at the point just above 0.3 on the vertical axis. Both the black and broken red lines barely ascend across the graph as they leave the graph at the point of 0.1 on the vertical axis.

Figure 24. Pressure versus volumetric strain showing the effects of the D subscript 1 parameter.

If D1 is not set to zero, the bulk modulus input should be the fully collapsed bulk modulus.

To simulate the loss of shear strength caused by excess pore-water effects, the model uses a standard soil mechanics technique(11) of reducing the total pressure, P, by the excess pore-water pressure, u, to get an "effective pressure," P' :

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

Figure 25. Effects on pressure caused by pore-water pressure. Graph.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. This figure is the same as figure 11.

Figure 25. Effects on pressure caused by pore-water pressure.

Figure 25 shows how pore-water pressure affects the algorithm for the plasticity surface. The excess pore-water pressure reduces the total pressure, which lowers the shear strength, square root of J subscript 2. Significant excess pore-water pressure can cause the effective pressure to become zero. To calculate the pore-water pressure, u, the model uses an equation similar to the equation used for the moisture effects on the bulk modulus:

 
Equation 26. U equals the total volumetric strain, E subscript V, times the quotient of the bulk modulus for soil without air voids (skeletal bulk modulus), K subscript SK, divided by the sum of 1 plus the product of K subscript SK times the material constant controlling the pore-water pressure before the air voids are collapsed, D subscript 2, times the current porosity, N subscript cur.
(26)

where:

Ksk = bulk modulus for soil without air voids (skeletal bulk modulus)

ncur = current porosity = Max[0, (w-ev)]

w = volumetric strain corresponding to the volume of air voids

=n(1- S)

epsilon symbol v = total volumetric strain

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

are collapsed D2greater than or equal to 0

n = porosity of the soil = e over 1 plus e

e = void ratio = gamma subscript sp (1+m subscript c) over rho minus 1

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

rho, gamma subscript sp, m subscript c= soil density, specific gravity, and moisture content, respectively

The increment pore-water pressure is zero if the incremental mean strain is negative (tensile).

Figure 26 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 strain behavior

The parameter D2 is found from Skempton pore-water pressure parameter B, where B is defined as:(7)

 
Equation 27. 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.
(27)

 

 
Equation 28. 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.
(28)

Figure 26. Effects of D subscript 2 and K subscript SK parameters on pore-water pressure. Graph. This graph shows three lines; the first is black and labeled K subscript SK equals 15, D subscript 2 equals 1500; the second is red and labeled K subscript SK equals 15, D subscript 2 equals 15; and the third is green and labeled K subscript SK equals 30, 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. From that point, they both 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 26. Effects of D subscript 2 and K subscript SK parameters on pore-water pressure.

To simulate strain softening behavior, the FHWA soil model uses a continuum damage algorithm. The strain-based damage algorithm is based on the work of J.W. Ju and J.C. Simo. They proposed a strain-based damage criterionthat is uncoupled from the plasticity algorithm.(12,13)

For the damage criterion , equation28a, where equation 28b = pressure and epsilon symbolpv = plastic volumetric strain, the damaged stress is found from the undamaged stresses:

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

where: d = isotropic damage parameter (diso)

The damage parameter is found at step j + 1 as:

 
Equation 30. 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.
(30)

where

r subscrip j plus 1 = damage threshold surface

r subscrip j plus 1 equals max (r subscript j, xi subscript j plus 1)

xi subscript zero equals r subscript zero (Dint)

The mesh-sensitivity parameter, alpha symbol, is described below.

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

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

 
Equation 31. 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.
(31)

 

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

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 the softening in the model is concentrated in one element. As the element size is reduced, the failure becomes localized in smaller volumes, which causes less energy to be dissipated by the softening. This can lead 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. The FHWA soil model uses an input parameter, "void formation," G f , that is like the fracture energy material property for metals. 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 :

 
Equation 32. 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.
(32)

with xi symbol0 as the volumetric strain at peak pressure (strain at initial damage (Dint)). Then, alpha symbol can be found as a function of the volume of the element V:

 
Equation 33. 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.
(33)

If Gf is made very small relative to K zeta , then the softening behavior will be brittle.

Strain-rate-enhanced strength is simulated by a two-parameter Devaut-Lions viscoplastic update algorithm developed by Y. Murray.(15) This algorithm interpolates between the elastic trial stress (beyond the plasticity surface) and the inviscid stress. The inviscid stresses (sigma symbol) are on the plasticity surface equation 33b, with equation 33cand equation 33d

As zeta symbol becomes 1, then the viscoplastic stress becomes the elastic trial stress. Setting the input value gamma symbolr = 0 (gamma) eliminates any strain-rate-enhanced strength effects.

The model allows element deletion if needed. As the strain softening (damage) increases, the effective stiffness of the element can become very small, causing severe element distortion and "hourglassing." The element can be "deleted" to remedy this behavior. There are two input parameters that affect the point of element deletion. Damlev is the damage threshold where element deletion will be considered. Epsmax is the maximum principal strain where the element will be deleted. Both d greater than or equal to Damlev and epsilon symbol pr max Epsmax are required for element deletion to occur. If Damlev is set to zero, there is no element deletion. Care must be taken when employing element deletion to ensure that the internal forces are very small (element stiffness is zero) or significant errors may be introduced into the analysis.

MAT_FHWA_SOIL_NEBRASKA

This option gives the soil parameters that were used to validate the material model with experiments performed at the University of Nebraska at Lincoln. The units of these default inputs are milliseconds, kilograms, and millimeters. There are no required input parameters except for material ID (MID). If different units are desired, the appropriate unit conversion factors can be input.

DISCUSSION OF SOIL MODEL USE

Material models for geomaterials (soils, concrete, rock, etc.) tend to be complex. The determination of the input parameters for the models is complicated. In addition, modeling different loading conditions and accurate simulation of boundary conditions add to the complexity involved in using these material models.

There are two methods that are typically used to determine the material input variables for soils. The most accurate method is to perform laboratory tests that include both triaxial compression and uniaxial strain tests. These tests can be used to determine the elastic moduli, yield surface parameters, and softening parameters. Typically, these tests use drained soil conditions. Laboratory tests with undrained soil conditions can be used to determine the pore-water effects.

A second method is to use full-scale testing of the specific application (e.g., a bogie impacting a steel post) to fit the parameters in a trial-and-error method. This method requires more time by the analyst. Since the soil model is nonlinear, there may not be a set of unique input parameters that can be determined.

Compaction of the soil is typically used to remove some of the air voids that exist in disturbed soils. However, the density, pore-water effects, stiffness, and strength are also changed upon compacting the soil. To simulate compaction in highway safety applications where the soil is exposed, we recommend that the values for the soil density, pore-water effects, stiffness, and strength be modified. Applying pressure to the ground surface to account for the effects of compaction is a less accurate method that will incorrectly simulate how the soil is deformed at the surface.

In full-scale testing or applications, the soil typically extends to infinity. Analyses typically do not extend to infinity, so some type of boundary condition must be applied to the exterior surfaces of a soil analysis model (except for soil surfaces exposed to atmospheric pressure). Standard boundaries reflect dynamic disturbances (stress waves), which does not happen in the real applications. Such reflections can cause serious contamination of the analysis results. Exterior boundaries for analyses involving soil need a nonreflecting boundary. A partial nonreflecting boundary exists in LS-DYNA. This boundary is an impedance-matching boundary, which is only good for high-frequency (highly transient) behavior. At this time, there is no nonreflecting boundary that matches both low- (quasi-static) and high- (highly transient) frequency behaviors. Also, only linear behavior is assumed. Thus, to use the current nonreflecting boundary, the material near the boundaries must only behave linearly. Also, the nonreflecting boundary should only experience high-frequency behavior.

Previous | Table of Contents| Next

Federal Highway Administration | 1200 New Jersey Avenue, SE | Washington, DC 20590 | 202-366-4000
Turner-Fairbank Highway Research Center | 6300 Georgetown Pike | McLean, VA | 22101