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-094
Date: November 2004

Evaluation of LS-DYNA Soil Material Model 147

PDF Version (2.87 MB)

PDF files can be viewed with the Acrobat® Reader®

CHAPTER 4. MATERIAL INPUT PARAMETER STUDY

One of the most difficult tasks associated with finite element modeling is the selection of appropriate material properties to accurately represent physical behavior. In many cases, this has led to the arbitrary tweaking of models by matching simulation results to a known physical test, sometimes without regard to the reasonableness of the material input variables.

Of course, it is desirable to be able to determine material input variables directly through physical testing. In order to determine the appropriate material input variables, the physical implications must be completely understood. Thus, a complete discussion of the variables and, when applicable, their physical meaning follows.

Many of the developer's initial soil parameters were based on triaxial testing performed by the U.S. Army Corps of Engineers (USACE) at the Waterways Experiment Station ( WES). These data are contained in the Material Property Query ( MPQ) database and are available only to documented government contractors.(9)

The objectives of this chapter are to: (1) provide an understanding of the soil parameters, (2) provide guidelines for appropriate parameter values, and (3) determine recommended parameter values for NCHRP 350 soil used at the user's facility. The following discussions attempt to provide an engineering interpretation of the material such that appropriate parameter values can be determined more readily.

Soil Density, RHO

RHO is the density of the soil material. It has the units of mass/volume (for this case, kilograms per cubic millimeter (kg/mm3)). The developer's initial models included a density of 2.35E-6 kg/mm3, almost the maximum dry density of the soil (rho symbolmax = 2.37E-6 kg/mm3) found using a Modified Procter test. The soil tested by USACE had a density of rho symbolmax = 2.27E-6 kg/mm3, which is more appropriate for densely compacted aggregate.

NCHRP 350 strong soil routinely has densities of about 2.114E-6 kg/mm3. However, densities between 2.082E-6 kg/mm3 and 2.242E-6 kg/mm3 are reasonable.

Mass Verification

From the developer's hydraulic tension test file, hydten.k, a 1-mm cube of soil material was modeled. Calculating the mass of the sample from the density and volume:

 
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.
(1)

This was verified in the D3HSP output file from LS-DYNA. The FHWA soil model produces the correct mass from an inputted density. It is critical to verify the behavior of the density function of the soil model because significant inertial effects involving the soil-post interaction have been documented.(5)

Plotting Options, NPLOT

Nplot allows the plotting of information concerning the soil model. This includes information on the effective strain, damage criterion threshold, isentropic damage, current damage criterion, and current friction angle. The Nplot options are as follows:

  1. Effective strain.
  2. Damage criterion threshold.
  3. Isentropic damage parameter, (diso).
  4. Current damage criterion.
  5. Not implemented.
  6. Current friction angle.

For post-processing using LS-POST, the value specified in Nplot is stored as the Plastic Strain variable for plotting the fringe component stresses. When post-processing, the information stored for plastic strain will not contain the plastic strain data, but rather the data specified under Plotting Options, Nplot.

Specific Gravity, SPGRAV

The specific gravity of soil solids, usually designated by Gs, is the ratio of the soil solids to the density of water. While it is possible to have a range of values from 2.2 to 3.5, most soils have a specific gravity from 2.60 to 2.80. Any values outside of this latter range should be viewed skeptically, and the soil should be retested to verify the value. Where specific values are not available, the following can be assumed for local soils:(10)

Sand and gravels: Gs= 2.65

Silts and clay: Gs = 2.78

It should be noted that the specific gravity, Gs , is not the bulk specific gravity, Gblk , of the material, but rather the specific gravity of only the soil solids. This is to say that air voids internal to and between soil particles are not considered when calculating the specific gravity. Since NCHRP 350 strong soil is designated as a gravel, a specific gravity of Gs  = 2.65 is appropriate. However, laboratory testing should be performed to accurately identify the exact specific gravity associated with the Nebraska crushed limestone used at the user's facility.

In the FHWA soil model, the specific gravity of the soil is used to calculate the porosity. Soil porosity is the ratio of the total volume of voids, Vv , to the total volume, Vt , of a sediment.

A discrepancy exists in the calculation of the Spgrav input variable as it is read by LS‑ DYNA. In the baseline model, a value of Spgrav = 2.79 was used. This value came from the USACE WES testing. In the D3HSP output file, produced by LS-DYNA, SPGRAV was reported as 0.1854, a variation of a factor of 15. These values are shown in bold in Table 3. The other values of the input deck matched exactly, as shown in table 3. The D3HSP output from LS-DYNA is shown in table 4. This discrepancy was observed through the normal verification of input decks with output parameters.

The ramifications of this discrepancy are not fully understood. It is impossible to determine whether the material model is reading the correct values and merely outputting incorrect values only to the D3HSP file or whether the material model is, in fact, altering Spgrav incorrectly. This can only be verified through a careful examination of the source code.

Table 3. Comparison between LS-DYNA input deck and D3HSP output reveals the Spgrav discrepancy (shown in bold).

Parameter

RO

Nplot

Spgrav

Rhowat

Vn

Gammar

Input Deck

2.350E-6

3

2.79

1.0E-6

1.1

0.0

D3HSP File

2.350E-6

3

0.1854

1.0E-6

1.1

0.0

 

Parameter

Itermax

K

G

Phimax

Ahyp

Coh

Input Deck

10

0.465

0.186

1.1

1.0E-7

6.2E-6

D3HSP File

10

0.465

0.186

1.1

1.0E-7

6.2E-6

 

Parameter

Eccen

An

Et

Mcont

Pwd1

PwKsk

Input Deck

0.7

0.0

0.0

0.034

0.0

0.0

D3HSP File

0.7

0.0

0.0

0.034

0.0

0.0

 

Parameter

Pwd2

Phires

Dint

Vdfm

Damlev

Epsmax

Input Deck

0.0

0.0E-0

2.5E-3

5.0

1.00

1.00

D3HSP File

0.0

0.0E-0

2.5E-3

5.0

1.00

1.00


Table 4. D3HSP output file (truncated).

SOIL

Part ID ........................ 1

material type .............. 147

equation-of-state type ..... 0

hourglass type ............. 1

bulk viscosity type ........ 1

den .............................. = 2.35000E-06

hourglass coefficient ............ = 1.00000E-01

quadratic bulk viscosity ......... = 1.50000E+00

linear bulk viscosity ............ = 6.00000E-02

element type ..................... = 0

eq.0: 8-node solid element

eq.1: 2-node beam/truss element

eq.2: 4-node membrane/shell element

eq.3: 8-node thick shell element

flag for bulk viscosity in shells = 0

flag for RBDOUT/MATSUM output .... = 0

eq.0: RBDOUT and MATSUM

eq.1: RBDOUT only

eq.2: MATSUM only

eq.3: no output

static coefficient of friction ... = 0.00000E+00

kinetic coefficient of friction .. = 0.00000E+00

exponential decay coefficient .... = 0.00000E+00

viscous friction coefficient ..... = 0.00000E+00

optional contact thickness ....... = 0.00000E+00

optional thickness scale factor .. = 0.00000E+00

local penalty scale factor ....... = 0.00000E+00

flag for adaptive remeshing ...... = 0

eq.0: inactive

eq.1: h-adaptive only

eq.2: r-adaptive only

rayleigh damping coefficient ..... = 0.00000E+00

nplot ............................ = 3.0000E+00

spgrav ........................... = 1.8540E-01

rhowat ........................... = 1.0000E-06

vn ............................... = 1.1000E+00

gammar ........................... = 0.0000E+00

itermax .......................... = 1.0000E+01

bulk modulus ..................... = 4.6500E-01

shear modulus .................... = 1.8600E-01

phimax ........................... = 1.1000E+00

ahyp ............................. = 1.0000E-07

coh .............................. = 6.2000E-06

eccen ............................ = 7.0000E-01

an ............................... = 0.0000E+00

et ............................... = 0.0000E+00

mcont ............................ = 3.4000E-02

pwd1 ............................. = 0.0000E+00

pwksk ............................ = 0.0000E+00

pwd2 ............................. = 0.0000E+00

phires ........................... = 0.0000E+00

dint ............................. = 2.5000E-03

vdfm ............................. = 5.0000E+00

damlev ........................... = 1.0000E+00

epsmax ........................... = 1.0000E+00

Density of Water, RHOWAT

RHOWAT is the density of water (1.0 x 10-6 kg/mm3 (62.4 lb/ft3)). This is used to determine the air void strain when calculating pore-water effects.

Viscosity Parameters, V n and GAMMAR

GAMMAR (gamma symbolr)and Vn are viscosity parameters used to develop the strain-rate-enhanced strength of the material model. The algorithm interpolates between the elastic trial stress (beyond the yield surface) and the inviscid stress (stresses where the material viscosity effects are so small that they can be neglected). The inviscid stresses are on the yield surface. In equation form, this is written:

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

where:

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

and

 
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.
(4)

Setting GAMMAR to 0.0 eliminates any strain-rate-enhanced strength effects, regardless of any values that remain for Vn . Additional work must be performed to determine the appropriate values for these strain-rate parameters.

Bulk Modulus, K

The bulk modulus, K, is an elastic constant that reflects the resistance of the material to an overall gain or loss of volume under conditions of hydrostatic stress. If the hydrostatic stress increases, then the volume will decrease and the volume change will be negative. If the hydrostatic stress decreases, then the volume will increase.

The relationship between the elastic modulus, E, and the bulk modulus, K, is:

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

The determination of the elastic modulus, E, is difficult because volume changes in the soil require the precise determination of pore pressures. Moreover, the value of Poisson's ratio may influence the results. Poisson's ratio, nu symbol, has the following approximate values in soils:(11)

nu symbol = 0.5 (saturated impervious soils)

nu symbol = 0.25 (pervious coarse materials)

For NCHRP 350 strong soil, it is appropriate to use nu symbol = 0.25. On this basis, the equation above simplifies to:

 

 
E equals 1.5K.
(6)

Rewritten,

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

In 1975, Penman performed tests on gravel, finding an elastic modulus of E = 15.8 MPa.(12) Penman also found a Poisson's ratio of nu symbol = 0.27. This corresponds well with known data as shown above. This would imply a bulk modulus of K = 10.5 MPa. The initial developer's models used a value of K = 465 MPa. As mentioned previously, the developer's revised suggested value for the bulk modulus was K = 3.25 MPa.

Shear Modulus, G

Just as the modulus of elasticity, E, is a measure of the relationship of the stress to the strain below the proportional limit, the shear modulus of elasticity, G, relates shear stress to shear strain. The shear modulus is also referred to as one of the two Lamé constants, G and lambda symbol. Using conventional engineering mechanics, the shear modulus can be expressed as a function of the modulus of elasticity, E, and Poisson's ratio, nu symbol :

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

Using the values found by Penman,(12) a shear modulus of G = 6.22 MPa is found. The initial developer's models used a value of G = 186 MPa. As mentioned previously, the developer's revised suggested value for the shear modulus was 1.3 MPa.

Poisson's Ratio

In 1987, Trautmann and Kulhawy found general ranges of Poisson's ratio for granular soils.(13) These values are shown in table 5. It is important to ensure that appropriate values of Poisson's ratio are used in the FHWA soil material model.

Table 5. General range of Poisson's ratio for granular soils.

Soil Type

Range of Poisson's Ratio

Loose Sand

0.20-0.40

Medium Dense Sand

0.25-0.40

Dense Sand

0.30-0.45

Silty Sand

0.20-0.40

Sand and Gravel

0.15-0.35

NCHRP 350 strong soil is most similar to the "sand and gravel" listed in table 5. It is believed that the most reasonable values would lie at approximately nu symbol = 0.25, as noted previously.

The bulk modulus, K; the shear modulus, G; Poisson's ratio, nu symbol; and the modulus of elasticity, E, are all interrelated. These relationships can be used to determine the value of Poisson's ratio, nu symbol, that the developer used in both the initial models and the subsequent recommended values. This is accomplished by solving for the elastic modulus, E, and combining equations 5 and 8 as follows:

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

Rearranging, it is seen that:

 
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.
(10)

Substituting with nu symbol= 0.25 yields:

 
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.
(11)

Solving through, this yields:

 
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.
(12)

Equation 9 can also be solved for Poisson's ratio as a function of the bulk and shear modulus. For the developer's initial values of K = 465 MPa and G = 186 MPa, a Poisson's ratio of 0.32 can be calculated. Similarly, with the developer's subsequent recommendations via email, with K = 3.25 MPa and G = 1.3 MPa, a Poisson's ratio of 0.32 can also be calculated. This is not an unreasonable value.

It is important that the end user understand the relationship between the shear modulus and the bulk modulus. In order to maintain consistency with the laws of physics and conventional engineering mechanics, reasonable and appropriate values of Poisson's ratio must exist.

Appropriate Values for Bulk and Shear Moduli

It is critical to associate material model parameters with physical test results that can be performed in the field. The modulus of elasticity has been correlated to the standard penetration number, N, and also the cone penetration resistance, qc, by various investigators.(14) These values and their corresponding bulk and shear moduli are shown in table 6.

Table 6. General range of bulk and shear moduli for nu symbol = 0.25.

Soil Type

Modulus of Elasticity, E (MPa)

Bulk Modulus, K (MPa)

Shear Modulus, G (MPa)

Loose Sand

10.35-24.15

6.90-16.10

4.14-9.66

Medium Dense Sand

10.35-17.25

6.90-11.50

4.14-6.90

Dense Sand

17.25-27.60

11.50-18.40

6.90-11.04

Silty Sand

34.50-55.20

23.00-36.80

13.80-22.08

Sand and Gravel

69.00-172.50

46.00-115.00

27.60-69.00

Nuclear densometer readings from field testing of roadbed materials have given values for the modulus of elasticity between 26.2 MPa and 193 MPa.(15) These values correspond well to values found by Das (10) and do not seem unreasonable compared to Penman.(12)

With these values, it would seem that selecting the median value for a "sand and gravel" soil would be appropriate. Median values would be K = 80.5 MPa and G = 48.3 MPa for the bulk and shear moduli, respectively. However, these values produced too stiff of a soil response in the direct shear test simulation, relative to physical testing. Using Penman's recommended values of 10.5 MPa and 6.22 MPa for the bulk and shear moduli, respectively, the direct shear test simulation still produced an unreasonably stiff soil response. When the bulk and shear moduli were adjusted to the developer's recommended values of 3.25 MPa and 1.3 MPa, respectively, the model appeared to provide a more reasonable prediction of soil stiffness.

Angle of Internal Friction, PHIMAX, aND Cohesion, COH

In 1900, Mohr presented a theory for rupture in materials that contended that a material fails because of a critical combination of normal stress and shearing stress, not from either maximum normal or shear stress alone.(16) For most soil mechanics problems, it is sufficient to approximate the shear stress on the failure plane as a linear function of the normal stress.(17) Hence, the linear function can be written as:

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

where:

c = Cohesion

sigma symbol = Normal stress

phi symbol = Angle of internal friction

The preceding relationship is called the Mohr-Coulomb failure criteria. Values for cohesion, c, and the angle of internal friction, phi symbol, can be determined through direct shear testing or triaxial compression tests. A graphical representation of the Mohr-Coulomb failure criteria is shown in figure 8.

View Alternative Text

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

Cohesion

When soil is removed from a bed of dry or completely immersed sand, the material at the sides of the excavation slides toward the bottom. This behavior indicates the complete absence of a bond between the individual sand particles. The sliding material does not come to rest until the angle of inclination of the slopes becomes equal to a certain angle known as the angle of repose (the angle of friction, which is the angle of repose in a cohesionless material such as sand).(18) The angle of repose of dry sand is independent of the height of the slope.

On the other hand, a trench 6.1 to 9.1 meters (m) (20 to 30 feet (ft)) deep with unsupported vertical sides can be excavated in stiff plastic clay. This fact indicates the existence of a firm bond between the clay particles. However, as soon as the depth of the trench exceeds a certain value, dependant on the intensity of the bond between the clay particles, the sides of the trench fail and the slope of the debris that covers the bottom of the cut after failure is far from vertical. The bond between the soil particles is called cohesion. No definite angle of repose can be assigned to a soil with cohesion, because the steepest slope at which such a soil can stand decreases with the increasing height of the slope. However, even sand, if it is moist, has apparent cohesion because of matrix suction between the grains of sand.

Coulomb's paper (1773) quoted Musschenbroek's idea (1729) that for construction materials, tensile strength (adhesion) is about equal to shear strength with no overburden (cohesion).(19) Coulomb found that, for physical tests on 1290-square millimeter (mm2) (2-square inch (inch2)) cross section specimens of limestone, the tensile failure load was 1.91 kN (430 pounds force (lbf)) and the shear failure load was 1.96 kN (440 lbf). These and other tests on brick and wood confirmed Musschenbroek's idea. Hence, if adhesion is known to be small or negligible for some material, then the cohesion of that material must also be taken to be zero.

For Coulomb, the fracture of intact bodies of undisturbed soil and rock involved both friction and cohesion; the flow of ground that has been broken up and is newly disturbed does not involve cohesion. Placement of fill behind a wall involves breaking up ground with picks, shoveling soil or broken rock into barrows, wheeling it to the site, and tipping it behind the wall. Coulomb states three times in his design calculations for such fill that there is no adhesion in newly disturbed soil.

For soil consistent with crushed limestone, such as NCHRP 350 strong soil, cohesion is, by definition, zero, because it is a cohesionless soil. This can be verified by relating the adhesion of the soil to the cohesion―tensile tests on the strong soil would show that there is no adhesion since there is no bond between the individual pieces of aggregate.

Rather than introduce an apparent cohesion of soil that is itself a function of strain, it is better to characterize peak strength as the sum of the critical state angle of repose plus a dilation angle. This interlocking strain rate depends on effective pressure and relative density.(20)

However, direct shear testing shows the presence of some cohesion―the failure envelope has a positive value as it passes through the pressure axis. It is noted, however, that the strong soil has no cohesion. This is caused by the dilation of the soil specimen during testing. The work caused by this effect, designated as interlocking by Taylor, is a distinctly different phenomenon than the work caused by friction.(21) While interlocking can be treated, in a general sense, as particle cohesion, it is important to differentiate between the two physical phenomena.

The concept of Taylor's aggregate interlock explains how soil can exhibit apparent cohesion as it flows, exhibit apparent cohesion at peak strength, and still satisfy Coulomb's law. From a mechanics standpoint, it would be best to represent the material with aggregate interlock; from a simplicity standpoint, it may be best to have artificial cohesion to represent the aggregate interlock.

Cohesion in the soil model was varied in a parameter study to determine the appropriate values. The results are provided in table 7 and figure 9. As a reminder, the results throughout this chapter are limited to the point where the force required to shear the soil reached a minimum force (referred to as "Valley" in the tables). Simulations continued after that point; however, forces began to rise unrealistically (as discussed in chapter 3).

It did not appear that significant differences existed between "small" (i.e., any value less than 6.2e-7 gigapascals (GPa)) and smaller values, as shown in table 7. The model successfully ran even with cohesion set to zero; however, the plasticity routines were limited by the parameter Itermax, the maximum number of iterations that allow the plasticity algorithm to converge.

Table 7. Cohesion parameter study.

Test Identifier

Angle of Internal Friction, phi symbol

Cohesion, Coh (GPa)

Itermax

CPU* Time (seconds (s))

Peak

Valley

Disp (mm)

Force (kN)

Internal Energy

Disp (mm)

Force (kN)

Internal Energy

T4

0.3°

+6.2E-7

10

2,183

18.8

805

8529

24.5

473

13,044

T5

0.3°

+6.2E-6

10

2,177

18.8

805

8529

24.5

473

13,044

T33

63°

+6.2E-5

10

1,962

19.4

741

7257

25.7

395

11,907

Baseline

63°

+6.2E-6

10

2,644

18.2

718

6795

24.6

385

11,450

T6

63°

+6.2E-7

10

3,608

17.1

689

6358

23.9

387

11,289

T7

63°

+6.2E-8

10

4,038

17.0

708

6450

23.4

389

11,049

T8

63°

+6.2E-9

10

4,018

16.9

703

6335

24.4

396

11,975

T9

63°

+6.2E-10

10

3,775

16.9

708

6350

23.3

376

10,865

T10

63°

+0.0E+0

10

3,727

16.9

704

6351

24.0

393

11,609

T11

63°

+0.0E+0

50

12,717

16.7

685

5972

24.0

378

11,241

T12

63°

+0.0E+0

100

21,647

17.1

691

6205

24.1

382

11,286

* All CPU (central processing unit) times reported in this report are for simulation runs of approximately 30 ms on an SGI ® Origin ® 300, R14000 TM 500 megahertz (MHz).

1 degree = 0.1592 radians

As shown in figure 9, decreasing the cohesion past 6.2E-7 GPa did not appear to have any significant effects. However, at values above this, the curve was shifted to the right, indicating a delay in the initial yielding of the soil material. This agrees with conventional soil mechanics. The cohesion of the soil would delay the failure caused by chemical attraction at the molecular level.

View Alternative text
(a) Force.

View Alternative text
(b) Internal Energy

Figure 9. Cohesion parameter (Coh) variations

It is recommeded that the values of cohesion for cohesionless soil be placed at approximately 6.2D-6 GPa. This value appears to be close enough to zero, but still allows the plasticity routines to converge relatively rapidly. Larger values of cohesion rapidly digress from the zero value for cohesion. However, larger values do allow for more rapid convergence in the plasticity algorithms. It should also be examined whether the value of cohesion should be increased to compensate for the soil dilation caused by Taylor's aggregate interlock.

Angle of Internal Friction

The angle of internal friction,phi symbol, is also the slope of the shear strength envelope and, therefore, represents the effect that increasing effective normal stress has on the shear strength of the soil. For a given soil, holding all other parameters constant, an increase in the angle of internal friction should increase the shear force required to fail the soil.

For a cohesionless soil, the angle of internal friction is equal to the angle of repose--the angle at which the soil will settle into naturally. Visually, if one pours dry sand into a pile, there is a maximum angle that is achieved. This angle is the angle of repose.

A parameter study was performed varying the angle of internal friction. These values are shown in table 8. Cohesion was maintained at a constant throughout the variations of the angle of internal friction. It is critical to note that this parameter is input into LS DYNA in radians, not degrees.

As shown in figure 10, shear forces increase for decreasing angles of internal friction. This result is counter to conventional soil mechanics theories, where shear forces are known to decrease with decreasing angles of internal friction. The user was unable to determine why increasing the angle of internal friction decreased the force levels.

The baseline value of 63 degrees (1.1 radians) for an angle of internal friction was determined through physical testing performed by the user.(6) This value is recommended for NCHRP 350 soil of crushed limestone. Other types of soil that meet NCHRP 350 specifications should be tested separately.

Table 8. Internal angle of friction parameter study.

Test Identifier Angle of Internal Friction, phi symbol CPU Time (s) Peak Valley
Disp(mm) Force(kN) Internal Energy Disp(mm) Force(kN) Internal Energy
Baseline
63°
2644
18.2
718
6795
24.6
385
11,450
T1
43°
2692
18.6
768
7131
24.8
420
11,597
T2
23°
3077
18.9
769
7605
24.6
455
11,757
T3
2179
18.8
805
8529
24.5
473
13,044
1 degree = 0.1592 radians

View Alternative Text

(a) Force.

View Alternative Text

(b) Internal Energy.

Figure 10. Angle of internal friction variations.

 

DruCker-Prager Coefficient, AHYP

The Mohr-Coulomb failure criterion can be represented as a straight line in space (sigma symbolmsigma symbol-), as shown in figure 11. The point where the line cuts the sigma symbolm-axis corresponds to the tip of the hexagonal Mohr-Coulomb pyramid; it is here that the gradient of the yield surface is undefined.

 

View Alternative Text

Figure 11. Hyperbolic approximation of Mohr-Coulomb.

To avoid such angularity, Drucker-Prager introduced an inscribed cone that still possesses a vertex, but in which the "ridge" corners have been smoothed.(22) Combinations of the Mohr-Coulomb and Drucker-Prager yield surfaces can give better approximations of real failure conditions than the Drucker-Prager alone (while still avoiding the singularity of the Mohr-Coulomb yield criterion).

The developer implemented a hyperbolic approximation of the plasticity surface based on the work of Abbo and Sloan.(23) The modified yield surface is given as:

 
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.
(14)

where:

sigma symbol y = Yield surface

P = Pressure

phi symbol = Angle of internal friction

J2 = Second invariant of the stress deviator

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

Ahyp = Drucker-Prager coefficient

c = Cohesion

The elimination of the vertex singularity is also extremely useful in speeding the convergence of numerical computation, particularly where large angles of internal friction, phi symbol, and small cohesion conditions exist. This is predominantly the case with respect to NCHRP 350 crash criteria, since American Association of State Highway and Transportation Officials (AASHTO) soil specifications stipulate exactly this variety of soil.

Selecting an appropriate value for the hyperbolic coefficient, Ahyp, is important for stability in the simulation. The Drucker-Prager coefficient can be chosen as a function of the angle of internal friction and cohesion.(24) A reasonable approximation that has been found to yield good results is:

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

For values of Ahypc divided by 4cot ( phi symbol), the hyperbolic surface closely represents the Mohr-Coulomb surface. At Ahyp = 0, the original Mohr-Coulomb surface is recovered. This also restores the vertex singularity. At larger values of Ahyp, the hyperbolic surface becomes increasingly disparate from the Mohr-Coulomb surface. For numerical considerations, Ahyp should be set to values of less than c cot ( phi symbol). A graphical representation of the influence that Ahyp has on the yield surface is shown in figure 12.


View Alternative Text

Figure 12. Ahyp influence on yield surface.

For the initial simulations, Ahyp was not changed in relation to the changing cohesion and angle of internal friction parameters. It was desirable to vary each parameter separately, thus Ahyp remained constant during the cohesion and angle of internal friction parameter studies.

It was found that significant increases in forces were found when Ahyp was increased, as shown in table 9 and figure 13. The parameter Ahyp does not have significant variations from Mohr-Coulomb when it is set to very small values (1.0E-7). Keep in mind, however, that increasing Ahyp to larger values significantly deviates from the Mohr-Coulomb failure envelope.

In order to maintain similarity to the original Mohr-Coulomb failure envelope, values on the order of 1.0E-7 are recommended. For an internal angle of friction, Phimax, equal to 63 degrees (1.1 radians) and a cohesion of 6.2E-6 GPa, equation 15 yields a value of Ahyp = 1.58E-7 GPa.

Table 9. Effects of material parameter Ahyp.

Test Identifier

Ahyp (GPa)

CPU Time (s)

Peak

Valley

Disp (mm)

Force (kN)

Internal Energy

Disp (mm)

Force (kN)

Internal Energy

T16

-1.0E-7

2636

18.2

718

6795

24.6

385

11,450

T15

0.0

2636

18.6

730

7148

24.5

392

11,369

Baseline

1.0E-7

2644

18.2

718

6795

24.6

385

11,450

T14

1.0E-3

3017

17.3

839

8181

28.6

653

20,416

View Alternative Text

(a) Force.

  View Alternative Text

(b) Internal Energy.

Figure 13. Ahyp parameter variations.

Plasticity Iterations, ITERMAX

To solve the global system of nonlinear plasticity equations, an iterative approach is frequently required. Although non-iterative methods exist, such as radial return, these algorithms may lead to inaccurate results. The FHWA soil model implements an iterative plasticity scheme to solve the plasticity equations.

Iterative approaches are required because the solution to the nonlinear system is not necessarily in an equilibrium state. Strain hardening or softening may have placed the current stress state beyond the yield surface and iterative schemes, such as the Newton-Ralphson or other methods, must be used to ensure that the plasticity model converges to the true plasticity surface.

Itermax controls the number of iterations for the plasticity routine. In cases where the cohesion, Coh, is extremely small, obtaining convergence can take several iterations. This is a quality inherent in most plasticity routines and can significantly affect the accuracy and precision of a simulation.

A parameter study was performed manipulating the cohesion, Coh, and the maximum number of iterations for the plasticity routines, Itermax. The results are shown in table 10 and figure 14.

Table 10. Examination of Itermax.

Test Identifier

Cohesion, Coh

Itermax

CPU Time (s)

Peak

Valley

Disp (mm)

Force (kN)

Internal Energy

Disp (mm)

Force (kN)

Internal Energy

T10

0.0

 10

3,727

16.9

704

6351

24.0

393

11,609

T11

0.0

50

12,717

16.7

685

5972

24.0

378

11,241

T12

0.0

100

21,647

17.1

691

6205

24.1

382

11,286

T29

6.2E-6

1

814

18.5

737

7412

24.3

397

11,856

T30

6.2E-6

3

1,050

18.8

743

7449

24.5

388

11,652

T31

6.2E-6

5

1,679

18.5

732

7135

24.7

386

11,695

T32

6.2E-6

7

1,772

18.7

739

7255

25.6

403

12,677

Baseline

6.2E-6

10

2,644

18.2

718

6795

24.6

385

11,450

T27

6.2E-6

15

4,507

17.7

696

6331

24.2

376

10,984

T28

6.2E-6

100

20,770

17.4

677

5818

24.2

368

10,473

(a) Cohesion = 6.2E-06.


View Alternative Text

(b) Cohesion = 0.0

View Alternative Text

Figure 14. Itermax parameter variations.

The user was unable to determine whether there was a convergence criteria within the soil material model (increasing Itermax always generated a longer run time, with no apparent check for some type of convergence criteria being satisfied).

When cohesion is set to the recommended value of 6.2E-6, low Itermax numbers (1 through 20) give roughly the same response. Although the developer has recommended using Itermax = 10, it appears that significant CPU time can be saved with lower values of Itermax without a loss of accuracy. However, with Itermax = 100, a significant difference is seen in the results. This parameter definitely needs revisiting when the soil model is able to handle larger deformation situations (such as being able to correctly handle the direct shear test up to 100 or 200 mm of deflection, rather than only 25 mm in the current implementation).

Eccentricity Parameter, ECCEN

Eccen is the eccentricity parameter for the third stress invariant effects. To generalize the shape of the yield surface in the deviatoric plane, the developer changed the standard Mohr-Coulomb K(t) function to one developed by Klisinski.(25-26) Klisinski's yield function takes the form:

 
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.
(16)

where:

 
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.
(17)

J3 = Third invariant of the stress deviator

e = Ratio of extension strength to compression strength (Eccen)

The ratio of the extension strength to the compression strength, Eccen, is the eccentricity parameter responsible for third invariant (J3) effects. If Eccen is set to 1.0, then a circular cone surface is formed. If Eccen is set to 0.55, then a triangular surface is formed. The function K(theta) is defined for values of Eccen ranging as follows: 0.5 < eless than or equal toless than or equal to1.0.

Initial models included a value of Eccen = 0.7. This creates a relatively smooth yield surface without over-smoothing the corners of the yield surface. The authors are unaware of any physical testing or theoretical means for determining the recommended values for Eccen.

Strain Hardening Parameters, An and Et

To simulate nonlinear strain hardening behavior, the angle of internal friction, phi symbol, is increased as a function of effective plastic strain, epsiloneff plastic. It is increased as a function of Et, the amount of the nonlinear strain hardening effects desired, and An, the percentage of Phimax where nonlinear behavior begins. The increase in the angle of internal friction is given by the equation:

 
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.
(18)

For input in LS-DYNA, An is expressed as a decimal, with values between 0 and 1.0 (0 percent and 100 percent). Et affects the rate at which nonlinear hardening occurs. The authors are unaware of any physical testing or theoretical means for determining the recommended values for An or Et..

Moisture Content, MCONT

Increasing moisture content significantly reduces soil shear strength.(27) Additionally, it has been reported that marked reductions in Poisson's ratios occur because of increases in moisture content.(28) The developer addresses this by reducing the second stress invariant, J2, to produce the resulting loss in shear strength. It is critical to note that for Mcont to have any effect, the parameters Pwd1, Pwd2, and PwKsk must also be active.

Although NCHRP 350 does not give a specific moisture content criterion, for the performance of crash tests, material specification requires compaction at or near optimum moisture content. In general, optimum moisture content is around 4 percent to 5 percent, based on dry weight. Generally, the moisture content after compaction and prior to crash testing does not vary significantly. For the direct shear testing discussed in this report, moisture content was 0 percent.

As of this writing, the moisture effects are not operable within the FHWA soil model.

Pore-Water Effects on the Bulk Modulus, PWD1

To simulate the effects of moisture and air voids, the FHWA soil material model modifies the nonporous bulk modulus by using a constant relating the stiffness of the soil material before the air voids are collapsed. In equation form, this is:

 
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.
(19)

where:

Ki = Nonporous bulk modulus

ncur = Current porosity (the maximum of either 0 or (w - epsilonv))

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

epsilonv = Total volumetric strain

D1 = Parameter controlling the stiffness before the air voids are collapsed (Pwd1)

n = Porosity of the soil = (e, 1 + e)

e = Void ratio = Void Ratio is the ratio of the void volume divided by the ratio of the solid

S = Degree of saturation = Degree of Saturation is the percentage of the void volume that contains water.

rhogammasp, mc , rhow = soil density, specific gravity, moisture content, and soil density.

Appropriate values for Pwd1 must be larger than zero, but no appropriate upper limit is known. At Pwd1 = 0, the standard linear bulk modulus, Ki, is used. If Pwd1 is not set to 0.0, the bulk modulus, K, should be the fully collapsed bulk modulus. Increasing this value reduces the stiffness of the response of the soil. Information provided by the developer included the values of Pwd1, ranging from 0.0 to 10.0.

The development of excess pore pressure in a soil matrix is dependent on the portion of the pore space occupied by fluid, the rate at which the fluid can move through the soil matrix, and the driving force moving the fluid. The dissipation of excess pore pressure is a key parameter in understanding the dynamic performance of soils. In partially saturated situations, consideration of both the movement of air and fluid is necessary to define the effects of excess pore pressure on soil strength properties. In order to determine the rate of pore-pressure dissipation, the permeability of the soil matrix (both in terms of fluid and air) is a key parameter. Excess pore pressure is created by consolidation of the soil pore space, leading to localized increases in fluid/air pressure. This pressure dissipates at a time rate dependent on pressure magnitude and resistance to fluid/air flow in the soil matrix. In NCHRP 350 strong soil, the relative permeability is high, meaning that excess pore-pressure effects tend to be localized and short lived. The criteria for the decay of pore pressure relative to the soil fabric are not clear from the summary of the developer's engineering report. Without consideration of permeability, there would be no way to rationally address excess pore pressure.

Pore-Water Effects on Pore-Water Pressure, PWD2

Excess pore-water pressure reduces the total pressure and will lower the shear strength of the soil. Large pore-water pressures can cause the effective stress to disappear, causing liquefaction of the soil. To simulate the effects of excess pore-water pressure, the FHWA soil material model calculates the pore-water pressure, u, in a similar manner to that of the moisture effects on the bulk modulus:

 
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.
(20)

where:

Ksk = Nonporous bulk modulus

ncur = Current porosity (the maximum of either 0 or (w - epsilonv))

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

epsilonv = Total volumetric strain

D2 = Parameter for pore-water pressure before the air voids are collapsed (Pwd2)

n = Porosity of the soil

e = Void ratio = Void Ratio is the ratio of the void volume divided by the ratio of the solid

S = Degree of saturation = Degree of Saturation is the percentage of the void volume that contains water.

rho,gammasp, mc, rhow = soil density, specific gravity, moisture content, and water density.

Pore-water pressure is not allowed to become negative. If Pwd2 is set relatively high compared to Ksk, there is no pore-water pressure developed until the volumetric strain is greater than the strains associated with the air voids. As Pwd2 is lowered, the pore pressure starts to increase before the air voids are fully collapsed.

For an initial porosity and bulk moduli, the parameter Pwd2 can be calculated using the Skempton pore-water pressure parameter, B, as defined below:

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

This allows for the calculation of the pore-water parameter Pwd2 directly, as follows:

 
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.
(22)

Again, the comments of the previous section apply to this input. Additionally, it is assumed that excess pore pressure is used to reduce effective stress, with the commensurate influence on shear strength. In terms of shear strength, negative pore pressures, generated from the capillary rise evidenced in many soil matrixes (NCHRP 350 strong soil would not be included on this list), are important to consider in developing reasonable failure criteria. Negative pore pressure, the source of apparent cohesion in sands, can influence the peak shear strength. In practice, however, the authors are unaware of any physical testing or theoretical means for determining specific recommended values for Pwd2.

Skeleton Bulk Modulus, PWKSK

The skeleton bulk modulus is a parameter that also determines the amount of the effect that pore-water pressure has on the bulk modulus. To eliminate pore-water effects, this parameter is set to zero.

For sands, Stephen found that the dry skeleton bulk modulus was two orders of magnitude lower than the grain bulk modulus. The units of measurement (stress) for the bulk modulus are gigapascals. In practice, however, the authors are unaware of any physical testing or theoretical means for determining specific recommended values for PwKsk.

Residual Shear Strength, PHIRES

This is the angle, in radians, of the slope of the failure envelope, phi symbolult. This failure envelope defines the residual strength after the initiation of shear failure. The implementation of this value is material-dependent. In other words, there is no fixed strain at which this value is appropriate. As evidenced in the direct shear tests performed by the user, there is a gradual decrease in shear strength after the peak. The slope of this decrease is dependent on particle shape and, particularly, on density. The dilatancy and confinement of the material play important roles in this value. The residual shear strength is defined as:

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

where:

sresidual = Residual shear strength

sigma' = Effective stress

phi symbolult = Residual angle of internal friction

This strength is easily defined for most materials; however, the current limitation of the model to calculate beyond peak shear strength in the trials makes the evaluation of this parameter impossible. The rate of change from phi symbol to phi symbolult is less available, but could be determined for soils of interest and appropriate confining conditions.

Void Formation Energy, V DFM, and Volumetric Strain, DINT

When material models include strain softening, special techniques must be used to prevent mesh sensitivity. Mesh sensitivity is the tendency of a finite element model to produce significantly different results as the element size is reduced. Mesh sensitivity occurs because the softening in the model becomes concentrated in one or in a few elements.

To reduce the effects of strain softening on mesh sensitivity, the softening parameter, alpha (the strain at full damage), must be modified as the element size changes. The FHWA soil model uses an input parameter, Vdfm (Gf), that is analogous to fracture energy in 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⅓:

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

where:

xo = Initial damage threshold strain, Dint

If Gf is made increasingly small relative to K times initial damage threshold strain, Xi subscript O, times V raised to the one-third power, the softening will become progressively more brittle. Conversely, larger ratios of Gf to K times initial damage threshold strain, Xi subscript O, times V raised to the one-third power will cause the softening to become more ductile.

Dint is the volumetric strain at the peak pressure. Physically, this is the point where damage effects begin to occur, such that Dint can be conceived as the strain at the initial damage.

The authors are unaware of any physical testing or theoretical means for determining the recommended values for Vdfm or Dint.

Deletion Damage, D AMLEV, and PRINCIPAL Failure Strain, E PSMAX

As strain softening (damage) increases, the effective stiffness of the element can get very small, causing severe element distortion. One solution to this problem is deleting these distorted elements. Damlev is the percentage of damage, expressed as a decimal, that causes the deletion of an element. Epsmax is the maximum principal failure strain at which the element is deleted. It is important to note that both Damlev and Epsmax must be exceeded in order for element deletion to occur. If it is desired to turn off element deletion, Damlev should be set to zero.

In the current application, erosion of the soil elements is an unstable process and is not recommended. This is discussed further in chapter 5.

The authors are unaware of any physical testing or theoretical means for determining the recommended values for Damlev or Epsmax.

Previous | Table of Contents | Next

ResearchFHWA
FHWA
United States Department of Transportation - Federal Highway Administration