Manual for LS-DYNA 

Wood Material Model 143 

Publication No. FHWA-HRT-04-097 August 2007

Research, Development, and Technology

Turner-Fairbank Highway Research Center

6300 Georgetown Pike

McLean, VA 22101-2296
Foreword 

This report documents a wood material model that has been implemented into the dynamic 

finite element code, LS-DYNA, beginning with version 970. This material model was developed 

specifically to predict the dynamic performance of wood components used in roadside safety 

structures when undergoing a collision by a motor vehicle. This model is applicable for all 

varieties of wood when appropriate material coefficients are inserted. Default material 

coefficients for two wood varieties―southern yellow pine and Douglas fir―are stored in the 

model and can be accessed for use. 

This report is one of two that completely documents this material model. This report, Manual for 

LS-DYNA Wood Material Model 143 (FHWA-HRT-04-097), completely documents this material 

model for the user. The companion report, Evaluation of LS-DYNA Wood Material Model 143

(FHWA-HRT-04-096), completely documents the model’s performance and the accuracy of the 

results. This performance evaluation was a collaboration between the model developer and the 

model evaluator. Regarding the model performance evaluation, the developer and the evaluator 

were unable to come to a final agreement regarding the model’s performance and accuracy. 

These disagreements are itemized and thoroughly discussed in section 17 of the second report. 

This manual will be of interest to research engineers associated with the evaluation and 

crashworthy performance of roadside safety structures, particularly those engineers responsible 

for the prediction of the crash response of such structures when using the finite element code 

LS-DYNA. 

Michael F. Trentacoste 

Director, Office of Safety 

 Research and Development 

Notice 

This document is disseminated under the sponsorship of the U.S. Department of Transportation 

in the interest of information exchange. The U.S. Government assumes no liability for the use of 

the information contained in this document. This report does not constitute a standard, 

specification, or regulation. 

The U.S. Government does not endorse products or manufacturers. Trademarks or 

manufacturers' names appear in this report only because they are considered essential to the 

objective of the document. 

Quality Assurance Statement 

The Federal Highway Administration (FHWA) provides high-quality information to serve 

Government, industry, and the public in a manner that promotes public understanding. 

Standards and policies are used to ensure and maximize the quality, objectivity, utility, and 

integrity of its information. FHWA periodically reviews quality issues and adjusts its programs 

and processes to ensure continuous quality improvement. 
Technical Report Documentation Page 

1. Report No.

FHWA-HRT-04-097 

2. Government Accession No. 3. Recipient’s Catalog No. 

5. Report Date

August 2007 

4. Title and Subtitle

MANUAL FOR LS-DYNA WOOD MATERIAL MODEL 143

6. Performing Organization Code

7. Author(s)

Yvonne D. Murray

8. Performing Organization Report No.

9. Performing Organization Name and Address 10. Work Unit No. (TRAIS)

APTEK, Inc. 

1257 Lake Plaza Drive 

Colorado Springs, CO 80906-3558

11. Contract or Grant No.

DTFH61-98-C-00071 

13. Type of Report and Period Covered

Final Report 

Sept. 28, 1998–Sept. 13, 2002

12. Sponsoring Agency Name and Address

Office of Safety Research and Development 

Federal Highway Administration 

6300 Georgetown Pike 

McLean, VA 22101-2296 

14. Sponsoring Agency Code

15. Supplementary Notes

Contracting Officer’s Technical Representative (COTR): Martin Hargrave, HRDS-04

16. Abstract

An elastoplastic damage model with rate effects was developed for wood and was implemented into LS-DYNA, 

a commercially available finite element code. This manual documents the theory of the wood material model, 

describes the LS-DYNA input and output formats, and provides example problems for use as a learning tool. 

Default material property input options are provided for southern yellow pine and Douglas fir. The model was 

developed for roadside safety applications, such as wood guardrail posts impacted by vehicles; however, it 

should be applicable to most dynamic applications. 

The companion report to this manual is: 

Evaluation of LS-DYNA Wood Material Model 143 (FHWA-HRT-04-096) 

17. Key Words

Wood, LS-DYNA, orthotropic, material model, damage, rate effects, 

guardrail. 

18. Distribution Statement

No restrictions. This document is 

available to the public through the 

National Technical Information 

Service, Springfield, VA 22161.

19. Security Classif. (of this report) 

Unclassified 

20. Security Classif. (of this page)

Unclassified 

21. No. of Pages

163 

22. Price 

Form DOT F 1700.7 (8-72) Reproduction of completed page authorized
ii

SI* (MODERN METRIC) CONVERSION FACTORS 

APPROXIMATE CONVERSIONS TO SI UNITS

Symbol When You Know Multiply By To Find Symbol 

LENGTH 

in inches 25.4 millimeters mm 

ft feet 0.305 meters m 

yd yards 0.914 meters m 

mi miles 1.61 kilometers km 

AREA 

in2 square inches 645.2 square millimeters mm2

ft2 square feet 0.093 square meters m2

yd2 square yard 0.836 square meters m2

ac acres 0.405 hectares ha 

mi2 square miles 2.59 square kilometers km2

VOLUME 

fl oz fluid ounces 29.57 milliliters mL 

gal gallons 3.785 liters L 

ft3 cubic feet 0.028 cubic meters m3 

yd3 cubic yards 0.765 cubic meters m3 

NOTE: volumes greater than 1000 L shall be shown in m3

MASS 

oz ounces 28.35 grams g

lb pounds 0.454 kilograms kg

T short tons (2000 lb) 0.907 megagrams (or "metric ton") Mg (or "t") 

TEMPERATURE (exact degrees) o

F Fahrenheit 5 (F-32)/9 Celsius oC 

or (F-32)/1.8 

ILLUMINATION 

fc foot-candles 10.76 lux lx 

fl foot-Lamberts 3.426 candela/m2 cd/m2

FORCE and PRESSURE or STRESS 

lbf poundforce 4.45 newtons N 

lbf/in2 poundforce per square inch 6.89 kilopascals kPa 

APPROXIMATE CONVERSIONS FROM SI UNITS 

Symbol When You Know Multiply By To Find Symbol 

LENGTH

mm millimeters 0.039 inches in 

m meters 3.28 feet ft 

m meters 1.09 yards yd 

km kilometers 0.621 miles mi 

AREA mm2

 square millimeters 0.0016 square inches in2 

m2 square meters 10.764 square feet ft2 

m2 square meters 1.195 square yards yd2 

ha hectares 2.47 acres ac 

km2 square kilometers 0.386 square miles mi2 

VOLUME 

mL milliliters 0.034 fluid ounces fl oz 

L liters 0.264 gallons gal 

m3 cubic meters 35.314 cubic feet ft3 

m3 cubic meters 1.307 cubic yards yd3 

MASS 

g grams 0.035 ounces oz

kg kilograms 2.202 pounds lb

Mg (or "t") megagrams (or "metric ton") 1.103 short tons (2000 lb) T 

TEMPERATURE (exact degrees) o

C Celsius 1.8C+32 Fahrenheit oF 

ILLUMINATION 

lx lux 0.0929 foot-candles fc 

cd/m2 candela/m2 0.2919 foot-Lamberts fl

FORCE and PRESSURE or STRESS 

N newtons 0.225 poundforce lbf 

kPa kilopascals 0.145 poundforce per square inch lbf/in2

*SI is the symbol for th International System of Units. Appropriate rounding should be made to comply with Section 4 of ASTM E e 380. 

(Revised March 2003)
iii

Preface 

The goal of the work performed under this program, Development of DYNA3D Analysis 

Tools for Roadside Safety Applications, is to develop wood and soil material models, 

implement the models into the LS-DYNA finite element code, and evaluate the 

performance of each model through correlations with available test data.(1) 

This work was performed under Federal Highway Administration (FHWA) Contract No. 

DTFH61-98-C-00071. The FHWA Contracting Officer’s Technical Representative 

(COTR) was Martin Hargrave. 

Two reports are available for each material model. One report is a user’s manual; the 

second report is a performance evaluation. This user’s manual, Manual for LS-DYNA 

Wood Material Model 143, thoroughly documents the wood model theory; reviews the 

model input; and provides example problems for use as a learning tool. It is written by 

the developer of the model. The performance evaluation for the wood model, Evaluation 

of LS-DYNA Wood Material Model 143, documents LS-DYNA parametric studies and 

correlations with test data performed by the model developer and by a potential end 

user.(2) The reader is urged to review this user’s manual before reading the evaluation 

report. A user’s manual and evaluation report are also available for the soil model.(3,4) 

The development of the wood model was conducted by the prime contractor. The 

associated wood model evaluation effort to determine the model’s performance and the 

accuracy of the results was a collaboration between two contractors, with each 

evaluation intended to be independent of the other. The prime contractor developed and 

partially evaluated the wood model. The subcontractor performed a second independent 

evaluation of the wood model, provided finite element meshes for the evaluation 

calculations, and provided static post and bogie impact test data for correlations with the 

model. Others provided valuable material property data for clear wood pine, and static 

compression and bending test data for correlations. A final company implemented the 

wood model into the LS-DYNA finite element code. 

The developer and the evaluator were unable to come to a final agreement regarding 

several issues associated with the model’s performance and accuracy during the 

second independent evaluation of the wood model. These issues are itemized and 

thoroughly discussed in section 17 of the wood model evaluation report.(2)

v

Table of Contents 

INTRODUCTION............................................................................................................. 1 

1. THEORETICAL MANUAL ...................................................................................... 3

1.1 BEHAVIOR OF WOOD................................................................................. 3

1.2 OVERVIEW OF FORMULATION................................................................ 11

1.3 ELASTIC CONSTITUTIVE EQUATIONS.................................................... 13

1.3.1 Measured Clear Wood Moduli......................................................... 13

1.3.2 Review of Equations ....................................................................... 13

1.3.3 Default Elastic Stiffness Properties ................................................. 15

1.3.4 Orientation Vectors ......................................................................... 16

1.4 FAILURE CRITERIA................................................................................... 16

1.4.1 Measured Clear Wood Strengths.................................................... 16

1.4.2 Wood Model Failure Criteria ........................................................... 18

1.4.3 Default Strength Properties............................................................. 23

1.5 PLASTIC FLOW ......................................................................................... 23

1.5.1 Consistency Parameter Updates .................................................... 24

1.5.2 Elastoplastic Stress Updates .......................................................... 25

1.6 HARDENING.............................................................................................. 26

1.6.1 Model Overview .............................................................................. 27

1.6.2 Default Hardening Parameters........................................................ 30

1.6.3 Hardening Model Theory................................................................. 30

1.6.4 Implementation Aspects.................................................................. 32

1.7 POSTPEAK SOFTENING .......................................................................... 33

1.7.1 Degradation Model.......................................................................... 33

1.7.2 Regulating Mesh-Size Dependency................................................ 38

1.7.3 Default Damage Parameters........................................................... 44

1.7.4 Modeling Breakaway....................................................................... 47

1.8 RATE EFFECTS......................................................................................... 47

1.8.1 High Strain-Rate Data..................................................................... 48

1.8.2 Shifted Surface Model Theory......................................................... 50

1.8.3 Viscoplastic Model Theory .............................................................. 53

1.8.4 Default Rate-Effect Parameters ...................................................... 56

1.9 MODEL INPUT ........................................................................................... 56

1.10 MOISTURE EFFECTS ............................................................................... 58

1.10.1 Southern Yellow Pine...................................................................... 58

1.10.2 Douglas Fir...................................................................................... 60

1.11 TEMPERATURE EFFECTS ....................................................................... 61

1.12 VARIABILITY BY GRADE........................................................................... 65

2. USER’S MANUAL .................................................................................................... 67

2.1 LS-DYNA INPUT ........................................................................................ 67

2.2 DESCRIPTION OF PROPERTIES ............................................................. 72

2.3 FITTING THE MODEL TO THE DATA ........................................................ 75

2.4 MODEL FORMULATION............................................................................ 77

2.5 WOOD MODEL OUTPUT........................................................................... 85
vi

3. EXAMPLES MANUAL.............................................................................................. 89

3.1 SINGLE-ELEMENT SIMULATIONS ........................................................... 89

3.2 BOGIE IMPACT SIMULATION ................................................................... 92

APPENDIX A. MEASURED VARIABILITY OF SOUTHERN YELLOW PINE.............. 95

APPENDIX B. QUADRATIC EQUATIONS FIT TO MOISTURE CONTENT DATA ..... 97

APPENDIX C. ANALYTICAL FORM OF CANDIDATE FAILURE CRITERIA............ 103

APPENDIX D. GRAPHICAL COMPARISON OF CANDIDATE FAILURE 

CRITERIA ........................................................................................................... 115

APPENDIX E. DERIVATION OF CONSISTENCY PARAMETER FOR 

PLASTICITY ALGORITHM................................................................................. 137

APPENDIX F. DERIVATION OF LIMITING FUNCTION FOR HARDENING 

MODEL ............................................................................................................... 141

APPENDIX G. SINGLE-ELEMENT INPUT FILE ........................................................ 143

REFERENCES............................................................................................................ 147
vii

List of Figures 

Figure 1. Wood material properties vary with orientation. The wood material 

coordinate system does not necessarily coincide with the board coordinate 

system. Source: American Society of Civil Engineers.............................................. 3 

Figure 2. Ultimate tensile strength of Douglas fir measured in off-axis tests drops 

rapidly as the load is oriented at increasing angles to the grain. Source: Society 

of Wood Science and Technology............................................................................ 4 

Figure 3. Measured stress-strain relationships of southern yellow pine depend on 

load direction (parallel or perpendicular), load type (tensile or compressive), 

and moisture content................................................................................................ 5 

Figure 4. Temperature affects the dynamic behavior of wood posts impacted by 

bogies at 9.4 m/s...................................................................................................... 6 

Figure 5. Wood exhibits progressive softening. Source: Forest Products Laboratory. .... 7 

Figure 6. Wood exhibits modulus reduction and permanent deformation (splitting test 

data for spruce wood from Stanzi-Tschegg, et al.). Source: Kluwer Academic 

Publishers, with the permission of Springer Science and Business Media. ............. 7 

Figure 7. Variability of southern yellow pine clear wood data at 12-percent moisture 

content depends on load direction and type............................................................. 8 

Figure 8. Wood material properties vary with position. Board strength depends on 

position and size of knot. Source: Society of Wood Science and Technology. ........ 9 

Figure 9. Dynamic strength of wood increases with impact velocity in Hopkinson bar 

tests and is most pronounced in the perpendicular direction. Source: Pergamon, 

Elsevier Science Ltd............................................................................................... 10 

Figure 10. Organization of wood material model........................................................... 12 

Figure 11. Failure criteria for wood depend on four of the five invariants of a 

transversely isotropic material................................................................................ 21 

Figure 12. Failure criteria for wood produce smooth surfaces in stress space.............. 22 

Figure 13. Prepeak nonlinearity is modeled in compression with translating yield 

surfaces that allow user to specify the hardening response................................... 28 

Figure 14. Postpeak hardening is modeled in compression with positive values of 

the parameter Ghard. ............................................................................................... 29 

Figure 15. Damage d accumulates with energy τ once an initial threshold τ0 is 

exceeded................................................................................................................ 35 

Figure 16. Softening depends on the values of the damage parameters C and D 

(calculated with dmax = 1). .................................................................................... 36 

Figure 17. Softening response modeled for parallel modes of southern yellow pine..... 42 

Figure 18. Softening response modeled for perpendicular modes of southern yellow 

pine. ....................................................................................................................... 43 

Figure 19. Hopkinson bar tests indicate that the measured strength of pine 

increases with impact velocity. Source: Pergamon, Elsevier Science Ltd.(11)......... 49 

Figure 20. Hopkinson bar data indicate that strength and stiffness increase with 

strain rate. Source: EDP Sciences. ........................................................................ 50 

Figure 21. These single-element simulations demonstrate the rate-effect behavior 

of the shifted surface formulation at 500/s.............................................................. 53 

Figure 22. Two-parameter viscoplastic model is flexible in fitting data.......................... 54 
viii

Figure 23. These single-element simulations demonstrate the rate-effect behavior 

of the viscoplastic formulation at 500/s. ................................................................. 56 

Figure 24. Effect of temperature and moisture interaction on longitudinal modulus. ..... 62 

Figure 25. Temperature effects are more pronounced for the strength parallel to the 

grain than for the modulus parallel to the grain. Source: Forest Products 

Laboratory.............................................................................................................. 64 

Figure 26. Yield criteria for wood produce smooth surfaces in stress space................. 80 

Figure 27. Prepeak nonlinearity in compression is modeled with translating yield 

surfaces that allow user to specify hardening response......................................... 81 

Figure 28. Softening response modeled for parallel modes of ...................................... 83 

southern yellow pine...................................................................................................... 83 

Figure 29. Example wood model input for selection of default input parameter 

(option MAT_WOOD_PINE)................................................................................... 89 

Figure 30. Example wood model input for user specification of input parameters 

(option MAT_WOOD)............................................................................................. 90 

Figure 31. Example single-element stress-strain results for clear wood pine................ 91 

Figure 32. Deformed configuration of post at 40 ms, including fringes of damage........ 92 

Figure 33. Post deflection and cross-sectional force histories....................................... 93 

Figure 34. Measured load displacement curves of southern yellow pine exhibit 

variability in tension parallel to the grain. Source: Forest Products 

Laboratory(14).......................................................................................................... 95 

Figure 35. Measured load displacement curves of southern yellow pine exhibit 

variability in compression perpendicular to the grain. Source: Forest Products 

Laboratory.............................................................................................................. 96 

Figure 36. Effect of moisture content on tensile modulus parallel to the grain. 

Source: Forest Products Laboratory ...................................................................... 97 

Figure 37. Effect of moisture content on tensile modulus perpendicular to the grain. 

Source: Forest Products Laboratory. ..................................................................... 97 

Figure 38. Effect of moisture content on tensile strength parallel to the grain. 

Source: Forest Products Laboratory.(14) ................................................................. 98 

Figure 39. Effect of moisture content on tensile strength perpendicular to the grain. 

Source: Forest Products Laboratory.(14) ................................................................. 98 

Figure 40. Effect of moisture content on compressive strength parallel to the grain. 

Source: Forest Products Laboratory.(14) ................................................................. 99 

Figure 41. Effect of moisture content on compressive strength perpendicular to the 

grain. Source: Forest Products Laboratory.(14) ....................................................... 99 

Figure 42. Effect of moisture content on shear strength parallel to the grain. Source: 

Forest Products Laboratory.(14)............................................................................. 100 

Figure 43. Effect of moisture content on mode I fracture intensity. Source: Forest 

Products Laboratory............................................................................................. 101 

Figure 44. Effect of moisture content on mode II fracture intensity. Source: Forest 

Products Laboratory............................................................................................. 101 

Figure 45. Compressive strength variation of clear wood is readily modeled by a 

sinusoidal correction in the R-T plane. Source: Krieger Publishing 

Company.(16) ........................................................................................................ 101 
ix

Figure 46. Geometry of an off-axis test specimen. Source: Krieger Publishing 

Company.(16) ........................................................................................................ 118 

Figure 47. Most of the interactive failure criteria are in agreement with Hankinson’s 

formula for the off-axis strength of southern yellow pine in the L-T plane. ........... 120 

Figure 48. Effect of ring angle variation at 90-degree grain angle on the relative 

compression strength of four wood species. Source: Society of Wood Science 

and Technology.................................................................................................... 122 

Figure 49. Failure criteria comparison for perpendicular modes as a function of the 

ring angle. ............................................................................................................ 123 

Figure 50. Predicted effect of perpendicular confinement and extension on the 

longitudinal strength of southern yellow pine in tension and compression. .......... 125 

Figure 51. Predicted effect of parallel shear and tangential stresses on the 

longitudinal strength of southern yellow pine in tension and compression. .......... 128 

Figure 52. Predicted effect of parallel shear invariant on the longitudinal strength of 

southern yellow pine in tension and compression. ............................................... 129 

Figure 53. Predicted strength of southern yellow pine perpendicular to the grain 

(no perpendicular shear stress applied). .............................................................. 131 

Figure 54. Shape of the failure surface is sensitive to perpendicular shear strength 

if the criteria are transversely isotropic................................................................. 133 

Figure 55. Combinations of perpendicular and shear stresses that satisfy the 

failure criteria in the isotropic plane...................................................................... 135 

Figure 56. Single-element input file............................................................................. 143 

Figure 57. First continuation of single-element input file. ............................................ 144 

Figure 58. Second continuation of single-element input file. ....................................... 145 
x

List of Tables 

Table 1. Average elastic moduli of southern yellow pine............................................... 13

Table 2. Average elastic moduli of Douglas fir. ............................................................. 13

Table 3. LS-DYNA default values for the room-temperature moduli (graded or clear 

wood) of southern yellow pine and Douglas fir at saturation. ................................. 16

Table 4. Average strength data for southern yellow pine. ............................................. 17

Table 5. Average strength data for Douglas fir.............................................................. 18

Table 6. LS-DYNA default values for room-temperature clear wood strengths of 

southern yellow pine and Douglas fir at fiber saturation.* ...................................... 23

Table 7. Default hardening parameters for clear wood southern yellow pine and 

Douglas fir.............................................................................................................. 30

Table 8. LS-DYNA default values for room-temperature clear wood softening 

parameters for southern yellow pine and Douglas fir at saturation......................... 44

Table 9. Average fracture intensity data for southern yellow pine measured 

perpendicular to the grain....................................................................................... 45

Table 10. Room-temperature clear wood fracture energies for southern yellow pine 

and Douglas fir as a function of moisture content (derived from measured 

fracture intensities)................................................................................................. 46

Table 11. Strength ratios versus strain rate derived from.............................................. 50

compressive rate-effect data. ........................................................................................ 50

Table 12. Default LS-DYNA rate-effect parameters that provide the dynamic-tostatic compressive strength ratios listed in table 11 (based on units of 

milliseconds for time) for pine at 12-percent moisture content. .............................. 53

Table 13. User-supplied parameters for wood material model. ..................................... 57

Table 14. Default material property requests for wood material model. ........................ 58

Table 15. Equations fit to moisture content data for southern yellow pine..................... 59

Table 16. Equations fit to stiffness moisture content data for Douglas fir...................... 60

Table 17. Input options for modeling strength reductions by grade............................... 65
xi

Symbols 

B, D Softening parameters (parallel and perpendicular)

CI, CII Constants that relate fracture intensity to fracture energy 

Cijkl Material stiffness tensor (elastic moduli) 

Cij Material stiffness components 

c||, c⊥ Hardening-rate parameters (parallel and perpendicular) 

d, d||, d⊥, dm Scalar damage parameters (general, parallel, perpendicular, 

and max(d||,d⊥)) 

dmax||, dmax⊥ Maximum damage allowed (parallel and perpendicular) 

E11, E22, E33 Normal moduli of an orthotropic material 

EL, ET Normal moduli (wood notation) 

FM, FS Factors to scale moduli and strengths with temperature 

f||, f⊥

 Yield surface functions (parallel and perpendicular) 

f||

*

, f⊥

* Trial elastic yield surface functions (parallel and 

perpendicular) 

G||, G⊥ Hardening model translational limit functions (parallel and 

 perpendicular) 

Gf ||, Gf ⊥ Fracture energies (tension and shear) 

Gf I ||, Gf II || Parallel fracture energies (tension and shear) 

Gf I ⊥, Gf II ⊥ Perpendicular fracture energies (tension and shear) 

Ghard Hardening parameter to override perfect plasticity 

G12, G13, G23 Shear moduli of an orthotropic material 

GLT, GLR, GTR Shear moduli (wood notation) 

I1, I2, I3, I4 Stress invariants of a transversely isotropic material 

I1

*

, I2

*

, I3

*

, I4

*

 Trial elastic stress invariants 

KI, KII Fracture intensities (tension and shear) 
xii

L Element length 

MC Moisture content 

n||, n⊥ Rate-effect power parameters (parallel and perpendicular) 

N||, N⊥ Hardening initiation parameters (parallel and perpendicular) 

QT, QC Quality factors (tension/shear and compression) 

R Stress enhancement factors (ratio of dynamic to static 

strength) 

Sij Compliance coefficients (reciprocals of elastic moduli) 

S||, S⊥ Shear strengths (parallel and perpendicular) 

T Temperature 

V Impact velocity in Hopkinson pressure bar tests 

X, XT, XC Parallel wood strengths (general, tension, and compression) 

Y, YT, YC Perpendicular wood strengths (general, tension, and 

compression) 

αij Backstress tensor (and incremental backstress for hardening 

model

γ, γ||, γ⊥ Viscoplastic interpolation parameters (general, parallel, and 

perpendicular) 

εij, Δεij Strain tensor and strain increments 

ε11, ε22, ε33, Strain components of an orthotropic material 

ε12, ε13, ε23 

ε1, ε2, ε3, Strain components (shorthand notation) 

ε4, ε5, ε6 

εL, εT, εR, Strain components (wood notation) 

εLT, εLR, εTR 

L Δε& , LR Δε& , LT Δε& Strain-rate increments parallel to the grain (wood notation) 

T Δε& , R Δε& , TR Δε& Strain-rate increments perpendicular to the grain (wood 

 notation) 
xiii

|| ε& , ε ⊥

& Scalar effective strain rates (parallel and perpendicular) 

|| Δε& , Δ ⊥ ε& Scalar effective strain-rate increments (parallel and 

perpendicular) 

Δλ||, Δλ ⊥ Plasticity consistency parameters (parallel and 

perpendicular) 

Δt Time-step increment 

η General rate-effect fluidity parameter 

η||, η⊥ Tension/shear rate-effect fluidity parameters (parallel and 

perpendicular) 

ηc||, ηc⊥ Compression rate-effect fluidity parameters (parallel and 

 perpendicular) 

ρ, ρs Density of wood and of wood solid phase 

F σ 11 , F σ 22 Ultimate yield surfaces (parallel and perpendicular) 

* σ ij , σ ij ˆ , σ ij

~ , σ ij , σ ij Stress tensors (trial elastic, inviscid, inviscid with backstress, 

 viscid, and viscid with damage) 

σ11, σ22, σ33 Stress components of an orthotropic material 

σ12, σ13, σ23

σ1, σ2, σ3 Stress components (shorthand notation) 

σ4, σ5, σ6

σL, σT, σR Stress components (wood notation) 

σLT, σLR, σTR

τ||, τ⊥ Instantaneous strain energy type term for damage 

accumulation 

τ0||, τ0⊥ Initial strain energy type value for damage initiation 

νij Poisson’s ratios (indicial notation) 

νLT, νLR, νTR Major Poisson’s ratios (wood notation) 
xiv

Subscripts 

L Longitudinal or parallel 

T Transverse or perpendicular 

R Radial 

|| Parallel 

⊥ Perpendicular 

 
1

INTRODUCTION 

This manual is divided into three main sections. Section 1, “Theoretical Manual,” begins 

with a general description of wood behavior, then continues with a detailed theoretical 

description of the model as implemented in LS-DYNA, version 970. It contains 

equations for all the formulations that are implemented. These include: 

• Transversely isotropic constitutive equations. 

• Yield surfaces with plastic flow. 

• Prepeak and late-time hardening. 

• Damage-based softening with erosion. 

• Rate effects for high strain-rate applications. 

One feature of the wood material model that is different from all other LS-DYNA material 

models implemented to date are the initialization routines that provide the user with 

default input parameters for southern yellow pine and Douglas fir. These initialization 

routines set the required strengths, stiffnesses, hardening, softening, and rate-effect 

parameters as a function of moisture content, temperature, and grade. Section 1 

describes all test data and assumptions that are used to set the default input 

parameters for pine and fir. 

Section 2, “User’s Manual,” describes the wood model input and output in LS-DYNA 

format. It includes descriptions of all input parameters, methods of fitting the model to 

the data, and a brief theoretical description of the model. This section is intended to be 

a manual for users who want to apply the model without delving deeply into the theory 

of the model. Most of the information contained in this section was forwarded for 

inclusion in the LS-DYNA User’s Manual.

(5) 

Section 3, “Examples Manual,” provides example problems, including input files and 

output plots. These include four single-element simulations and one bogie impact 

simulation. These example problems are intended to help the user become familiar with 

the application of the wood material model.

3

1. THEORETICAL MANUAL 

This section documents the theory of the wood material model in detail. It begins with 

an overview of wood behavior, followed by an overview of the formulation. Then 

detailed equations are provided for each feature of the model (elasticity, plasticity, 

hardening, damage, and rate effects). Data are also tabulated for southern yellow pine 

and Douglas fir that are needed to fit the model parameters. Throughout this report, 

numerous figures, particularly those of test data, are reproduced from the various 

references cited at the end of each caption. 

1.1 BEHAVIOR OF WOOD 

Wood is a variable material; however, certain trends are evident. Stiffness and strength 

properties vary as a function of orientation between the longitudinal, tangential, and 

radial directions. Figure 1 helps to illustrate this point. The longitudinal direction is the 

fiber or grain direction. Stiffness and strength are greatest in the fiber direction. The 

tangential and radial directions are transverse to the fiber direction, and tangential and 

perpendicular to the growth rings. For modeling purposes, the distinction between the 

tangential and radial directions is not always significant. Therefore, this manual uses the 

term perpendicular to the grain when no distinction is made between the radial and 

tangential directions, and parallel to the grain to describe the longitudinal direction. 

Figure 1. Wood material properties vary with orientation. The wood material 

coordinate system does not necessarily coincide with the board coordinate 

system. Source: American Society of Civil Engineers.(6) 

Cut-Away Tree 

Showing Board 

 Location 

Enlargement 

Of Board 

With Knot 
4

Loading wood at angles to the grain has a significant effect on strength, as 

demonstrated in figure 2 for Douglas fir. The data are indicated by the dots. Hankinson’s 

formula is indicated by the surface. This formula is discussed in section 1.4. 

Figure 2. Ultimate tensile strength of Douglas fir measured in off-axis tests drops 

rapidly as the load is oriented at increasing angles to the grain. 

Source: Society of Wood Science and Technology.(7)

The failure modes and measured stress-strain relationships of wood depend on the 

direction of the load relative to the grain and the type of load (tension, compression, or 

shear). The stress-strain relationships of wood in parallel tension, perpendicular tension, 

and shear are typically linear to brittle failure, while the stress-strain relationships of 

wood in parallel compression and perpendicular compression are typically nonlinear 

and ductile. 

Another factor that affects the measured stress-strain relationships is moisture content. 

The stress-strain behavior of southern yellow pine in tension and compression is plotted 

in figure 3 as a function of moisture content. The data indicate a factor of up to three 

variations in strength with moisture content. The data also demonstrate brittle behavior 

in tension versus ductile behavior in compression. Saturation refers to the fiber 

saturation point, which is approximately 23 percent. The behavior in shear was not 

measured. 

It is important to distinguish between the modes of failure because the effect of each 

mode on the ultimate strength of the wood posts may be quite different. For example, 

impacted wood posts have been observed to fail by parallel tensile and shear 

mechanisms. Thus, ultimate wood post failure occurs in the brittle modes (parallel to the 

Ring Angle 

Relative Compression 

Grain Angle Strength 
5

grain), not the ductile modes. It is possible for perpendicular yielding to precede parallel 

failure, but not be catastrophic. 

(a) Tension parallel (b) Compression perpendicular 

Figure 3. Measured stress-strain relationships of southern yellow pine 

depend on load direction (parallel or perpendicular), load type 

(tensile or compressive), and moisture content. 

Temperature also affects the behavior of wood. This is demonstrated in figure 4 for 

wood posts impacted by bogie vehicles at 9.4 meters per second (m/s). There is a 

difference of a factor of 3 in measured response (force and velocity) between the frozen 

and unfrozen posts, which are made of southern yellow pine. The measured response 

of the Douglas fir post lies between that of the frozen and unfrozen southern yellow pine 

posts. 

Wood exhibits progressive softening, as demonstrated by the splitting test data shown 

in figure 5 for southern yellow pine. In addition, wood exhibits modulus reduction and 

permanent plastic deformation, as demonstrated by the cyclic load curve shown in 

figure 6. Modulus reduction is indicated by the decrease in the elastic loading/unloading 

slopes as strain softening progresses. Permanent deformation is evident because the 

data unloads to zero stress at nonzero values of strain. The fracture area is the area 

under the load-displacement curve following peak stress. The data were obtained from 

splitting tests conducted by Stanzl-Tschegg, et al.(8) Although the data are for spruce 

wood, similar behavior is expected for pine and fir. 
6

(a) Impact force (b) Bogie velocity reduction 

(c) Post displacement 

Figure 4. Temperature affects the dynamic behavior of wood posts 

impacted by bogies at 9.4 m/s. 
7

 (a) Softening curve (b) Testing apparatus 

Figure 5. Wood exhibits progressive softening. 

Source: Forest Products Laboratory.(9)

 (a) Softening curve (b) Testing apparatus 

Figure 6. Wood exhibits modulus reduction and permanent deformation (splitting 

test data for spruce wood from Stanzi-Tschegg, et al.). Source: Kluwer Academic 

Publishers, with the permission of Springer Science and Business Media.(8) 
8

Wood is a variable material. The variability of the clear wood properties of southern 

yellow pine is given in figure 7 for tension parallel and compression perpendicular at 

12-percent moisture content. The variability at other moisture contents is given in 

appendix A. Note that strength measurements vary by about a factor of 2 at each 

moisture content. In general, material properties vary as a function of position in a post, 

board, or test specimen. This is usually caused by natural variations in density; the 

presence of latewood and earlywood growth rings; and defects and growth 

characteristics, such as knots, checks, and shakes. Latewood is typically denser and 

stronger than earlywood. Knots, in particular, reduce the strength of wood. The 

reduction in strength depends on the knot size relative to the board size; the knot 

position (edge, center); and the wood parallel tensile strength, as shown in figure 8. 

Bogie impact tests indicate that the peak force in DS-65 posts is about 40 percent 

greater than that in grade 1 posts.(10) 

Finally, wood exhibits an increase in strength with strain rate. This is demonstrated in 

figure 9 for various wood species.(11) The stress ratio increases with impact velocity and 

is most pronounced in the perpendicular direction. The stress ratio is the dynamic-tostatic ratio measured in the Hopkinson bar tests. Rate effects are important when 

modeling vehicle collisions into wooden roadside structures. 

 (a) Tension parallel (b) Compression perpendicular 

Figure 7. Variability of southern yellow pine clear wood data at 12-percent 

moisture content depends on load direction and type. 
9

1000 pounds force per square inch (lbf/inch2

) = 6.895 megapascals (MPa), 

1 inch/inch = 1 millimeter per millimeter (mm/mm) 

Figure 8. Wood material properties vary with position. Board strength depends on 

position and size of knot. Source: Society of Wood Science and Technology.(12) 
10

(a) Parallel direction 

(b) Perpendicular direction 

Figure 9. Dynamic strength of wood increases with impact velocity in Hopkinson 

bar tests and is most pronounced in the perpendicular direction. 

Source: Pergamon, Elsevier Science Ltd.(11) 

0° Oak 

0° Redwood 

0° Pine 

0° Balsa 

6.0 

5.5 

5.0 

4.5 

4.0 

3.5 

3.0 

2.5 

2.0 

1.5 

1.0 

0.5 

0.0 

0 10 80 120 160 200 240 280 320 360

Stress ratio (Sr)

Impact velocity (m/s) 

0 10 80 120 160 200 240 280 320 360

30 

28 

26 

24 

22 

20 

18 

16 

14 

12 

10 

 8 

 6 

 4 

 2 

0

Stress ratio (Sr)

Impact velocity (m/s) 

90° Ekki 

90° Oak 

00° Redwood 

90° Pine 

90° Balsa 
11

1.2 OVERVIEW OF FORMULATION 

The wood model consists of a number of formulations that are merged together to form 

a comprehensive model: 

• Elastic constitutive equations. 

• Failure criteria. 

• Plastic flow. 

• Hardening 

• Postpeak softening. 

• Strain-rate enhancement. 

Each of these formulations is discussed separately. The flowchart in figure 10 shows 

how each formulation interacts with the others. 

Each formulation requires specific input parameters, such as stiffness, strength, and 

fracture energy. The main source of the material properties listed in this section and 

used as default input parameters is the data measured by Forests Products Laboratory 

(FPL) for southern yellow pine as a function of moisture content.(13,14) FPL data are 

plotted in appendix A. 

It is important to note that the FPL data are for clear wood (small specimens without 

defects such as knots), whereas real-world posts are graded wood (grades 3, 2, 1, or 

DS-65). Clear wood is stronger than graded wood. Nevertheless, the clear wood data 

are used as the basis for the default material properties. To account for strength 

reductions as a function of grade, reduction factors are applied to the clear wood 

strengths. Our reduction factor methodology is thoroughly discussed in section 1.12. 
12

 

Figure 10. Organization of wood material model. 

f 0 *

|| > f 0 * Elastic || ≤ Plastic 

f 0 * > ⊥ f 0 * ≤ ⊥ Elastic Plastic 

Parallel 

Perpendicular 

Exit Wood Material Model 

+1 n+1

ij nij σ α

Update Stresses With Viscoplastic Rate Effects 

22 33 23

|| 11 12 13

1 1 * 1

, ,

, ,

~ (1 )

γ γ σ σ σ

γ γ σ σ σ

σ γ σ γ σ

for

for

n

ij nij nij

⊥

+ + +

=

=

= − +

Update Stresses With Damage 

|| 22 33 23

|| 11 12 13

1 1

max( , ) , ,

, ,

(1 )

σ σ σ

σ σ σ

σ σ

d d d for

d d for

d n

ij nij

⊥

+ +

=

=

= −

0

0

0

33

22

Δ =

Δ =

Δ ⊥ =

α

α

λ

Update Elastoplastic Stresses 

ij

n ij

ijkl nij nij C α

σ

σ σ λ +

∂

∂ = − Δ ~ +1 * +1 f

0

0

11

||

Δ =

Δ =

α

λ

Δ = Δ Δt

∂

∂

∂

∂

+

∂

∂

∂

∂

− Δ =

11 || || 11 11 ||

n ||

4

4

||

n ||

1

1

||

*

||

||

c G ( - )

I

I

I f

I

f

f

α σ α ε

λ λ

λ

&

n

3

3 n

2

2

*

I

I

I f

I

f

f

⊥

⊥

⊥

⊥

⊥

⊥

∂

∂

∂

∂

+

∂

∂

∂

∂

− Δ =

λ λ

λ

t

t

Δ = Δ Δ

Δ = Δ Δ

⊥ ⊥ ⊥

⊥ ⊥ ⊥

α σ α ε

α σ α ε

&

&

c G ( - )

c G ( - )

33 33 33

22 22 22

Enter Wood Material Model 

ij

n

ij

n

ij σ α Δε

Update Trial Elastic Stresses Update Invariants 

ijkl kl

n n ij

ij C

d

ε

σ

σ + Δ − = +

1

* 1

 *

4 *3 *2 *1I I I I 

Calculate Value of Each Initial Yield Function 

* * || f f⊥
13

1.3 ELASTIC CONSTITUTIVE EQUATIONS 

1.3.1 Measured Clear Wood Moduli 

Wood guardrail posts are commonly made of southern yellow pine or Douglas fir. The 

clear wood moduli of southern yellow pine are given in table 1 in terms of the parallel

and perpendicular directions. The parallel direction refers to the longitudinal direction. 

The perpendicular direction refers to the radial or tangential direction, or any normal 

stress measurement in the R-T plane. These properties were measured as a function of 

moisture content.(13) Moisture content is more thoroughly discussed in section 1.9. 

Table 1. Average elastic moduli of southern yellow pine. 

The clear wood moduli of Douglas fir are given in table 2 as a function of the 

longitudinal (parallel), tangential, and radial directions. These properties were measured 

as a function of moisture content.(15) Only compressive moduli were measured. 

Table 2. Average elastic moduli of Douglas fir. 

Moisture Modulus of Elasticity (MPa) Poisson’s Ratio 

Content 

(%) 

Compression 

Longitudinal 

Compression 

Tangential 

Compression 

Radial LT LR TR 

 7 

13 

20 

16,345 

16,414 

15,193 

993 

779 

483 

959 

1,062 

166 

0.441 

0.449 

0.496 

0.295 

0.292 

0.274 

0.368 

0.374 

0.396 

1.3.2 Review of Equations 

Wood materials are commonly assumed to be orthotropic because they possess 

different properties in three directions―longitudinal, tangential, and radial. The elastic 

stiffness of an orthotropic material is characterized by nine independent constants. The 

nine elastic constants are E11, E22, E33, G12, G13, G23, ν12, ν13, and ν23, where E = 

Young’s modulus, G = shear modulus, and ν = Poisson’s ratio. 

The general constitutive relationship for an orthotropic material, written in terms of the 

principal material directions(16) is: 

Modulus of Elasticity (MPa) 

Poisson’s 

Moisture Ratio 

Content 

(%) 

Tension 

Parallel 

Tension 

Perpendicular 

Compression 

Parallel 

Compression 

Perpendicular LT LR 

 4 

 7 

12 

18 

Saturated 

16,469 

15,559 

15,503 

12,283 

11,283 

959 

1,000 

- 

552 

297 

18,345 

17,469 

16,572 

11,690 

7,959 

848 

766 

593 

414 

221 

0.291 

0.270 

0.260 

0.183 

0.162 

0.158 

0.133 

0.126 

0.078 

0.138 
14

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

6

5

4

3

2

1

66

55

44

13 23 33

12 22 23

11 12 13

6

5

4

3

2

1

0 0 0 0 0 2

0 0 0 0 2 0

0 0 0 2 0 0

0 0 0

0 0 0

0 0 0

ε

ε

ε

ε

ε

ε

σ

σ

σ

σ

σ

σ

C

C

C

C C C

C C C

C C C

(1)

Subscripts 1, 2, and 3 refer to the longitudinal, tangential, and radial stresses and 

strains (σ1 = σ11, σ2 = σ22, σ3 = σ33, ε1 = ε11, ε2 = ε22, and ε3 = ε33, respectively). 

Subscripts 4, 5, and 6 are a shorthand notation that refers to the shearing stresses and 

strains (σ4 = σ12, σ5 = σ23, σ6 = σ13, ε4 = ε12, ε5 = ε23, and ε6 = ε13). As an alternative 

notation for wood, it is common to substitute L (longitudinal) for 1, R (radial) for 2, and T

(tangential) for 3. The components of the constitutive matrix, Cij, are listed here in terms 

of the nine independent elastic constants of an orthotropic material: 

C11 = E11(1−ν 23ν 32 )/ Δ (2)

C22 =E22 (1−ν 31ν 13 )/ Δ (3)

C33 =E33 (1−ν 12ν 21)/ Δ (4)

C12 = (ν 21+ν 31ν 23 )E11 / Δ (5)

C13 = (ν 31+ν 21ν 32 )E11 / Δ (6)

C23 = (ν 32 +ν 12ν 31)E22 / Δ (7)

C44 =G12 (8)

C55 =G23 (9)

C66 =G13 (10)

Δ =1−ν 12ν 21 −ν 23ν 32 −ν 31ν 13 − 2ν 21ν 32ν 13 (11)
15

The following identity, relating the dependent (minor Poisson’s ratios ν21, ν31, and ν32) 

and independent elastic constants, is obtained from symmetry considerations of the 

constitutive matrix: 

= for i, j =1, 2, 3

E E j

ji

i

ν ij ν

(12) 

Another common assumption is that wood materials are transversely isotropic. This 

means that the properties in the tangential and radial directions are modeled the same 

(i.e., E22 = E33, G12 = G13, and ν12 = ν13). This reduces the number of independent 

elastic constants to five: E11, E22, ν12, G12, and G23. Furthermore, the Poisson’s ratio in 

the isotropic plane, ν23, is not an independent quantity. It is calculated from the isotropic 

relationship ν = (E – 2G)/2G, where E = E22 = E33 and G = G23. Transverse isotropy is a 

reasonable assumption if the difference between the tangential and radial properties is 

small in comparison with the difference between the tangential and longitudinal 

properties. 

The wood model formulation is transversely isotropic because the clear wood data in 

table 1 for southern yellow pine do not distinguish between the tangential and radial 

moduli. In addition, the clear wood data for Douglas fir in table 2 indicate that the 

difference between the tangential and radial moduli is less than 2 percent of the 

longitudinal modulus. 

1.3.3 Default Elastic Stiffness Properties 

Room-temperature moduli at saturation (23-percent moisture content) are listed in 

table 3. The same stiffnesses are used for graded wood as for clear wood. For southern 

yellow pine, the default Young’s moduli are average tensile values obtained from 

empirical fits to the clear wood data given in table 1 and shown in appendix B. These 

fits were published by FPL.(13) For Douglas fir, the default Young’s moduli and Poisson’s 

ratios are also obtained from empirical fits, made by the contractor, to the clear wood 

data given in table 2. The shear moduli were not measured. For both wood species, the 

parallel shear modulus is estimated from predicted elastic parameter tables for 

softwoods found in Bodig and Jayne.(16) The perpendicular shear modulus, G23, is 

calculated from the isotopic relationship between G23, E22, and ν23. For both wood 

species, the value of the perpendicular Poisson’s ratio used to estimate G23 is that 

measured for Douglas fir (because no values are listed in table 1 for pine). 
16

Table 3. LS-DYNA default values for the room-temperature moduli (graded 

or clear wood) of southern yellow pine and Douglas fir at saturation. 

 

1.3.4 Orientation Vectors 

Because the wood model is transversely isotropic, the orientation of the wood specimen 

must be set relative to the global coordinate system of LS-DYNA. The transversely 

isotropic constitutive relationships of the wood material are developed in the material 

coordinate system (i.e., the parallel and perpendicular directions). The user must define 

the orientation of the material coordinate system with respect to the global coordinate 

system. Appropriate coordinate transformations are formulated in LS-DYNA between 

the material and global coordinate systems. Such coordinate transformations are 

necessary because any differences between the grain axis and the structure axis can 

have a great effect on the structural response. 

Keep in mind that the wood grain axis may not always be perfectly aligned with the 

wood post axis because trees do not always grow straight. It is up to the analyst to set 

the alignment of the grain relative to the wood post in LS-DYNA simulations. 

1.4 FAILURE CRITERIA 

Strength variations are readily modeled with failure criteria, which are also called yield 

criteria. Failure criteria relate critical combinations of stresses or strains to failure in a 

material. Two types of failure criteria are limit and interactive criteria. With limit criteria, 

like the Maximum Stress criterion, there is no interaction between the stresses, so 

failure depends on one component of stress or strain. With interactive criteria, like the 

Hashin criterion, the stresses interact, so failure depends on more than one component 

of stress or strain. The failure stresses/strains for the interactive and maximum 

stress/strain criteria typically agree in the material principal directions (uniaxial stress 

states). The criteria disagree on what constitutes failure in off-axis directions (biaxial 

and triaxial stress states). 

1.4.1 Measured Clear Wood Strengths 

Failure criteria are formulated with coefficients that are obtained from fits to measured 

strengths (peak strength in tension and shear, yield strength in compression). The clear 

wood strengths of southern yellow pine are given in table 4 in terms of the parallel, 

perpendicular, and shear directions. The shear strength refers to the parallel-to-the-

 Southern 

Yellow Pine Douglas Fir 

E11 or EL

E22 or ET

G12 or GLT

G23 or GTR

ν12 or νLT

11,350 MPa 

247 MPa 

715 MPa 

88 MPa 

 0.16 

15,190 MPa 

 324 MPa 

784 MPa 

116 MPa 

 0.39 
17

grain direction. No shear strength was reported for the perpendicular-to-the-grain 

direction because it is difficult to measure and interpret. The modulus of rupture is 

calculated from the beam-bending test results, in which the grain runs parallel to the 

length of the beam. It is not an input parameter of the wood material model. These 

strengths were measured as a function of moisture content.(13) The saturated data are 

measured at the fiber saturation point of approximately 23-percent moisture content. 

Table 4. Average strength data for southern yellow pine. 

Moisture 

Content 

(%) 

Tension

Parallel 

(MPa) 

Tension 

Perpendicular 

(MPa) 

Compression 

Parallel 

(MPa) 

Compression 

Perpendicular 

(MPa) 

Shear 

Parallel 

(MPa) 

Modulus 

of 

Rupture 

(MPa) 

 4 

 7 

12 

18 

Saturated 

119 

136 

146 

134 

101 

3.96 

4.26 

4.50 

3.38 

1.86 

76.7 

66.8 

52.0 

33.1 

21.5 

14.8 

13.0 

10.0 

7.3 

4.0 

19.9 

19.2 

16.8 

13.5 

 8.9 

129 

121 

107 

 76 

 49 

The clear wood strengths of Douglas fir are given in table 5. These strengths were 

obtained from a variety of sources. No single source provides a complete set of 

strengths. Some sources distinguish between the radial and tangential directions, while 

others report strengths in the perpendicular direction. Whenever perpendicular 

strengths were reported, they were listed under the subheading Tangential for the 

normal strengths and LT for the shear strengths. 
18

Table 5. Average strength data for Douglas fir. 

Source 

Strength 

(MPa) 

Goodman 

and 

Bodig(7)

(12%) 

Woodward 

and 

Minor(17)

(12%) 

USDA 

Handbook(18)

(Green) 

USDA 

Handbook(18)

(12–13%) 

PattonMallory, 

et al.(19)

Tension 

Longitudinal 

Tangential 

Radial 

— 

— 

— 

123.1 

— 

3.8 

107.6 

2.3 

— 

— 

2.7 

— 

156.6 

3.2 

— 

Compression 

Longitudinal 

Tangential 

Radial 

51.9 

 5.1 

 4.3 

— 

— 

— 

23.9 

2.5 

— 

47.6 

5.3 

— 

45.2 

— 

— 

Shear 

LT 

LR 

RT

 5.4 

 7.5 

 9.0 

— 

7.7 

— 

6.6 

— 

— 

9.7 

— 

— 

8.1 

— 

— 

1.4.2 Wood Model Failure Criteria 

The strength of wood is modeled as transversely isotropic for a number of reasons. 

First, the data measured by FPL do not distinguish between the strengths in the 

tangential and radial directions. Second, the data from Goodman and Bodig suggest 

that Douglas fir is about 15 percent weaker (compressively) in the radial direction than 

in the tangential direction.(7) However, this difference in strength is small in comparison 

with the difference between the parallel and perpendicular directions. Table 4 indicates 

that the tensile strength measured parallel to the grain is about 30 to 50 times greater 

than that measured perpendicular to the grain. The compressive strength measured 

parallel to the grain is about five times greater than that measured perpendicular to the 

grain. 

The wood model failure criterion is formulated from six ultimate strength measurements 

obtained from uniaxial and pure-shear tests on wood specimens: 

XT Tensile strength parallel to the grain 

XC Compressive strength parallel to the grain 

YT Tensile strength perpendicular to the grain 

YC Compressive strength perpendicular to the grain 

S|| Shear strength parallel to the grain 

S⊥ Shear strength perpendicular to the grain 

Here, X and Y are the strengths parallel and perpendicular to the grain, respectively, 

and S is the shear strength. 
19

Seven criteria were evaluated for modeling the failure of wood. The theoretical form of 

each candidate criterion and the graphical comparisons are given in appendix C. A 

reduced form of the Modified Hashin criterion was chosen for implementation for the 

following reasons: 

• Fits off-axis and uniaxial test data well. 

• Identifies mode of failure. 

• Allows wood to fail or yield in the perpendicular modes prior to catastrophic 

failure in the parallel modes. 

• Produces a smooth surface in stress space for the parallel modes and a 

separate smooth surface for the perpendicular modes. 

• Failure strength predictions in the parallel modes are moderate in 

comparison with the extreme strengths predicted by some of the other 

criteria. 

• Failure strength predictions in the perpendicular (isotropic) plane are realistic 

under transformation of stress. 

• Provides the greatest flexibility (compared with other failure criteria) in 

modeling failure and yielding in the perpendicular modes. 

The analytical form of the Hashin criterion is different for the parallel and perpendicular 

modes. 

Parallel Modes 

For the parallel modes, the failure criterion is composed of two terms involving two of 

the five stress invariants of a transversely isotropic material. These invariants are 

1 = σ 11 I and 2

13

2 I4 = σ12 +σ . This criterion predicts that the normal and shear stresses 

are mutually weakening (i.e., the presence of shear stress reduces the strength below 

that measured in the uniaxial stress tests). This form is equivalent to that discussed in 

appendix C under Modified Hashin or Extended Yamada-Sun. Failure occurs when 

f|| ≥ 0, where: 

⎩

⎨

⎧

<

> − = + = +

0

0

1 ( )

11

11

2

||

2

13

2

12

2

2

11

|| σ

σ σ σ σ

X for

X for X

X S

f

C

T (13) 
20

Perpendicular Modes

For the perpendicular modes, the failure criterion is also composed of two terms 

involving two of the five stress invariants of a transversely isotropic material. These 

invariants are I2 = σ22 + σ33 and 22 33

2

3 23 I = σ − σ σ . This form is similar to that discussed 

in appendix C under Modified Hashin, except that two of the three reported terms are 

retained (the parallel shear stress invariant term (I4) in equation 172 is neglected). This 

is because its effect on perpendicular failure was not evaluated in appendix C and no 

test data are available to aid in the evaluation. It is desirable to keep the failure criterion 

as simple as possible unless measured data suggest otherwise. Failure occurs when f⊥

≥ 0, where: 

⎩

⎨

⎧

+ <

+ > − = − +

+ =

⊥

⊥ 0

0

1 ( ) ( )

C 22 33

T 22 33

2

22 33

2

23

2

2

22 33

σ σ

σ σ σ σ σ σ σ

Y for

Y for Y

Y S

f (14) 

Failure Surface Plots 

Four modes of failure are predicted: tension and compression failure parallel to the 

grain, and tension and compression failure perpendicular to the grain. Parallel shear 

failure is a subset of the parallel modes and perpendicular shear failure is a subset of 

the perpendicular modes. 

Each failure criterion is plotted in two dimensions in figure 11 in terms of the stress 

invariants of a transversely isotropic material. Separate plots are drawn for failure or 

yielding in the parallel and perpendicular modes. Each failure criterion is plotted in three 

dimensions in figure 12 in terms of the parallel and perpendicular stresses. Each 

criterion is a smooth surface (no corners). 
21

(a) Parallel modes 

(b) Perpendicular modes 

Figure 11. Failure criteria for wood depend on four of the five invariants 

of a transversely isotropic material. 
22

s L

s LT

s LR

(a) Parallel modes 

(b) Perpendicular modes 

Figure 12. Failure criteria for wood produce smooth surfaces in stress space. 

σL

σLT

σLR

σR σT

σTR
23

1.4.3 Default Strength Properties 

Room-temperature clear wood strengths at fiber saturation are listed in table 6. 

Strengths for southern yellow pine are average values obtained from empirical fits to the 

data previously reported by FPL in table 4 and reproduced in appendix B. Those for 

Douglas fir are based on the U.S. Department of Agriculture (USDA) Wood Handbook

strengths previously reported in table 6.(18) The shear strength perpendicular to the grain 

has been included as an input parameter even though it was not measured for southern 

yellow pine. This is because it is included in the failure criterion that was selected. Here 

it is assumed that the shear strength perpendicular to the grain is 140 percent of the 

shear strength measured parallel to the grain. This percentage was chosen because the 

perpendicular shear strength measured by Goodman and Bodig for Douglas fir is 140 

percent greater than the parallel shear strength (average of σLT and σLR).(7)

Table 6. LS-DYNA default values for room-temperature clear wood strengths 

of southern yellow pine and Douglas fir at fiber saturation.* 

 Southern 

Yellow Pine 

Douglas Fir 

XT 85.2 MPa 107.6 MPa 

XC 21.2 MPa 23.9 MPa 

YT 2.1 MPa 2.3 MPa 

YC 4.1 MPa 2.5 MPa 

S|| 9.1 MPa 6.6 MPa 

S⊥ 12.7 MPa 9.3 MPa 

 *Fiber saturation point is 23 percent for southern yellow 

 pine and 20 percent for Douglas fir. Perpendicular shear 

 strength is 140 percent of the parallel shear strength. 

The strength of graded wood posts is less than that of clear wood posts; therefore, the 

clear wood strengths in table 6 must be scaled down according to grade. One scale 

factor, QT, reduces the tensile and shear strengths as a function of grade. A second 

scale factor. QC, reduced the compressive strengths as a function of grade. Default 

scale factors for grade 1 are QT=0.43 and QC=0.63 for pine, and QT=0.40 and QC= 0.70 

for fir. Default scale factors for DS-65 are QT=0.80 and QC=0.93 for both pine and fir. 

Scale factors by grade are more thoroughly discussed in section 1.12. 

1.5 PLASTIC FLOW 

The plasticity algorithms limit the stress components once the failure criterion in 

equations 13 or 14 is satisfied. This is done by returning the stress state back to the 

yield surface.1

 

 

1 Two simple types of plasticity algorithms are those that reduce the moduli directly and those that scale 

back the stresses directly. Although simple to implement, such methods do not allow one to control plastic 

strain generation and do not necessarily satisfy the second law of thermodynamics.
24

Our traditional approach for modeling plasticity is to partition the stress and strain 

tensors into elastic and plastic parts. Partitioning is done with return mapping algorithms 

that enforce the plastic consistency conditions. Such algorithms allow one to control 

plastic strain generation. In addition, return mapping algorithms with associated flow

satisfy the second law of thermodynamics. Associated flow is discussed in appendix D. 

1.5.1 Consistency Parameter Updates 

Separate plasticity algorithms are modeled for the parallel and perpendicular modes by 

enforcing separate consistency conditions. Each consistency condition is derived in 

appendix D. The solution of each consistency condition determines the consistency 

parameters Δλ⎥⎥ and Δλ⊥. The Δλ solutions, in turn, determine the stress updates. 

Parallel Modes 

The parallel failure criterion from equation 13 is redefined as f|| (I1, I4) ≥ 0, with: 

(15)

where I1 and I4 are two of the five invariants of a transversely isotropic material. In this 

case, the expression Δλ|| is: 

n n n n

I

I

f I

I

f

f

||

4

4

||

||

1

1

||

*

||

||

λ

∂

∂

∂

λ

∂

∂

∂

λ

∂

+

∂

− Δ =

(16)

Δλ|| is calculated from specification of the total strain increments and the yield function 

( , ) *

4 * || 1 *|| f = f I I . The trial elastic stress invariants, *

1I and *

4I , are calculated from the trial 

elastic stresses (see section 1.5.2). For the failure surface function in equation 15, the 

partial derivatives in equation 16 are: 

2

1

1

|| 2

X

I

I

f = ∂

∂

(17)

2

|| 4

|| 1

I S

f = ∂

∂

(18)

2

1

11

||

1 2

X

I C I = − ∂∂

λ (19)

( , ) 1 2

||

4

2

2

1 || 1 4 = + −

S

I

X

I f I I
25

2

||

4

12

||

4 4

S

I G I = −

∂λ

∂

(20)

 

Perpendicular Modes 

The perpendicular failure criterion from equation 14 is redefined as f⊥ (I3, I4) ≥ 0, with: 

( , ) 1 2

3

2

2

2

2 3 = + −

⊥

⊥

S

I

Y

I

f I I (21)

where I2 and I3 are two of the five invariants of a transversely isotropic material. In this 

case, the expression for Δλ⊥ is: 

n n n n

I

I

I f

I

f

f

⊥

⊥

⊥

⊥

⊥

⊥

+

− Δ =

∂λ

∂

∂

∂

∂λ

∂

∂

∂

λ

3

3

2

2

*

(22)

Δλ⊥ is calculated from specification of the total strain increments and the yield function 

( , ) *

3 *2 * f⊥ = f⊥ I I , where *

2I and *

3I are the trial elastic stress invariants. For the failure 

surface function in equation 21, the partial derivatives in equation 22 are: 

2

2

2

2

Y

I

I

f = ∂∂ ⊥

(23)

2 3

1

⊥

⊥ = ∂

∂

I S

f (24)

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

+ − = +

∂

∂

⊥ ⊥

2 2 22 23 2

2 4 1 ( )

Y S

C C I I

λ (25)

2

3

23 2

2

2

22 2

2

2

22 23

3 2( ) 4

⊥ ⊥ ⊥

= + − − ∂

∂

S

I G

S

I C

Y

I C C I

λ (26)

1.5.2 Elastoplastic Stress Updates 

The stresses are readily updated from the total strain increments and the consistency 

parameters: 

n l

ijkl

n

ij

n

ij

k

f C

∂σ

∂

σ =σ − Δλ +1 * +1 ˆ (27)
26

ijkl kl nij n σ ij =σ +C Δε * +1 (28)

where: 

n = nth time step in the finite element analysis 

σ*

ij = trial elastic stress updates calculated from the total strain increments, Δεij, prior 

to application of plasticity 

Total strain increments are calculated by the finite element code from the dynamic 

equations of motion and the time step. 

Each normal stress update depends on the consistency parameters and failure surface 

functions for both the parallel (Δλ = Δλ|| and f = f||) and perpendicular (Δλ = Δλ⊥ and 

f = f⊥) modes. Each shear stress update depends on just one consistency parameter 

and failure surface function. If neither parallel nor perpendicular yielding occurs 

( 0 * f|| < and 0) * f⊥ < , then Δλ = 0 and the stress update is trivial: 1 * 1 ˆ + + = n

ij

n σ ij σ . 

For the stress updates in equation 27, the partial derivatives are: 

2

11

11

|| 2

X

f σ

∂σ

∂ = (29)

2

33

2

2

22

2

⊥

⊥ = −

Y S

f I σ

∂σ

∂

(30)

2

22

2

2

33

2

⊥

⊥ = −

Y S

f I σ

∂σ

∂

(31)

2

||

12

12

|| 2

S

f σ

∂σ

∂ = (32)

2

||

13

13

|| 2

S

f σ

∂σ

∂ = (33)

2

23

23

2

⊥

⊥ = S

f σ

∂σ

∂

(34)

1.6 HARDENING 

Wood exhibits prepeak nonlinearity in compression parallel and perpendicular to the 

grain. Perpendicular hardening was previously demonstrated in figure 3(b) for clear 
27

pine. A choice was made between two possible approaches for modeling prepeak 

nonlinearity. The first approach was to model bilinear hardening behavior by specifying 

initial and hardening moduli. This approach simulates a sudden change in modulus. 

However, real data, such as that previously shown in figure 3(b), exhibits a gradual 

change in modulus. The second approach was to model a translating yield surface that 

simulates a gradual change in modulus. This is the chosen approach. 

1.6.1 Model Overview 

Our approach is to define initial yield surfaces that harden (translate) until they coincide 

with the ultimate yield surfaces, as demonstrated in figure 13 for the parallel modes. 

Separate formulations are modeled for the parallel and perpendicular modes. The initial 

location of the yield surface determines the onset of plasticity. The rate of translation 

determines the extent of nonlinearity. 
28

(a) Initial and ultimate yield surfaces 

XC

(1-NII)XC

NIIXC

Initial Yield Strength

Ultimate Yield Strength

Increasing Rate of Translation

CII small

CII large

Strain

 | Longitudinal Stress |

(b) Stress-strain behavior 

Figure 13. Prepeak nonlinearity is modeled in compression with translating yield 

surfaces that allow user to specify the hardening response. 

Strain 

Longitudinal Stress 

Ultimate Yield Strength

Initial Yield Strength

(1-N||)Xc Increasing Rate of Translation 

C|| small 

C|| large 

N||Xc 
29

For each mode (parallel and perpendicular), two hardening parameters must be 

defined. The first parameter, N, determines the onset of nonlinearity of the loaddeflection or stress-strain curves. For example, consider the case where the user wants 

prepeak nonlinearity to initiate at 60 percent of the peak strength. The user inputs 

N = 0.4 so that 1 – N = 0.6. The second parameter, c, determines the rate of hardening

(i.e., how rapidly or gradually the load-deflection or stress-strain curves harden). If the 

user wants rapid hardening, then a large value of c is input (e.g., c = 1000). If the user 

wants gradual hardening, then a small value of c is input (e.g., c = 10). The value of c

needed for any particular fit depends on the properties of the material being modeled. It 

is selected by running simulations (single-element simulations are fastest) with various 

values of c and comparing those simulations with data. Hardening model equations are 

given in section 1.6.3. 

In addition to modeling prepeak nonlinearity as shown in figure 13, a separate 

formulation models postpeak hardening, as shown in figure 14. Instead of coinciding 

with the ultimate yield surface, the initial yield surface passes through the ultimate yield 

surface. The larger the value of the input parameter Ghard, the more pronounced the 

postpeak hardening. A zero value for Ghard will produce perfectly plastic behavior. The 

default value is zero. Currently, Ghard controls both the parallel and perpendicular 

behaviors simultaneously. 

YDM 13-MAR-02 FHWA

0.0 0.5 1.0 1.5 2.0

0

5

10

15

20

25

30

35

40

45

50

Strain (%)

Stress (MPa)

Ghard = 0.00

Ghard = 0.05

Ghard = 1.00

Figure 14. Postpeak hardening is modeled in compression 

with positive values of the parameter Ghard. 

A 

B

C

A 

B 

C 
30

1.6.2 Default Hardening Parameters 

Default hardening parameters for clear wood pine and fir are given in table 7. These 

hardening parameters were set by the contractor by correlating LS-DYNA simulations 

with compression tests conducted on 15- by 15- by 305-centimeter (cm) timbers.(2)

Values of N|| = 0.5 and N⊥ = 0.4 were fit to the timber compression data for parallel and 

perpendicular behavior, respectively. These values for N|| and N⊥ are independent of 

grade, temperature, and moisture content. Values for c|| and c⊥ depend on the wood 

grade, but are independent of temperature and moisture content. Good fits to the data 

were obtained with the following formulas: 

Perpendicular

Q

Parallel c

Q

c 2

c 2c

400 100

|| = ⊥ = (35)

Qc = compression strength-reduction factor for graded wood discussed in section 1.12. 

Table 7. Default hardening parameters for clear wood southern 

yellow pine and Douglas fir. 

Parallel Perpendicular 

Wood N|| c|| N⊥ c⊥ 

Southern Yellow Pine 

Douglas Fir 

0.5 

0.5 

400 

400 

0.4 

0.4 

100 

100 

1.6.3 Hardening Model Theory 

Parallel Modes 

The state variable that defines the translation of the yield surface is known as the 

backstress and is denoted by αij. Prepeak nonlinearity is modeled in compression, but 

not shear, so the only backstress required for the parallel modes is α11. The value of the 

backstress is α11 = 0 upon initial yield and α11 = –N|| Xc at ultimate yield (in uniaxial 

compression). The maximum backstress occurs at ultimate yield and is equal to the 

total translation of the yield surface in stress space. 

The hardening rule defines the growth of the backstress. Hardening rules are typically 

based on stress or plastic strain. Hardening is based on stress. This is accomplished by 

defining the incremental backstress: 

Δ = c G − Δ Δt || ( ) 11 || || 11 11 α σ α ε& (36)

where: 

c|| = user input parameter that determines rate of translation 

G|| = function that properly limits the increments 
31

σ11 – α11 = reduced stress that determines direction of translation (longitudinal) 

|| Δε& = effective strain-rate increment parallel to the grain 

t = time step 

These terms are internally calculated by the material model and LS-DYNA, and are 

included to keep the hardening response independent of time step, time-step scale 

factor, and strain increment. The effective strain rate is a scalar value that is equal to: 

2

LR

2

LT

2

L Δε& = Δε& + 2Δε& + 2Δε& || (37)

The function G|| restricts the motion of the yield surface so that it cannot translate 

outside the ultimate surface.(20) The functional form of G|| is determined from the 

functional form of the yield surface and the longitudinal stress definition. A brief 

derivation is given in appendix E. Thus, it is defined as: 

1 0

|| 11

11

|| = − = F N

G

σ

α

(38)

The value of the limiting function is G|| = 1 at initial yield (because α11 = 0) and G|| = 0 at 

ultimate yield (because α11 = N|| F σ11 ). Thus, G|| limits the growth of the backstress as 

the ultimate surface is approached. If postpeak hardening is active, then the minimum 

value is maintained at G|| = Ghard rather than G|| = 0. The ultimate yield surface is 

defined from equation 13 as: 

2

||

4

11 1

S

I Xc

F σ = − (39)

For the case of uniaxial compressive stress (no shear), the ultimate yield surface 

reduces to F σ11 = XC. 

Perpendicular Modes

Prepeak nonlinearity is modeled in compression, but not shear, so the backstress 

components required are α22 and α33. The value of the backstress sum is α22 + α33 = 0 

upon initial yield and α22 + α33 = –N⊥Yc at ultimate yield (biaxial compression without 

shear). The backstress increments are defined as follows: 
32

( )

( )

33 33 33

22 22 22

α σ α

α σ α

Δ = −

Δ = −

⊥ ⊥

⊥ ⊥

c G

c G

ε t

ε t

Δ Δ

Δ Δ

⊥

⊥

&

&

(40)

where: 

c⊥ = user input parameter that determines rate of translation 

Δε⊥& = effective strain-rate increment perpendicular to the grain: 

2

TR

2

R

2 Δε& = Δε&T + Δε& + 2Δε& ⊥ (41)

The functional form of G⊥ is determined from the functional form of the yield surface in 

equation 14 as: 

F N

G

22

22 33 ( ) 1

σ

α α

⊥

⊥

+ = − (42)

The value of G⊥ ranges from 1 at initial yield to 0 at ultimate yield. If postpeak hardening 

is active, then the minimum value is G⊥ = Ghard rather than G⊥ = 0. The ultimate yield 

surface is defined from equation 14 as: 

2

3

22 1

⊥

= −

S

I

Yc

F σ (43)

Consider the case of biaxial compressive stress (I3 = 0). Initially, α22 + α33 = 0, so 

G⊥ = 1. Ultimately, α22 + α33 = –N⊥ Yc and FI2 = σ22 + σ33 = Yc, so G⊥ = 0. 

1.6.4 Implementation Aspects 

The plasticity model discussed in the section 1.5 is modified to account for compressive 

hardening. Modifications are made to the failure surface definitions and the stress 

updates. 

Failure Surface Definitions With Hardening 

For the parallel modes, initial yielding occurs when f|| > 0, with: 
33

1 0 ( )

(1 )

11 2

||

2

13

2

12

2

||

2

c

2

11

|| − <

+

+ − = σ

σ σ σ

X N S

f (44)

For the perpendicular modes, initial yielding occurs when f⊥ > 0, with: 

1 0 ( )

(1 )

( )

22 33 2

22 33

2

23

2 2

c

2

22 33 − + < − + −

+ =

⊥ ⊥

⊥ σ σ

σ σ σ σ σ

Y N S

f (45)

No modifications are necessary for the tensile modes. 

Stress and Backstress Updates

Total stress is updated from the sum of the initial yield stress plus the backstress: 

ij nij n α ij =α + Δα +1 (46)

1 1 1 ˆ ~ + + + = + n

ij

n

ij

n σ ij σ α (47)

1.7 POSTPEAK SOFTENING 

In addition to predicting the critical combination of stresses at failure, modeling postfailure degradation of these stresses (softening) is particularly important. Postpeak 

degradation occurs in the tensile and shear modes of wood. This was previously 

demonstrated in figure 3. 

1.7.1 Degradation Model 

Degradation models are used to simulate postpeak softening. A choice was made 

between simple and sophisticated approaches for modeling degradation. Simple 

degradation models fit into one of three categories: instantaneous unloading, gradual 

unloading, and no unloading (constant stress after yielding, as modeled with plasticity). 

Although tensile and parallel shear failures are brittle, instantaneous unloading over one 

time step would cause dynamic instabilities. An alternative is to gradually unload over a 

number of time steps. Although simple to implement, such an ad hoc treatment will 

produce mesh-size dependency. This means that the same physical problem will 

produce different results for different mesh configurations. Two other disadvantages of 

these formulations are: (1) stiffness is not degraded in conjunction with strength and 

(2) progressive softening is independent of subsequent loading. Both of these behaviors 

are unrealistic. 
34

A more sophisticated approach is to model degradation with a damage formulation. A 

scalar damage parameter, d, transforms the stress tensor associated with the 

undamaged state, σ ij , into the stress tensor associated with the damaged state, σij : 

σ ij d σ ij = (1− ) (48)

The stress tensor σ ij is calculated by the plasticity and viscoplasticity algorithms prior 

to application of the damage model. The damage parameter ranges from zero for no 

damage to approaching unity for maximum damage. Thus, 1 – d is a reduction factor 

associated with the amount of damage. A distinct advantage of damage formulations 

over ad hoc formulations is that they degrade stiffness in conjunction with strength. 

Experimental evidence for stiffness degradation was previously shown in equation 16. 

In addition, progressive degradation does not occur over one time step; it depends on 

subsequent loading (stress spikes and transient waves will not cause spurious early 

failures because the load is not sustained). However, mesh-size dependency is still an 

issue and is discussed in section 1.7.2. This is the chosen approach. 

Damage Parameter Functional Form 

Two damage formulations are implemented for modeling degradation of wood: one 

formulation for the parallel modes and a separate model for the perpendicular modes: 

Parallel Modes 

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

−

+

+ = − − 1

1

1 ( ) ( )

max

|| || || 0

||

τ τ

τ A

d

Be

B

B

d (49)

Perpendicular Modes 

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡ −

+

+ = ⊥ ⊥

⊥

− − 1

1

1 ( ) ( )

max

|| 0 τ τ

τ C

d

De

D

D

d (50)

For each formulation, damage is specified by three user-supplied parameters. For the 

parallel modes, these parameters are A, B, and dmax||. For the perpendicular modes, 

these parameters are C, D, and dmax⊥ . The parameter dmax limits the maximum level 

of damage. It ranges between 0 and 1. No damage accumulates if dmax= 0. Typically, 

dmax≈ 1, which means that the maximum damage level attained is d ≈ 1. This, in turn, 

means that the stiffness and strength are ultimately degraded to zero. 
35

The evolution of the damage parameter d is shown in figure 15 as a function of τ. The 

strain-based energy term τ is calculated by the model. Its analytical form for both the 

parallel and perpendicular modes is discussed in subsequent paragraphs. Damage 

accumulates when τ exceeds an initial damage threshold τ0. Four curves are shown in 

figure 15 that correspond to four sets of softening parameters, C and D. The parameter 

C sets the midslope of the curve near d = 0.5 (larger values of C produce larger 

midslopes). The parameter D sets the initial slope near the threshold (smaller values of 

D produce larger initial slopes). Here, τ0 ≈ 0.055. 

Figure 15. Damage d accumulates with energy τ once 

an initial threshold τ0 is exceeded. 

Softening is demonstrated in figure 16 for tensile failure perpendicular to the grain. Four 

softening curves are shown, which correspond to the four sets of softening parameters 

previously used in figure 15. The parameters C and D shape the softening curve. Larger 

values of D produce a flatter peak. Larger values of C produce more severe softening. 
36

Figure 16. Softening depends on the values of the damage 

parameters C and D (calculated with dmax = 1). 

Damage Parameter Strain Basis 

Damage formulations are typically based on strain, stress, or energy. The wood model 

bases damage accumulation on the history of strains. One of the more famous strainbased theories is that proposed by Simo and Ju for modeling damage in isotropic 

materials such as concrete.(21) They base damage on the total strains and the 

undamaged elastic moduli. They do this by forming the undamaged elastic strain energy 

norm, . * Cijkl ij kl kl kl τ = ε ε = σ ε For simplicity,2 we have defined , *

kl Cijkl ij σ = ε where 

Cijkl is the linear elasticity tensor previously given in equation 1. One way of expanding τ

is: 

2( ) 23 * 13 23 * 12 13 * 33 12 * 22 33 * 11 22 *

11 τ = σ ε +σ ε +σ ε + σ ε +σ ε +σ ε (51)

Separate strain energy norms are implemented for modeling damage accumulation in 

the parallel and perpendicular modes. 

 

2

 The stress tensor * σ ij is not equal to the elasto-viscoplastic stress tensor ij σ nor to the damaged 

elastoplastic stress tensor ij σ . It is a fictitious stress that is based on the total strain and is defined for 

convenience. 
37

Parallel Modes 

Damage in the parallel modes is based on the following strain energy norm: 

⎪⎩

⎪

⎨

⎧

+ <

+ + ≥ =

2( ) 0

2( ) 0

13 11 * 12 13 *12

13 11 * 12 13 * 11 12 *11

||

σ ε σ ε ε

σ ε σ ε σ ε ε

τ (52)

Wood is being treated as transversely isotropic; therefore, τ|| was chosen to be the 

portion of the undamaged elastic strain energy that is associated with the parallel 

modes (specific terms from equation 51 were retained that contain the parallel normal 

and shear strains). Damage accumulates when τ|| exceeds 0||

τ . The initial threshold 0||τ

is not a user-supplied parameter. It is calculated and permanently stored by the wood 

material model at the time that the parallel failure criterion is first satisfied in tension. 

Damage accumulates in both the tensile and parallel shear modes, but not the 

compressive mode. 

Perpendicular Modes 

Damage in the perpendicular modes is based on the following strain energy norm: 

⎪

⎩

⎪

⎨

⎧

+ <

+ + + ≥

⊥ =

2 0

2 0

23 22 33 *23

23 22 33 * 33 23 * 22 33 *22

σ ε ε ε

σ ε σ ε σ ε ε ε

τ (53)

The term τ⊥ was chosen to be the portion of the undamaged elastic strain energy that is 

associated with the perpendicular modes (specific terms from equation 51 are retained 

that contain the perpendicular normal and shear strains). Damage accumulates when τ⊥

exceeds 0⊥

τ . The initial threshold 0⊥τ is not a user-supplied parameter; it is calculated 

and permanently stored by the wood material model once the perpendicular failure 

criterion is satisfied in tension. Damage accumulates in both the tensile and 

perpendicular shear modes, but not in the compressive mode. 

Strength Coupling 

Another issue is strength coupling, in which degradation in one direction affects 

degradation in another direction. If failure occurs in the parallel modes, then all six 

stress components are degraded uniformly. This is because parallel failure is 

catastrophic and will render the wood useless. The wood is not expected to carry load in 

either the parallel or perpendicular directions once the wood fibers are broken. If failure 

occurs in the perpendicular modes, then only the perpendicular stress components are 

degraded. This is because perpendicular failure is not catastrophic (the wood is 
38

expected to continue to carry the load in the parallel direction). Based on these 

assumptions, the following degradation model is implemented: 

max( ( ), ( )) dm = d τ || d τ ⊥ (54)

( ) || || d = d τ (55)

11 || 11 σ = (1− d )σ (56)

22 22 σ = (1− dm )σ (57)

33 33 σ = (1− dm )σ (58)

12 || 12 σ = (1− d )σ (59)

13 || 13 σ = (1− d )σ (60)

23 23 σ = (1− dm )σ (61)

1.7.2 Regulating Mesh-Size Dependency 

If a model is mesh-size dependent, then different mesh refinements produce different 

computational results. This is undesirable and is the result of modeling element-toelement variation in the fracture energy instead of modeling uniform fracture energy. 

Fracture energy is the area under the stress-displacement curve in the softening 

regime. If the fracture energy is not constant from element to element, then excess 

damage will accumulate in the smallest elements because the fracture energy is less in 

the smaller elements. Regulatory methods eliminate this variation and the excess 

damage accumulation. 

Fracture energy is a property of a material and special care must be taken to treat it as 

such. There are a number of approaches for regulating mesh-size dependency. One 

approach is to manually adjust the damage parameters as a function of element size to 

keep the fracture energy constant. However, this approach is not practical because the 

user would have to input different sets of damage parameters for each size element. A 

more automated approach is to include an element-length scale in the model. This is 

done by passing the element size through to the wood material model and internally 

calculating the damage parameters as a function of element size. Finally, viscous 

methods for modeling rate effects also regulate mesh-size dependency. However, if 

rate-independent calculations are performed, then viscous methods will be ineffective. 

The wood model regulates mesh-size dependency by explicitly including the element 

size in the model. The element size is calculated as the cube root of the element 

volume. Softening in the parallel and perpendicular modes is regulated separately 

because different fracture energies are measured in the parallel and perpendicular 

modes. 
39

Parallel Mode Regularization 

The relationship between the parallel fracture energy, Gf ||; the softening parameters, A

and B; and the element size, L, is: 

log(1 ) 1

|| 0||L B

AB

B Gf + ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + =τ (62)

This expression for the fracture energy was derived by integrating the analytical stressdisplacement curve: 

∫

∞ = − 0 || (1 ) T x G d X dx f (63)

where: 

x = displacement 

x0 = displacement at peak strength, XT

To accomplish the integration, the damage threshold difference is related to the 

stiffness and displacement as follows: ⎟

⎠

⎞ ⎜

⎝

⎛ − − = L

x x E 0

|| 0 11 ||

τ τ . In addition, the functional 

form of d is given by equation 49 with dmax = 1. 

To regulate mesh-size dependency, the wood model requires values for B, || I

G f , and 

|| II

G f , rather than A and B. When the parallel failure criterion is satisfied, the wood 

material model internally solves equation 62 for the value of A based on the initial 

element size, initial damage threshold, and fracture energy for the particular mode of 

parallel failure that is initiated (tensile, shear, or compressive): 

log(1 ) 1 L

||

|| 0 B

BG

B A

f

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + = τ (64)

The fracture energy varies with failure mode in the following manner: 
40

⎪

⎪

⎪

⎩

⎪

⎪

⎪

⎨

⎧

< ⎟

⎟

⎟

⎠

⎞

⎜

⎜

⎜

⎝

⎛

+

≥ ⎟

⎟

⎟

⎠

⎞

⎜

⎜

⎜

⎝

⎛ +

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

=

0

0

11 2

13

2

12

2

||

11 2

||

2

13

2

12

2

T

2

11

||

|| ||

||

σ

σ σ

σ

σ σ σ

S

G

S

G

X

G

G

II

I II

f

f f

f (65)

When failure is entirely tensile (σ11= XT, σ12 = σ13 = 0), then 

|| ||

G f = G f I and softening is 

brittle. When failure is entirely shear ( , 0) 11

2

||

2

13

2 σ 12 +σ = S σ = , then || || II

G f = G f and 

softening is more gradual. When failure is compressive (σ11 = XC, σ12 = σ13 = 0), then

Gf = ∞ and no softening occurs. 

Example stress-strain curves for these failure modes are given in figure 17. They were 

calculated for the default clear wood moduli and strengths at 12-percent moisture 

content. Default hardening and fracture energies, however, were not set at the time 

these figures were created. The parallel fracture energies simulated are seven times the 

perpendicular fracture energies. 

Perpendicular Mode Regularization 

An expression similar to equation 62 is readily obtained for the fracture energy 

perpendicular to the grain: 

log(1 ) 1

0 L D

CD

D G f ⎟ + ⎠

⎞ ⎜

⎝

⎛ + = ⊥ ⊥ τ (66)

To regulate mesh-size dependency, the wood model requires input values for D, G f I , 

and G f II , rather than C and D. 

When the perpendicular failure criterion is satisfied, the wood material model internally 

solves equation 66 for the value of C based on the initial element size, initial damage 

threshold, and fracture energy for the particular mode of perpendicular failure that is 

initiated (tensile, shear, or compressive): 

log(1 ) 1

0 L D

DG

D C

f

+ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ + =

⊥

⊥ τ (67)
41

Fracture energy varies with failure mode in the following manner: 

⎪

⎪

⎩

⎪

⎪

⎨

⎧

+ < ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

−

+ ≥ ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − +

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ +

=

⊥

⊥

⊥

⊥ ⊥

⊥

0

0 ( )

22 33

22 33

2

23

2

22 33 2

22 33

2

23

2

2

22 33

|

||

||

σ σ

σ σ σ

σ σ

σ σ σ σ σ

S G

S

G

Y

G

G

f

f

T

f

f (68)

When failure is entirely tensile (σ22 + σ33 = YT, σ23 = 0), then ⊥ ⊥ = I

G f G f and softening is 

brittle. When failure is entirely shear ( , 0) 22 33

2

22 33

2 σ 23 −σ σ = S⊥ σ +σ = , then 

⊥ ⊥ = II

G f G f and softening is more gradual. When failure is compressive 

( , 0) 22 33

2 σ 22 +σ 33 = XC σ 23 −σ σ = , then Gf = ∞ and no softening occurs. 

Stress-strain curves for these failure modes are given in figure 18. They were calculated 

for the default moduli, strengths, and fracture energies at 12-percent moisture content. 
42

(a) Tensile softening 

(b) Shear softening 

(c) Compressive yielding 

Figure 17. Softening response modeled for parallel modes 

of southern yellow pine. 
43

(a) Tensile softening 

(b) Shear softening 

(c) Compressive yielding 

Figure 18. Softening response modeled for perpendicular modes 

of southern yellow pine. 
44

1.7.3 Default Damage Parameters 

The damage model requires input of eight damage parameters―four for the parallel 

modes (B, || G f I , || G f II , and dmax||) and four for the perpendicular modes (D, ⊥ G f I , 

⊥ G f II , and dmax⊥). Default values for the maximum damage parameters are 

dmax|| = 0.9999 and dmax⊥ = 0.99. These values are slightly less than 1. This is to avoid 

potential computational difficulties associated with zero stiffness for dmax = 1. In 

addition, the parallel damage parameter is closer to 1 than the perpendicular parameter. 

This is because elements erode with maximum parallel damage, but not with maximum 

perpendicular damage (see section 1.7.4). One percent of the original perpendicular 

stiffness and strength is retained to avoid computational difficulties. Default shape 

parameters are B = D = 30. Data are not available to set the shape of the softening 

curves, so they have been arbitrarily chosen. Default fracture energies are based on the 

FPL data reported in appendix A for fracture intensities. Fracture energies are derived 

from the fracture intensities, as discussed in subsequent paragraphs. Default values for 

all eight input parameters are given in table 8 for room-temperature clear wood pine and 

fir at saturation. 

Table 8. LS-DYNA default values for room-temperature clear wood softening 

parameters for southern yellow pine and Douglas fir at saturation. 

Perpendicular Parallel 

D 

G f I ⊥

(MPa⋅cm)

G f II ⊥

(MPa⋅cm) dmax⊥ B 

G f I ||

(MPa⋅cm)

G fII||

(MPa⋅cm)

dmax||

Pine 

Fir 

30 

30 

0.0210 

0.0210 

0.0788 

0.0788 

0.99 

0.99 

30 

30

2.2344 

2.2344 

8.3843 

8.3843 

0.9999 

0.9999 

Fracture energies for Douglas fir are set equal to those for southern yellow pine. This is 

because fracture intensity data are not available for Douglas fir, and grade 1 bogie-post 

impact simulations correlated with test data suggest that this assumption is 

reasonable.(2) 

Measured Fracture Intensities 

The effect of moisture content on the mode I and mode II fracture intensities of southern 

yellow pine is given in table 9 for the perpendicular modes. These data were measured 

with the load applied in the tangential direction and the crack propagation in the 

longitudinal direction.(13) Fracture intensities were measured from compact-tension 

specimens (7.62 by 8.26 by 2.0 cm) and center-split beams (65 by 6 by 2 cm). No data 

are reported for the parallel modes. Published data indicate that the mode I fracture 
45

intensities measured parallel to the grain are about seven times those measured 

perpendicular to the grain.3(16) 

Table 9. Average fracture intensity data for southern yellow pine 

measured perpendicular to the grain. 

Moisture 

Content 

(%) 

KI

(MPa-cm1/2) 

KII

(MPa-cm1/2) 

4 

7 

12 

18 

Saturated 

4.69 

5.03 

4.66 

3.85 

2.85 

18.54 

20.41 

20.65 

18.39 

13.71 

Derived Fracture Energies 

The mode I and mode II fracture energies are related to the fracture intensities through 

the following analytical expressions: 

2

2

f II II

f I I

G C K

G C K

II

I

=

= (69)

where: 

CI and CII = constants related to the moduli of an orthotropic material:(22) 

1/ 2

11

12 66

11

11 22

1/ 2

11

12 66

11

22

1/ 2

11 22

2

2

2

2

2

2

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ +

+⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ =

⎥

⎦

⎤ ⎢

⎣

⎡

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ +

+⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ⎟

⎠

⎞ ⎜

⎝

⎛ =

S

S S

S

S S C

S

S S

S

S S S C

II

I

(70)

Here, the compliance coefficients, Sij, are the reciprocals of the elastic moduli: 

12

66

22

11

12

22

22

11

11

1 1 1

G

S

E

S

E

S

E

S = − = = = ν (71)

 

3

 Personal communication with Dr. David Kretschmann of FPL. 
46

The values of CI and CII vary with moisture content because moduli vary with moisture 

content (see appendix B). 

Tensile and shear perpendicular fracture energies as a function of moisture content are 

derived from equation 69 and FPL’s quadratic equations for fracture intensity as a 

function of moisture content. Default fracture energies at saturation are given in table 10 

for the perpendicular modes of pine and fir. These values are default values regardless 

of the grade or temperature of the wood. 

No fracture intensity or energy data are available for the parallel modes, so default 

values are based on LS-DYNA bogie-post impact simulations correlated with test 

data.(2) Correlations were made for grade 1 pine and fir posts, DS-65 pine posts, and 

frozen grade 1 pine posts. 

Good room-temperature grade 1 correlations (pine and fir) are obtained when the 

fracture energy parallel to the grain is 50 times greater than the fracture energy 

perpendicular to the grain. This parallel-to-perpendicular factor of 50 for energy is 

consistent with a parallel-to-perpendicular factor of 7 for intensity once equation 69 is 

applied (72 ≈ 50).(16) In addition, good DS-65 correlations are obtained with a factor of 

85, which indicates that the parallel fracture energy depends on the grade. To 

accommodate variation with grade, the default parallel-to-the-grain fracture energies are 

modeled as: 

f II f II T

f I f I T

G G Q

G G Q

⊥

⊥

=

=

106

106

||

|| (72)

Table 10. Room-temperature clear wood fracture energies for southern 

yellow pine and Douglas fir as a function of moisture content 

(derived from measured fracture intensities). 

Perpendicular Parallel 

Moisture 

Content 

(%) 

G f I

(MPa⋅cm) 

G f II

(MPa⋅cm)

G f I

(MPa⋅cm)

G f II

(MPa⋅cm) 

4 

7 

12 

18 

Saturated 

0.0204 

0.0227 

0.0225 

0.0219 

0.0210 

0.0768 

0.0945 

0.1028 

0.1035 

0.0788 

2.171 

2.415 

2.394 

2.330 

2.234 

8.172 

10.055 

10.938 

11.012 

8.384 

These equations indicate that the default fracture energy of clear wood, parallel to the 

grain, is 106 times greater than the default fracture energy perpendicular to the grain. 
47

Default fracture energies derived from the above equations are given in table 10 for the 

parallel modes of pine and fir. 

1.7.4 Modeling Breakaway 

Complete failure and breakup of a wood post is simulated with element erosion.4

 The 

erosion location is determined by the wood model from the physics of the problem. 

Dynamic instability is not an issue because the element erodes after it loses all strength 

and stiffness. In addition, mesh-size sensitivity is regulated through the damage 

formulation. 

Parallel damage is catastrophic because cracking occurs across the grain, which breaks 

the fibers or tubular cells of the wood. An element will automatically erode if it fails in the 

parallel mode and the parallel damage parameter exceeds d|| = 0.99. Recall that default 

dmax|| = 0.9999, so the element erodes just prior to accumulating maximum damage. All 

six stress components are degraded with parallel damage, so the element loses nearly 

all strength and stiffness before eroding. If the user sets dmax|| < 0.99, then erosion will 

not occur. 

Perpendicular damage is not catastrophic because cracking occurs between the fibers, 

causing the wood to split at relatively low strength and energy levels. The fibers are not 

broken. An element does not automatically erode if it fails in the perpendicular mode. 

This is because only three ( R, T, and RT) of the six stress components are degraded 

with perpendicular damage. Since erosion does not automatically occur, dmax⊥ is set 

to 0.99 instead of 0.9999 so that the element retains 1 percent of its elastic stiffness and 

strength. In this way, computational difficulties associated with extremely low stiffness 

and strength will hopefully be avoided. As an option, a flag is available, which, when set, 

allows elements to erode when the perpendicular damage parameter exceeds 

d⊥ = 0.989. Setting this flag is not recommended unless excessive perpendicular 

damage is causing computational difficulties. As an additional precaution, an element 

will erode if d⊥ > 0.98, and the perpendicular normal and shear strains exceed a 

predefined value of 90 percent. 

1.8 RATE EFFECTS 

A number of approaches are available to accommodate strain-rate sensitivity. Simple 

approaches scale the input strengths or moduli as a function of strain rate. More 

sophisticated approaches use viscoplastic mechanisms. Viscoplastic mechanisms 

retard the development of plasticity, which allows the stress state to exceed the static 

 

4

 Typically, erosion is more computationally efficient than the alternative approach of modeling fracture 

surfaces using tied surfaces with failure. One drawback of the tied surface approach is that an interface 

model (criteria) must be developed and validated for the tied nodes because the wood model is for 

elements, not interfaces. Another drawback is that the user must mesh the entire model with tied 

interfaces or else guess the failure location prior to running the calculation in order to specify the tied 

surface location. This is not practical. Other drawbacks are dynamic instability caused by sudden tiednode failure and mesh-size sensitivity. 
48

failure surface. Viscoplastic mechanisms also increase the stiffness in conjunction with 

strength if plastic hardening is modeled in the prepeak regime. 

Both types of formulations are easy to implement and are computationally efficient. 

Therefore, two methods are implemented. One method scales each strength as a 

function of strain rate, effectively expanding the failure surface. The second method is 

viscoplastic. Both methods are discussed here. However, the shifted surface method is 

the easiest to fit, and it is the only method that has been validated to date, so it is the 

only active option in the wood model at this time. 

1.8.1 High Strain-Rate Data 

A common test for measuring strain-rate effects is the split Hopkinson pressure bar 

(SHPB). Hopkinson bar data for pine are shown in figure 19.(11) Here, the dynamic-tostatic stress ratio is plotted as a function of impact velocity for both the parallel and 

perpendicular directions. Rate effects are more pronounced in the perpendicular 

direction than in the parallel direction. 

Also plotted in figure 19 are theoretical curves derived from shock theory (i.e., 

conservation of mass and momentum across a shock front). The theoretical stress 

enhancement factors are: 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ −

R = +

s

sV

ρ

ρ α

ρ

150 1

(0 ) 1

2 o

(73)

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ −

R = +

s s

sV

ρ

ρ α

ρ

ρ

ρ

70 1

(90 ) 1

2 o

(74)

where: 

R = ratio of dynamic-to-static stress 

ρ = average wood density 

ρs = wood solid phase density 

α = empirical shape factor for the stress-strain curve 

Values of ρs = 1500 kilograms per square meter (kg/m2) and α = 1.35 are reported by 

Reid and Peng(11) for wood and pine. For southern yellow pine with a density of ρ = 

530 kg/m2

 at 12-percent moisture content, the stress enhancement curves reduce to: 

2

1000 (0 ) 1 19 ⎟

⎠

⎞ ⎜

⎝

⎛ R = + o V (75)
49

2

1000 (90 ) 1 116 ⎟

⎠

⎞ ⎜

⎝

⎛ R = + o V (76)

where V is measured in m/s. 

These curves indicate that the stress enhancement in the perpendicular direction is 

much greater than that in the parallel direction. However, the curves tend to 

underestimate the stress enhancement, particularly at low-impact velocities (see figure 

19). To fit the rate-effect model parameters, the stress enhancement ratios need to be 

derived as a function of strain rate. To do this, impact velocity must be converted to 

strain rate. However, Reid and Peng(11) did not report this conversion. 

Figure 19. Hopkinson bar tests indicate that the measured strength of pine 

increases with impact velocity. Source: Pergamon, Elsevier Science Ltd.(11)

Bragov and Lomunov also measured the dynamic properties of pine using the SHPB 

technique.(23) Stress-strain diagrams are depicted in figure 20. Data are measured 

parallel and perpendicular to the grain. In addition, two curve sets are shown for each 

direction, which correspond to two strain rates (approximately 1000 per second (s-1) and 

500 s-1). The data indicate that both strength and parallel stiffness increase with strain 

rate. 

Impact Velocity (m/s) 

Stress ratio Sr 
50

Figure 20. Hopkinson bar data indicate that strength and stiffness increase 

with strain rate. Source: EDP Sciences.(23) 

Parallel and perpendicular strength ratios versus strain rate are needed to fit the rateeffect model. Here, those ratios are estimated. Bragov and Lomunov report a static 

parallel compressive strength of XC = 41.4 MPa for pine.(23) This strength is used to 

calculate the parallel strength ratios from figure 20. Calculated values are R(0°) = 2.3 at 

1000 per second (s) and R(0°) = 1.7 at 500/s. Bragov and Lomunov did not report a 

static perpendicular compressive strength for YC.

(23) Therefore, one cannot directly 

estimate the perpendicular strength ratio from figure 20. Instead, the Reid and Peng 

velocity-strength ratio relationships in equation 76 are used to estimate the impact 

velocities and perpendicular strength ratios.(11) Impact velocities of 262 m/s and 192 m/s 

are derived from the parallel relationship in equation 76 for R(0°) = 2.3 and R(0°) = 1.7, 

respectively. Using these velocities, perpendicular strength ratios of R(0°) = 9.0 and 

R(0°) = 5.3 are derived from the perpendicular relationship in equation 76. These ratios 

are reported in table 11. 

Table 11. Strength ratios versus strain rate derived 

from compressive rate-effect data. 

Rate (s-1) R(0°) R(90°) 

0 

500 

1000 

1.0 

1.7 

2.3 

1.0 

5.3 

9.0 

1.8.2 Shifted Surface Model Theory 

Separate rate-effect formulations are modeled for the parallel and perpendicular modes 

because the Reid and Peng data indicate that dynamic strength enhancement is more 

pronounced in the perpendicular direction than in the parallel direction.(11) The rate-
51

effect formulations increase strength with increasing strain rate by expanding each yield 

surface: 

Parallel 

XT

DYNAMIC

 = XT + EL || ε& η|| Tension (77) 

XC 

DYNAMIC

 = XC + EL || ε& ηc|| Compression (78) 

S||

DYNAMIC

 = S|| + GLT || ε& η|| Shear (79) 

Perpendicular 

YT 

DYNAMIC

 = YT + ET ε⊥& η⊥ Tension (80) 

YC 

DYNAMIC

 = YC + ET ε⊥& ηc⊥ Compression (81) 

S⊥

DYNAMIC

 = S⊥ + GTR ε⊥& η⊥ Shear (82) 

where: 

X and Y = static strengths 

X DYNAMIC and Y DYNAMIC = dynamic strengths 

Cε& η = excess stress components 

Excess stress components depend on the value of the fluidity parameter, η; the 

stiffness, C; and the effective strain rate, ε& . When rate effects are requested via the flag 

IRATE = 1, dynamic strengths are used in place of the static strengths in the yield 

surface formulations (in equations 13 and 14, or in equations 44 and 45). 

Four effective fluidity parameters are used (η||, ηc||, η⊥, and ηc⊥). Two parameters are 

needed for the parallel modes―one for the tensile and shear strengths, and a second 

for the compressive strength. Two parameters are also needed for the perpendicular 

modes. 

Each effective fluidity parameter is formulated from two input parameters: 

||

||

||

|| n

ε

η

η &= (83)
52

||

||

||

|| n

c

c ε

η

η & = (84)

⊥

⊥

⊥ = n ε

η η & (85)

⊥

⊥

⊥

⊥ = n

c

c

ε

η η & (86)

The two-parameter formulation allows the user to model a nonlinear variation in 

dynamic strength with strain rate.(24) Setting n|| = 0 or n = 0 allows the user to model a 

linear variation in dynamic strength with strain rate. The total number of input 

parameters is six. Three parameters are input for the parallel modes (η||, ηc||, and n||) 

and three parameters are input for the perpendicular modes (η||, ηc||, and n ). 

Single-element simulations performed with and without rate effects (500/s) are shown in 

figure 21 for room-temperature clear wood at 12-percent moisture content. The 

dynamic-to-static strength ratio R(0°) = 1.7, previously listed in table 11, was used to set 

the rate-effect parameters for the compressive simulation. These parameters are given 

in table 12. The plots in figure 21 indicate that different dynamic-to-static strength ratios 

naturally occur for the tension simulations than for the compressive simulations, even 

though η|| = ηc||. This is because the tensile and compressive simulations have the same 

excess stress components, but different static strengths. 
53YDM 18-JUN-99 FHWA

0.0 0.5 1.0 1.5 2.0

0

50

100

150

200

Tensile Strain Parallel to Grain (%)

Tensile Stress Parallel to Grain (MPa)

Without Rate Effects

With Rate Effects

YDM 18-JUN-99 FHWA

0.0 0.5 1.0 1.5

0

10

20

30

40

50

60

70

80

90

100

Compressive Strain Parallel to Grain (%)

Compressive Stress Parallel to Grain (MPa)

Without Rate Effects

With Rate Effects

E11 ε η||

.

(a) Tension parallel (b) Compression parallel 

Figure 21. These single-element simulations demonstrate the rate-effect 

behavior of the shifted surface formulation at 500/s. 

Table 12. Default LS-DYNA rate-effect parameters that provide the dynamic-tostatic compressive strength ratios listed in table 11 (based on units of 

milliseconds for time) for pine at 12-percent moisture content.

Parallel Perpendicular 

Wood η|| ηc|| n|| η ηc n

Pine 

Fir 

0.0045 

0.0045 

0.0045 

0.0045 

0.107 

0.107 

0.0962 

0.0962 

0.0962 

0.0962 

0.104 

0.104 

1.8.3 Viscoplastic Model Theory 

Viscoplastic algorithms allow the stress state to exceed the yield surface. The excess 

stress depends on the value of the fluidity parameter, η; the stiffness; and the strain 

rate. The viscoplastic (dynamic) stresses attained in uniaxial simulations are 

approximately those given in equations 77 through 82. 

Single- and two-parameter viscoplastic formulations are available. The single-parameter 

formulation was developed by Duvaut and Lions and models a linear variation in 

strength ratio with strain rate.(25) This is accomplished by specifying a constant value for 

each fluidity parameter (independent of strain rate). However, a linear variation is not 

necessarily appropriate for wood, as demonstrated in figure 19. A two-parameter 

formulation was recently suggested by Murray to model a nonlinear variation in strength 

ratio with strain rate.(24) The two-parameter formulation was previously given in 

equations 83 through 86. The viscoplastic formulation uses the same input parameters 

as the shifted surface formulation. 

E11 ε& η||

A

B 

B

A

A

B 

B 

A
54

The behavior of the viscoplastic model is shown in figure 22 for η = 0.0000094 s and 

four values of n. Here, the strength enhancement ratio (dynamic-to-static) has been 

plotted as a function of the strain rate. The smaller the value of n, the greater the 

strength enhancement at a given strain rate. The behavior of the two-parameter 

formulation reduces to that of the single-parameter formulation for n = 0. Thus, the 

Duvaut-Lions single-parameter model is a subset of the two-parameter model. 

Figure 22. Two-parameter viscoplastic model is flexible in fitting data. 

The two-parameter viscoplastic update scheme is implemented as follows: 

n ε

η η & = (87)

/ 1

1

Δ + = η

γ t

(88)

~ * (1 ) ij ij ij σ = − γ σ + γσ (89)

η σ ij σ ij Inviscid Solution ~ → 0 = (90)

ij ij Elastic Solution * η → ∞ σ = σ (91)

where: 
55

Δt = time step 

σ ij and σ ij

~ = viscid (dynamic) and inviscid (static) stress tensors, respectively 

* σ ij = trial elastic stress 

At each time step, σ ij

~ and * σ ij are calculated and saved by the plasticity algorithm. 

Then, equations 87 through 91 are applied to update the viscous solution. The inviscid 

and instantaneous elastic solutions are obtained as limiting cases of the 

implementation. This algorithm is easy to implement because the plasticity algorithm is 

not modified. The viscoplastic formulation is separate from and follows the plasticity 

algorithm. 

Separate formulations are implemented for the parallel (σ11, σ12, and σ13) and 

perpendicular modes (σ22, σ33, and σ23), because the Reid and Peng test data suggest 

that the rate effects depend on grain orientation.(11) The parallel stress components are 

updated with the parallel fluidity parameters, while the perpendicular stress components 

are updated with the perpendicular fluidity parameters. At this time, viscoplastic 

implementation requires only four parameters as input. These are η|| = ηc|| and n|| for the 

parallel modes, and η = ηc and n for the perpendicular modes. Therefore, one cannot 

specify separate dynamic overstresses for the tensile and compressive modes. 

Enhancement from four to six parameters could be conducted at a later date. 

Single-element simulations performed with and without viscoplastic rate effects (500/s) 

are shown in figure 23 for room-temperature clear wood pine at 12-percent moisture 

content, using the input parameters specified in table 12.5 The primary difference 

between the viscoplastic simulations and the shifted surface simulations (figure 21) is 

that the viscoplastic tensile simulation is nonlinear to peak, whereas the shifted surface 

simulation is linear. This is because plastic strains are being calculated with the 

viscoplastic formulation, not the shifted surface formulation. The user has no control 

over the linearity or nonlinearity with either formulation. 

Also note that the compressive simulations calculated with viscoplasticity attain the 

same dynamic peaks as those calculated with the shifted surface formulation. Different 

hardening parameters were used in each simulation, so the hardening is not identical. 

On the other hand, the tensile simulations calculated with the viscoplastic formulation 

attain lower peak strengths than those calculated with the shifted surface formulation. 

This is because damage accumulates in tension (not compression) and is applied to the 

stresses at the same time as viscoplasticity. Damage does not accumulate in tension 

during application of the shifted surface formulation. 

 

5

 Default hardening and softening parameters were not finalized at the time that these plots were made, 

so these figures cannot be reproduced using the current default material property values. 
56

1.8.4 Default Rate-Effect Parameters 

Default rate-effect parameters for the shifted surface formulation were previously given 

in table 12. They provide the strength enhancement ratios listed in table 11. No rateeffect data are available for Douglas fir, so the model uses the same default parameters 

as southern yellow pine. The viscoplastic algorithm uses four of the six shifted surface 

parameters (the compressive fluidity parameters are ignored). 

 

(a) Tension parallel (b) Compression parallel 

Figure 23. These single-element simulations demonstrate the rate-effect 

behavior of the viscoplastic formulation at 500/s. 

1.9 MODEL INPUT 

The previous eight subsections reviewed the behavior and properties of wood and the 

theory of the model. Here, we examine the effects of moisture content, temperature, 

and grade on model input parameters. Two methods are available for wood model 

input. One method is to provide the following input parameters: 

• Five moduli for transversely isotropic constitutive equations. 

• Six strengths for the yield criteria. 

• Four prepeak hardening parameters. 

• Eight postpeak softening parameters. 

• Six rate-effect parameters. 

These 29 parameters are listed in table 13 and were discussed in detail throughout this 

section (“Theoretical Manual”). Any consistent system of units may be used. 

A second method is to request default material properties. Default material property 

requests are listed in table 14. This method is convenient because it allows the user to 

bypass the manual input of the material parameters. 
57

Default material properties are provided for southern yellow pine and Douglas fir as a 

function of moisture content, temperature, grade, and units. A discussion of the effects 

of moisture, temperature, and grade is given in the following three subsections. Four 

consistent sets of units are available. Most of the default material properties are 

obtained from clear wood samples tested by FPL.(14) 

In addition to the choice between user-supplied and default input parameters, the user 

may also specify control parameters. Control parameters include material identification, 

density, options to plot damage, increasing the maximum number of iterations 

performed by the plasticity algorithm, turning on/off the rate effects, overriding perfect 

plasticity, and controlling perpendicular erosion. 

Table 13. User-supplied parameters for wood material model. 

Moduli Parallel Normal Modulus EL

Perpendicular Normal Modulus ET

Parallel Shear Modulus GLT

Perpendicular Shear Modulus GTR

Parallel Major Poisson’s Ratio νLT

Strengths Parallel Tensile Strength XT

Parallel Compressive Strength XC

Parallel Shear Strength S||

Perpendicular Compressive Strength YC

Perpendicular Tensile Strength YT

Perpendicular Shear Strength S⊥

Hardening Parallel Hardening Initiation N||

Parallel Hardening Rate c||

Perpendicular Hardening Initiation N⊥

Perpendicular Hardening Rate c⊥

Softening Parallel Mode I Fracture Energy Gf I || 

Parallel Mode II Fracture Energy Gf II || 

Parallel Softening B 

Parallel Maximum Damage dmax||

Perpendicular Mode I Fracture 

Energy 

Gf I ⊥

Perpendicular Mode II Fracture 

Energy 

Gf II ⊥

Perpendicular Softening D 

Perpendicular Maximum Damage dmax⊥

Rate Effects Parallel Fluidity η||

Parallel Fluidity ηC||

 Parallel Power n||

 Perpendicular Fluidity η⊥

 Perpendicular Fluidity ηC⊥

 Perpendicular Power n⊥
58

Table 14. Default material property requests for wood material model. 

Property Choices Input 

Wood Species Southern Yellow Pine Pine 

Douglas Fir Fir 

Moisture Content (%) 30% (default) 

Any moisture content (MC) 

0 

0 < MC ≤ 100 

Temperature (°C) 20 °C (default) 

Any temperature (T) 

0 

-50 ≤ T ≤ 150 

Grade Grades 1, 1D, 2, or 2D QT = 0 

DS-65 or Structural Select QT = -1 

Clear Wood QT = -2 

Quality Factors (Q) 0.0 < QT 1.0 

0.0 < QC 1.0 

Apply Quality Parameters 

Perpendicular (IQUAL) 

Yes IQUAL = 0 

No IQUAL = 1 

Units GPa, mm, ms, kg/mm3

, kN 0 

MPa, mm, ms, g/mm3

, N 1 

MPa, mm, s, Mg/mm3, N 2 

lbf/inch2, inch, s, lb-s2/inch4, lbf 3 

1.10 MOISTURE EFFECTS 

Empirical equations are implemented that specify the clear wood moduli, strengths, and 

fracture energies as a function of moisture content. The user specifies the percent 

moisture content and the model uses the appropriate moduli, strengths, and fracture 

energies. If the user does not specify the moisture content, then a moisture content of 

30 percent is used as the default. 

1.10.1 Southern Yellow Pine 

Moisture content has a significant effect on the measured moduli, strengths, and 

fracture intensities of southern yellow pine. The effect of moisture content on the elastic 

moduli was given in table 1. The effect of moisture content on strength was given in 

table 4. The effect of moisture content on the mode I and mode II fracture intensities 

was given in table 10. Plots of clear wood measurements versus moisture content are 

reproduced in appendix B. 

The empirical equations implemented for southern yellow pine are given in table 15. 

Comparisons of the equations with measured data are given in appendix B. They were 

derived by plotting the moduli, strength, and fracture intensity data as a function of 

moisture content and then fitting quadratic curves through the data.(14) Note that the 

data are highly variable. Therefore, the equations represent average clear wood 

properties. The equations for the moduli were obtained from fits to the tensile data, 
59

rather than to the compressive data. The fiber saturation point is reported as 23 percent. 

This point is the moisture content at which the cell walls are saturated with water, but no 

water exists in the cell cavities. It is generally assumed that the material properties do 

not change above this saturation point. Therefore, all material properties are held 

constant above 23 percent and set equal to those calculated by the empirical equations 

at 23 percent. The label saturated in plots indicates a moisture content of 23 percent. 

Table 15. Equations fit to moisture content data for southern yellow pine. 

P = A(MC)2 Parameter + B(MC) + C

P A B C 

Moduli 

EL Parallel Normal (MPa) 

ET Perpendicular Normal (MPa) 

νLT Parallel Poisson’s Ratio 

 -8.50 

 -2.06 

 -0.00013 

-45.3 

 17.2 

-0.00354 

16774 

 944 

0.307 

Fracture Intensities 

KIc (kN/m3/2) 

KIlc (kN/m3/2) 

-0.79 

-4.80 

10.9 

104 

447 

1505 

Strengths 

XT Tension Parallel (MPa) 

YT Tension Perpendicular (MPa) 

XC Compression Parallel (MPa) 

YC Compression Perpendicular (MPa) 

S|| Shear Parallel (MPa) 

-0.448 

-0.016 

0.011 

0.000 

-0.0226 

 10.51 

0.33 

 -3.25 

-0.555 

 0.056 

80.57 

 2.82 

90.17 

16.93 

19.86 

No data are available for parameters not listed in table 15, such as shear moduli and 

the strength in the isotropic plane. Therefore, the following assumptions are made: 

• The shear modulus parallel to the grain (G12) varies linearly with the normal 

modulus parallel to the grain (E11), as follows: 

( ) 835 619

12000

6000 619 11

12 ⎟ − ⎠

⎞ ⎜

⎝

⎛ − = +

E G (92)

This linear relationship was obtained from the predicted elastic parameters 

for softwoods found in table 3.3 of Bodig and Jayne for softwoods.(15)

• Poisson’s ratio perpendicular to the grain (ν23) is obtained from a fit to the 

Douglas fir data (see section 1.10.2). 

• The shear modulus perpendicular to the grain is obtained from the isotropic 

relationship: 
60

( ) 23

22

23 2 1+ν = E G (93)

• The shear strength perpendicular to the grain (S23) is 140 percent of the 

shear strength parallel to the grain (S12). This approximate percentage was 

obtained from measurements for four wood species reported by Goodman 

and Bodig.(7) Conversely, the USDA Wood Handbook reports that the rolling

shear strength is only 18 to 28 percent of the parallel shear strength and is 

thus quite small.(18) Nevertheless, the larger value of 140 percent is 

implemented. The effect of large versus small perpendicular shear strength 

on the shape of the yield surface is evaluated in appendix D. 

1.10.2 Douglas Fir 

There is a lack of material property data for Douglas fir. The limited data documented by 

FPL were used, and the missing information was supplemented with handbook values 

or pine data. The effect of moisture content versus elastic moduli was previously given 

in section 1. Strength measurements from various sources were given in table 5. 

Updating of the default properties is suggested as more data become available in the 

future. 

The empirical equations implemented for the Douglas fir moduli are listed in table 16. 

They were derived by fitting quadratic curves through the data from table 2. The shear 

modulus parallel to the grain (G12) varies linearly with the normal modulus parallel to the 

grain (E11), according to equation 92. 

Table 16. Equations fit to stiffness moisture content data for Douglas fir. 

P = A(MC)2 + B(MC) + C Parameter 

P A B C 

Moduli 

EL Parallel Normal (MPa) 

ET Perpendicular Normal (MPa) 

νLT Parallel Poisson’s Ratio 

νTR Perpendicular Poisson’s Ratio 

-14.3 

 -5.88 

 -0.0001154 

 0.0001649 

297.4 

108.5 

 -0.001808 

 -0.002297

14959 

 508 

 0.375 

 0.376

The equations implemented for the Douglas fir strengths are based on the equations (P) 

implemented for southern yellow pine and listed in table 16: 

 Pine

fir fir P

P Strength Strength MC MC ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ = ( )

( ) ( ) ( ) 20

20 (94)
61

Douglas fir strengths vary with moisture content in the same manner as southern yellow 

pine strengths. The term in brackets on the right of equation 94 is a scale factor with a 

value of 1.0 at 20-percent moisture content (assumed fiber saturation point). The 

strengths implemented for Douglas fir at 20-percent moisture content are the green 

material strengths listed in the USDA Wood Handbook.

(18) They were previously listed in 

table 6. The shear strength perpendicular to the grain (S23) is 140 percent of the shear 

strength parallel to the grain (S12). 

No fracture intensity data are available for Douglas fir, so the same fracture intensity 

equations and values are used as implemented for southern yellow pine. In addition, all 

Douglas fir material properties are held constant above 20-percent moisture content. 

This is because our quadratic fit to the perpendicular modulus drops to zero stiffness 

just above 22 percent. 

1.11 TEMPERATURE EFFECTS 

Typically, the moduli and strength of wood decrease as temperature increases. In 

addition, temperature interacts with moisture to influence the mechanical properties, as 

shown in figure 24.(16) This figure indicates that the temperature effects are more 

pronounced at high moisture content than at low moisture content. If the temperature 

change is not sustained, then the effect is reversible―the moduli and strengths will 

return to their original value at their original temperature. If the temperature change is 

sustained, then the effect is permanent, at least at elevated temperatures. This is 

because chemical changes occur in wood with prolonged exposure to elevated 

temperatures that degrade the wood properties and produce large reductions in 

strength. Repeated exposure to elevated temperatures also has a cumulative effect. 

Guardrail posts may be exposed to extreme temperatures for prolonged periods of time. 

Guardrail posts in warm regions that are exposed to high temperatures for long periods 

of time will have different properties than those in cold regions of the country. In 

addition, bogie tests on guardrail posts indicate that energy absorption goes down when 

posts freeze. This was previously demonstrated in figure 4. Although these observations 

on guardrail posts suggest that mechanical properties are affected by temperature, no 

clear wood data are available for southern yellow pine or Douglas fir that document the 

effect. Thus, the data from Bodig and Jayne and the USDA Wood Handbook are 

used.(16,18)
62

 

 

(a) Data reproduced from Bodig and Jayne (b) Quadratic fits to data 

 Source: Krieger Publishing Company(16)

Figure 24. Effect of temperature and moisture interaction 

on longitudinal modulus. 

Here, three methods are suggested for modeling temperature effects. One method 

takes regional variations into account through the selection of a predefined temperature 

exposure parameter tabulated by region. This method would vary the moduli, strengths, 

and fracture energies as a function of the temperature exposure parameter. However, 

this approach is not practical at this time because neither the regional data nor the clear 

wood data are available to develop such a detailed model. 

A second method specifies temperature effects by range. For example, three broad 

ranges could be modeled: 

• Low temperatures (below freezing). 

• Intermediate temperatures. 

• High temperatures. 

The user would specify a flag that indicates one of the three ranges. This method is 

simple to implement and easy to use, but does not give the user many choices. 

A third method specifies the temperature directly. With this method, the roomtemperature (20 °C) properties are scaled up or down according to the temperature 

specified. If no temperature is specified, the temperature defaults to room temperature 

(20 °C) and no scale factor is applied. This method is currently used for the default 

material properties. 

Temperature (°C) 

Relative Modulus of Elasticity 

Temperature (°C)
63

The following factor (FM) is implemented to scale the clear wood moduli as a function of 

input temperature (T): 

( ) ( 20) ( 20) 1 2 FM T = a T − + b T − + (95)

2 3

2

1 a = a (MC) + a (MC) + a (96)

2 3

2

1 b = b (MC) + b (MC) + b (97)

where: 

Coefficient a1 (°C-2) = −0.0000000377625 

Coefficient a2 (°C-2) = −0.000001416 

Coefficient a3 (°C-2) = −0.0000003125 

Coefficient b1 (°C-1) = −0.000004817 

Coefficeint b2 (°C-1) = −0.000109895 

Coefficient b3 (°C-1) = −0.000875 

All six coefficients (a1, a2, a3, b1, b2, and b3) are obtained from fits to the data previously 

shown in figure 24(a) for six wood species. Equations 95 through 97 are plotted in figure 

24(b). The data indicate that the stiffness of the wood increases when frozen and 

decreases when heated. 

The following factor (FS) is implemented to scale the clear wood strengths as a function 

of input temperature: 

FS (T ) = 2[FM (T ) −1]+1 (98)

Equation 98 is based on the data shown in figure 25. Figure 25 indicates that 

temperature has a stronger effect on strength than it has on stiffness. For temperatures 

below 20 °C, the increase in strength is twice that modeled for the increase in moduli. 

For temperatures above 20 °C, the decrease in strength is twice that modeled for the 

decrease in moduli. These plots are reproduced from the USDA Wood Handbook and 

are a composite of the results obtained from several studies.(18) 
64

(a) Strength (b) Stiffness 

Figure 25. Temperature effects are more pronounced for the strength 

parallel to the grain than for the modulus parallel to the grain. 

Source: Forest Products Laboratory.(18) 

The wood data shown in figure 25 may seem inconsistent with the bogie test data 

previously shown in figure 4. The bogie data indicate that it takes less impact force to 

break a frozen wood post than a room-temperature post. This may be because frozen 

posts are more brittle than room-temperature posts, so temperature may affect fracture 

energy. 

The author is not aware of any fracture intensity or energy data for frozen pine or fir that 

demonstrate the effect of temperature, either parallel or perpendicular to the grain. 

Therefore, the perpendicular-to-the-grain energy is modeled independent of 

temperature (using the FPL quadratic equations). However, wood is expected to 

become more brittle as temperature decreases. Therefore, parallel to the grain, a 

reduction in fracture energy is modeled upon freezing and an increase in fracture 

energy is modeled upon heating. 

The variation in fracture energy with temperature is based on correlations of dynamic 

bogie impact with frozen and room-temperature grade 1 posts. Good frozen grade 1 

correlations are obtained with a parallel-to-perpendicular energy ratio of 5, which is very 

brittle behavior. To accommodate variation with temperature, the default parallel-to-thegrain energies are modeled with a linear variation with temperature (T) between frozen 

and room-temperature values: 

RoomTemp

f I G f I T G || || = ( 0.1+ / 22.2223) (99)

RoomTemp

f I G f I T G ⊥ ⊥ = (0.1+ / 22.2223) (100)

Here, T is the temperature in degrees Celsius (°C) between 0 and 20 °C, and room 

temperature is 20 °C. For temperatures lower than 0 °C, the frozen fracture energies 
65

are the default values. For temperatures higher than 20 °C, the room-temperature 

fracture energies are the default values. 

1.12 VARIABILITY BY GRADE 

Visual grading includes an assessment of wood defects. When analyzing structures 

such as guardrail posts, the exact position of the defects will not be known, except in 

highly controlled laboratory tests. Therefore, a practical approach for addressing defects 

is to modify the material properties globally as a function of visual grade. This approach 

requires the implementation of the grade as an input parameter. 

Four input options are available for modeling strength and stiffness reductions by grade. 

The four options are listed in table 17: 

Table 17. Input options for modeling strength reductions by grade. 

User Input Grade Reduction Factors 

QT QC Wood 

QT = 0 QC Ignored 1, 1D, 2, or 2D 0.47 0.63 Pine 

0.40 0.70 Fir 

QT = −1 QC Ignored DS-65 or SEL STR 0.80 0.93 Pine & Fir 

QT = −2 QC Ignored Clear Wood 1.00 1.00 Pine & Fir 

0.0 < QT 1.0 

0.0 < QC 1.0 User-Specified QT QC Pine & Fir 

The user may request clear, high-grade (DS-65 or Select Structural), or low-grade 

(grades 1, 1D, 2, or 2D) wood. This is because significant differences were noted 

between DS-65 posts and posts of all other grades (1, 1D, 2, and 2D) tested. The 

DS-65 posts are significantly stronger than all other posts. In addition, there is no 

statistically significant difference in response among the other posts, although the dense 

posts (grades 1D and 2D) tend to absorb more energy than the low-density (grades 1 

and 2) posts. Therefore, the posts can effectively be divided into two grades. 

Our approach is to implement global strength-reduction factors as a function of grade. 

Two factors must be applied. One factor, QT, reduces the tensile and shear strengths of 

clear wood; the other factor, QC, reduces the compressive strengths. The factors are 

simultaneously applied to the parallel and perpendicular strengths. In addition to 

applying the predefined quality factors by grade, the user may directly specify the 

tensile and compressive reduction factors between 0 and 1. 

The Standard Grading Rules Handbook for southern pine lumber indicates that the 

perpendicular-to-the-grain compressive strength of lumber does not always vary with 

grade (perpendicular tensile strengths are not reported).(26) Therefore, an input flag 
66

(IFAIL) is available to allow the user to apply the strength-reduction factors parallel to 

the grain, but not perpendicular to the grain. 

Although strength-reduction factors are a global approach for modeling strength and 

stiffness reductions caused by defects, the user may want to model local defects in 

detail. This can be done by the appropriate selection of mesh, local wood properties, 

and local grain orientation. The modeling of defects in an explicit manner is expected to 

be time-consuming to set up and may be computationally intensive to run. 
67

2. USER’S MANUAL 

This section is intended to be a brief user’s manual for those users who want to run the 

model with a cursory, rather than indepth, understanding of the underlying theory and 

equations. This section includes a description of the LS-DYNA wood model input, a brief 

parameter description, and methods of fitting the parameters to data. This section 

concludes with a brief description of the wood model theory and an example output file. 

2.1 LS-DYNA INPUT 

*MAT_WOOD_{OPTION} 

This is material type 143. This is a transversely isotropic material and is available for 

solid elements in LS-DYNA. The user has the option of inputting his or her own material 

properties (<BLANK> option) or requesting default material properties for southern 

yellow pine (PINE) or Douglas fir (FIR). 

Options include: 

PINE 

FIR 

<BLANK> 

such that the keyword cards appear: 

*MAT_WOOD_PINE 

*MAT_WOOD_FIR 

*MAT_WOOD 

Define the following card for all options: 

Card Format 

Card 1 1 2 3 4 5 6 7 8 

Variable MID RO NPLOT ITERS IRATE GHARD IFAIL 

Type I F I I I F I 
68

Define the following cards for the PINE and FIR options: 

Card 2 1 2 3 4 5 6 7 8 

Variable MC TEMP QT QC UNITS IQUAL

 

Type F F F F I I 

 

Define the following cards for the <BLANK> option (do not define for PINE or FIR): 

Card 3 1 2 3 4 5 6 7 8 

Variable EL ET GLT GTR νLT

 

Type F F F F F

 

Card 4 1 2 3 4 5 6 7 8 

Variable XT XC YT YC S|| S⊥

 

Type F F F F F F

 

Card 5 1 2 3 4 5 6 7 8 

Variable Gf I ⊥ Gf II || B dmax|| Gf I ⊥ Gf II ⊥ D dmax⊥

Type F F F F F F F F

Card 6 1 2 3 4 5 6 7 8 

Variable η|| ηC|| n|| η⊥ ηC⊥ n⊥

 

Type F F F F F F

 

Card 7 1 2 3 4 5 6 7 8 

Variable N|| c|| N⊥ c⊥

 

Type F F F F

 
69

Define for all options:

Card 8 1 2 3 4 5 6 7 8 

Variable AOPT 

 

Type I 

 

Card 9 1 2 3 4 5 6 7 8 

Variable XP YP ZP A1 A2 A3 

 

Type F F F F F F 

 

Card 10 1 2 3 4 5 6 7 8 

Variable D1 D2 D3 

 

Type F F F 

 

Variable Description

 MID Material identification (a unique number has to be chosen) 

 RO Mass density 

 NPLOT Plotting options: 

EQ. 1: Maximum of parallel and perpendicular damage (default) 

EQ. 2: Perpendicular damage 

 ITERS Number of plasticity algorithm iterations (default is one iteration; 

values greater than 1 are not recommended) 

 IRATE Rate-effect options: 

 EQ. 0: Rate-effect model turned off (default) 

 EQ. 1: Rate-effect model turned on 

 GHARD Perfect plasticity override (values greater than or equal to zero are 

allowed). Positive values model late-time hardening in compression (an 

increase in strength with increasing strain). A zero value models perfect 

plasticity (no increase in strength with increasing strain). The default is 

zero. 
70

 IFAIL Erosion perpendicular to the grain: 

 EQ. 0: No (default) 

 EQ. 1: Yes (not recommended except for debugging) 

Define for PINE and FIR options: 

Variable Description 

 MC Percent moisture content (if left blank, moisture content defaults to 

saturated at 30 percent) 

 TEMP Temperature in °C (if left blank, temperature defaults to room temperature 

at 20 °C) 

QT Quality factor options (these quality factors reduce the clear wood 

tension/shear and compression strengths as a function of grade): 

EQ. 0: Grades 1, 1D, 2, 2D 

Predefined strength-reduction factors are: 

Pine: QT = 0.47 in tension/shear 

QC = 0.63 in compression 

Fir: QT = 0.40 in tension/shear 

QC = 0.70 in compression 

EQ. -1:DS-65 or SEL STR 

Predefined strength-reduction factors are: 

QT = 0.80 in tension/shear 

QC = 0.93 in compression 

 EQ. -2:Clear wood 

QT = 1.0 in tension/shear 

QC = 1.0 in compression 

GT. 0: User-defined quality factor in tension (values greater than 0 and 

less than or equal to 1 are expected; values greater than 1 are 

allowed, but may not be realistic. 

QC User-defined quality factor in compression (This input value is used if 

QT > 0. Values greater than 0 and less than or equal to 1 are expected. 

Values greater than 1 are allowed, but may not be realistic. If left blank 

when QT > 0, a default value of QC = QT is used.) 
71

 UNITS Unit options: 

EQ. 0: gigapascals (GPa), mm, milliseconds (ms), kilograms per cubic 

 millimeter (kg/mm3), kilonewtons (kN) 

EQ. 1: MPa, ms, grams per cubic millimeter (g/mm3

), newtons (N) 

EQ. 2: MPa, mm, s, megagrams per cubic millimeter (Mg/mm3

), N 

EQ. 3: lbf/inch2

, inch, s, pound second squared per inch to the fourth 

power (lb-s2/inch4), pounds force (lbf) 

 IQUAL Apply quality factors perpendicular to the grain: 

 EQ. 0: Yes (default) 

 EQ. 1: No 

Remarks: Material property data are for clear wood (small samples without defects 

such as knots), whereas real structures are composed of graded wood. Clear wood is 

stronger than graded wood. Quality factors (strength-reduction factors) are applied to 

the clear wood strengths to account for reductions in strength as a function of grade. 

One quality factor (QT) is applied to the tensile and shear strengths. A second quality 

factor (QC) is applied to the compressive strengths. As an option, predefined quality 

factors are provided based on correlations between LS-DYNA calculations and test data 

for pine and fir posts impacted by bogie vehicles. By default, quality factors are applied 

to the parallel strengths and to the perpendicular strengths. An option is available 

(IQUAL) to eliminate application perpendicular to the grain. 

Define for <BLANK> option only: 

Variable Description 

EL Parallel normal modulus 

ET Perpendicular normal modulus 

GLT Parallel shear modulus (GLR = GLT) 

GTR Perpendicular shear modulus 

νLT Parallel major Poisson’s ratio 

XT Parallel tensile strength 

XC Parallel compressive strength 

YT Perpendicular tensile strength 

YC Perpendicular compressive strength 

S|| Parallel shear strength 

S⊥ Perpendicular shear strength 

 Gf I || Parallel fracture energy in tension 

 Gf II || Parallel fracture energy in shear 

B Parallel softening parameter 

dmax|| Parallel maximum damage 

 Gf I ⊥ Perpendicular fracture energy in tension 

 Gf II ⊥ Perpendicular fracture energy in shear 

D Perpendicular softening parameter 

dmax⊥ Perpendicular maximum damage 
72

η|| Parallel fluidity parameter in tension/shear 

ηc|| Parallel fluidity parameter in compression 

n|| Parallel power 

η⊥ Perpendicular fluidity parameter in tension/shear 

ηc⊥ Perpendicular fluidity parameter in compression 

n⊥ Perpendicular power 

N|| Parallel hardening initiation 

 c|| Parallel hardening rate 

 N⊥ Perpendicular hardening initiation 

c⊥ Perpendicular hardening rate 

Define for all options: 

AOPT material axes option (see MAT_OPTIONTROPIC_ELASTIC for a more complete 

description): 

EQ. 0: Locally orthotropic with material axes determined by element 

EQ. 1: Locally orthotropic with material axes determined by a point in 

 space and the global location of the element center; this is the 

a-direction 

EQ. 2: Globally orthotropic with material axes determined by vectors 

defined below, as with *DEFINE_COORDINATE_VECTOR 

XP, YP, ZP Coordinates of point p for AOPT = 1 

A1, A2, A3 Coordinates of vector a for AOPT = 2 

D1, D2, D3 Components of vector d for AOPT = 2 

Remarks: One common option is AOPT = 2. The user defines vectors a and d. 

Typically, a is the parallel-to-the-grain direction and d is one of the perpendicular-to-thegrain directions. Then, a x d = c and c x a = b, where a, b, and c are the principal 

material axes. 

2.2 DESCRIPTION OF PROPERTIES 

Wood is generally considered an orthotropic material with different properties in the 

longitudinal, tangential, and radial directions. This is a transversely isotropic model 

(simplification of orthotropic) in which the properties in the tangential and radial 

directions are modeled the same. For simplicity, the longitudinal direction is referred to 

as the parallel-to-the-grain direction and the tangential and/or radial directions as the 

perpendicular-to-the-grain direction. This is an elastoplastic damage model with rate 

effects. Separate elliptical yield surfaces and scalar damage formulations are modeled 

for the parallel and perpendicular modes. The damage formulations admit progressive 

degradation of the tensile and shear stress components while retaining perfect plasticity 

in compression. Hardening formulations are available, based on a translating yield 

surface, to model nonlinearity in compression. Evaluation of the model through 
73

correlations with wood post data (static and dynamic) is given in Murray and Reid.(2)

Clear wood data for pine, available as default material properties, are given in Green 

and Kretschmann.(14) 

1. EL is the normal modulus of the undamaged wood parallel to the grain. The subscript 

L refers to the longitudinal direction of the wood. Typical units for moduli are GPa, 

MPa, or lbf/inch2. 

2. ET is the normal modulus of the undamaged wood perpendicular to the grain. The 

subscript T refers to the tangential direction, which, for a transversely isotropic 

material, is modeled the same as the ‘R’ radial direction. 

 

3. GLT is the shear modulus of the undamaged wood parallel to the grain. 

4. GTR is the shear modulus of the undamaged wood perpendicular to the grain. 

5. νLT is the major Poisson’s ratio parallel to the grain. A transversely isotropic material 

has three ratios, only one of which is independent. The minor Poisson’s ratio parallel 

to the grain, νTL, is calculated internally within the model as νTL = νLT (ET / EL). The 

Poisson’s ratio perpendicular to the grain (in the isotropic plane) is calculated 

internally as νTR = νRT = (ET – 2GTR) / 2GTR. 

6. Inputs XT, XC, and S|| are the strengths, parallel to the grain, in uniaxial tensile stress, 

uniaxial compressive stress, and pure shear stress, respectively. Together, these 

three strengths form an irregular elliptical surface for modeling failure or yielding 

parallel to the grain. These three strengths form the ultimate yield surface. The 

compressive strength is also scaled back to form the initial yield surface when 

prepeak hardening is requested. Yielding is initiated once the stress state exceeds 

the yield surface. In this case, a plasticity algorithm with associated flow returns the 

stress state to the yield surface (normal return). Typical units for strength are GPa, 

MPa, or lbf/inch2

. 

7. Inputs YT, YC, and S⊥ are the strengths, perpendicular to the grain, in uniaxial tensile 

stress, uniaxial compressive stress, and pure shear stress, respectively. The yield 

surface for the perpendicular modes is separate from the yield surface for the 

parallel modes. 

8. Gf I || and Gf II || are the fracture energies in uniaxial tensile stress and pure shear 

stress parallel to the grain. The fracture energy is the area under the stressdisplacement curve as it softens from peak stress to zero stress. The softening 

model incorporates fracture energy and element size (passed to the wood model 

internally) in order to regulate mesh-size sensitivity. The result is a softening 

response that is independent of element size. This is accomplished by keeping the 

fracture energy constant from element to element, regardless of size. Typical units 

are GPa-mm, MPa-cm, or lbf/inch2

-inch. 
74

9. B is the parallel-to-the-grain softening parameter. It sets the shape of the softening 

curve plotted as stress-strain or stress-displacement. B is unitless. A typical value is 

B = 30. 

10. dmax|| is the maximum damage that can accumulate parallel to the grain. Damage 

accumulation is based on an undamaged elastic strain energy norm formulated from 

the parallel strains. No damage threshold is required as input. The threshold is 

internally calculated and stored once the parallel yield criterion is satisfied. Thus, 

damage initiates with plastic yielding and accumulates if the energy norm of the 

current time step exceeds the energy norm of the previous time step. dmax|| is 

unitless and ranges from 0 (no damage) to 1 (100-percent damaged). A value of 

dmax|| = 0.9999 is typical. 

11. Gf I ⊥ and Gf II ⊥ are the fracture energies in tension and shear perpendicular to the 

grain. Typical units are GPa-mm, MPa-cm, or lbf/inch2-inch. 

12. D is the perpendicular-to-the-grain softening parameter. It sets the shape of the 

softening curve plotted as stress-strain or stress-displacement. D is unitless. A 

typical value is D = 30. 

13. dmax⊥ is the maximum damage that accumulates perpendicular to the grain. 

Damage accumulation is based on an undamaged elastic strain energy norm 

formulated from the perpendicular strains. Perpendicular damage accumulates 

separately from parallel damage. dmax⊥ is unitless. A value of dmax⊥ = 0.99 is 

typical. 

14. Inputs η||, ηc||, and n|| are rate-effect parameters that increase the parallel strengths 

of wood as a function of strain rate via shifted surface formulations. Dynamic 

strength is equal to the static strength plus an overstress. The overstress is equal to 

EL ε

( − ) & 1 n|| η|| in tension, GLT ε( − ) & 1 n || η|| in shear, and EL ε( − ) & 1 n || ηc|| in compression, 

where ε& is the effective strain rate. The products ε

( − ) & 1 n || η|| and ε( − ) & 1 n || ηc|| are 

unitless. Values for n|| range from 0 to less than 1, with a typical value of n|| = 0.1. 

15. Inputs η⊥, ηc⊥, and n⊥ are rate-effect parameters that increase the perpendicular 

strengths of wood as a function of strain rate. Strength increases for the 

perpendicular stresses are modeled separately from strength increases for the 

parallel stresses. The overstress is equal to ET ε

( − ⊥ ) & 1 n η⊥ in tension, GTR ε( − ⊥ ) & 1 n η⊥

in shear, and ET ε

( − ⊥ ) & 1 n ηc⊥ in compression. 

16. Inputs N|| and c|| are the hardening parameters for modeling prepeak nonlinearity in 

compression parallel to the grain. The hardening formulation translates an initial 

yield surface until it coincides with the ultimate yield surface that is formulated from 

the parallel strengths. The parameter N|| determines the location of the initial yield 
75

surface. The parameter c|| determines the rate of translation. Both parameters are 

unitless. Typical values are N|| = 0.5 with c|| = 1000. 

17. Inputs N⊥ and c⊥ are the hardening parameters for modeling prepeak nonlinearity in 

compression perpendicular to the grain. Perpendicular hardening is formulated 

separately from parallel hardening. 

18. Input GHARD is an additional hardening parameter that allows each translating yield 

surface to surpass the ultimate yield surface. A positive value will produce continued 

hardening in excess of the yield strengths. A zero value will produce perfect 

plasticity in compression. A small positive value (e.g., GHARD = 0.02) is 

recommended if computational difficulties are suspected because of perfect 

plasticity. 

2.3 FITTING THE MODEL TO THE DATA 

1. Stiffnesses (EL, ET, GLT, and GTR) and the major Poisson’s ratio (νLT) are directly 

measured from the elastic portion of the stress-strain curves or are selected from 

data tabulated in the literature. If the parallel stiffness EL is known, then one can 

estimate the other elastic stiffnesses from tables documented in the literature for 

softwoods (e.g., pine).(16) Typically, ET is 5 to 15 percent of EL. Care must be taken 

to measure or select the correct Poisson’s ratio. The major ratio, νLT, is typically 

greater than 0.1 and is about 10 times larger than the minor ratio, νTL. 

2. Strengths (XT, XC, YT, YC, S||, and S⊥) are obtained from measurements of peak/yield 

stress from stress histories or are selected from data tabulated in the literature. 

Typically, the parallel tensile strength of pine is 30 to 50 times greater than the 

perpendicular tensile strength. The parallel compressive strength is four to five times 

greater than the perpendicular compressive strength. The parallel shear strength is 

about 10 to 15 percent of the parallel tensile strength. Data for perpendicular shear 

strength vary. One source suggests that the perpendicular shear strength, S⊥, is 140 

percent of the parallel shear strength, S||.

(7) Another source suggests that the 

percentage is only 18 to 28 percent.(18)

3. The prepeak hardening parameters (N||, c||, N⊥, and c⊥) are derived from the 

nonlinear portion of compressive stress-strain curves measured both parallel and 

perpendicular to the grain. Hardening is modeled separately for the parallel and 

perpendicular modes. The parameters N|| and N⊥ determine the onset of 

nonlinearity. If the user wants prepeak nonlinearity to initiate at 70 percent of the 

yield stress, then the user would input N = 0.3 (for 30 percent), which is derived from 

N = 1 – 0.7 (100 percent – 70 percent). The parameters c|| and c⊥ set the amount of 

nonlinearity and are selected through iteration: Pick a value for c, plot the results 

from a single-element simulation, and compare with the test data. Typical values of c

are between 100 and 1000 (unitless). Gradual hardening is accomplished with 

smaller values of c. Rapid hardening is accomplished with larger values of c. If no 
76

data are available, prepeak hardening can be neglected (giving a linear stress-strain 

curve to yield) by setting both values of N equal to zero. 

4. The fracture energy (Gf I ||, Gf II ||, Gf I ⊥, or Gf II ⊥) is the area under the stressdisplacement curve as it softens from peak stress to zero stress. The best approach 

is to measure fracture energy directly in both direct tension and simple shear. An 

alternative approach is to measure fracture intensity and convert it to fracture 

energy. 

5. The postpeak softening parameters (B and D) are derived from the softening portion 

of the stress-strain or stress-displacement curves and are set in conjunction with the 

fracture energies (Gf I ||, Gf II ||, Gf I ⊥, and Gf II ⊥) and the damage parameters (dmax||

and dmax⊥). The parameters B (for parallel modes) and D (for perpendicular modes) 

set the shape of the softening portion of the stress-strain or stress-displacement 

curves once the fracture energy is selected. Although each parameter (B or D) is 

intended to be simultaneously fit to both tension and shear data, shear data are 

rarely available. The procedure is iterative: Pick a value for D or B, plot the results 

from a single-element simulation, and compare with the test data. Typical values for 

B and D are 10 to 50 (unitless). The smaller the value, the more gradual the initial 

softening. 

6. The postpeak damage parameters (dmax|| and dmax⊥) are derived from the final 

softening portion of the stress-strain or stress-displacement curves and are set in 

conjunction with the fracture energies (Gf I ||, Gf II ||, Gf I ⊥, and Gf II ⊥) and the softening 

parameters (B and D). The parameters dmax|| (for parallel modes) and dmax⊥ (for 

perpendicular modes) set the maximum damage that can accumulate. Values 

between 0 and 1 are allowed. A value of 0 neglects softening, while a value of 1 

models complete softening to 0 stress and stiffness. A typical value for parallel 

damage is dmax|| = 0.9999. Erosion is automatically modeled when parallel damage 

exceeds 0.99. Set dmax|| > 0.99 to model erosion at high damage levels. Set 

dmax|| ≤ 0.99 to bypass erosion (not recommended because this may cause mesh 

entanglement and element inversion). A typical value for perpendicular damage is 

dmax⊥ = 0.99. Values greater than 0.99 are not recommended because elements 

have a tendency to tangle and invert with near zero stiffness, and erosion is not 

automatically modeled with perpendicular damage. Values less than 0.99 should be 

input if mesh entanglement and element inversion are a problem. 

7. The rate-effect parameters (η||, ηc||, n||, η⊥, ηc⊥, and n⊥) are obtained from fits to 

strength versus strain-rate measurements or are selected from data tabulated in the 

literature. Separate measurements must be made parallel and perpendicular to the 

grain to separately fit the parameters parallel (η||, ηc||, and n||) and perpendicular (η⊥, 

ηc⊥, and n⊥) to the grain. Separate measurements may also be made in 

tension/shear versus compression to separately fit the tensile/shear (η|| or η⊥) and 

compression (ηc|| or ηc⊥) parameters. If separate tension/shear versus compression 

measurements are not available, then the following ratios are recommended: η||/ηc|| 

= η⊥/ηc⊥ = QT / QC. Each set of parameters is obtained from dynamic strength 
77

measurements at two different strain rates (because there are two parameters), as 

well as from the static strength (zero strain rate). At each of two strain rates, the 

difference between the dynamic and static strength is Cε& η, where C is the stiffness 

(EL, ET, GLT, or GTR), ε& is the effective strain rate, and η is the effective fluidity 

parameter equal to η = η0/ε& n. Here, η0 (η|| or ηc||, or η⊥ or ηc⊥) and n (n|| or n⊥) are 

the two parameters to be fit (C and ε& are known). The two parameters are found by 

simultaneously solving two equations with two unknowns (see equation 60 in the 

Model Formulation section). Typical clear wood values are n|| = 0.107, with η|| = ηc|| = 

0.0045 for parallel rate effects, and n⊥ = 0.104, with η⊥ = ηc⊥ = 0.0962 for 

perpendicular rate effects (for time in milliseconds). Typical graded wood values are 

n|| = 0.107, with η|| = 0.0045QT and ηc|| = 0.0045QC (parallel), and n⊥ = 0.104, with η⊥

= 0.0962QT and ηc⊥ = 0.0962QC (perpendicular). 

2.4 MODEL FORMULATION 

Elastic Constitutive Equations 

The general constitutive relationship for an orthotropic material, written in terms of the 

principal material directions, is:(16) 

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

=

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

⎡

6

5

4

3

2

1

66

55

44

13 23 33

12 22 23

11 12 13

6

5

4

3

2

1

0 0 0 0 0 2

0 0 0 0 2 0

0 0 0 2 0 0

0 0 0

0 0 0

0 0 0

ε

ε

ε

ε

ε

ε

σ

σ

σ

σ

σ

σ

C

C

C

C C C

C C C

C C C

(101)

The subscripts 1, 2, and 3 refer to the longitudinal, tangential, and radial stresses and 

strains (σ1 = σ11, σ2 = σ22, σ3 = σ33, ε1 = ε11, ε2 = ε22, and ε3 = ε33), respectively. The 

subscripts 4, 5, and 6 are in a shorthand notation that refers to the shearing stresses 

and strains (σ4 = σ12, σ5 = σ23, σ6 = σ13, ε4 = ε12, ε5 = ε23, and ε6 = ε13). As an alternative 

notation for wood, it is common to substitute L (longitudinal) for 1, R (radial) for 2, and T

(tangential) for 3. The components of the constitutive matrix, Cij, are listed here in terms 

of the nine independent elastic constants of an orthotropic material: 
78

C11 = E11(1−ν 23ν 32 )/ Δ (102)

C22 =E22 (1−ν 31ν 13 ) / Δ (103)

C33 =E33 (1−ν 12ν 21) / Δ (104)

C12 =(ν 21+ν 31ν 23 )E11 / Δ (105)

C13 =(ν 31+ν 21ν 32 )E11 / Δ (106)

C23 =(ν 32 +ν 12ν 31)E22 / Δ (107)

C44 =G12 (108)

C55 =G23 (109)

C66 =G13 (110)

Δ =1−ν 12ν 21−ν 23ν 32−ν 31ν 13−2ν 21ν 32ν 13 (111)

The following identity, relating the dependent (minor Poisson’s ratios ν21, ν31, and ν32) 

and independent elastic constants, is obtained from symmetry considerations of the 

constitutive matrix: 

= for i, j =1, 2, 3

E E j

ji

i

ν ij ν

(112)

One common assumption is that wood materials are transversely isotropic. This means 

that the properties in the tangential and radial directions are modeled the same (i.e., 

E22 = E33, G12 = G13, and ν12 = ν13. This reduces the number of independent elastic 

constants to five (E11, E22, ν12, G12, and G23). Furthermore, the Poisson’s ratio in the 

isotropic plane, ν23, is not an independent quantity. It is calculated from the isotropic 

relationship ν = (E – 2G) / 2G, where E = E22 = E33 and G = G23. Transverse isotropy is 

a reasonable assumption because the difference between the tangential and radial 

properties of wood (particularly southern yellow pine and Douglas fir) is small in 

comparison with the difference between the tangential and longitudinal properties. 

Yield Surfaces 

The yield surfaces parallel and perpendicular to the grain are formulated from six 

ultimate strength measurements obtained from uniaxial and pure-shear tests on wood 

specimens: 

XT Tensile strength parallel to the grain 

XC Compressive strength parallel to the grain 
79

YT Tensile strength perpendicular to the grain 

YC Compressive strength perpendicular to the grain 

S|| Shear strength parallel to the grain 

S⊥ Shear strength perpendicular to the grain 

The formulation is based on the work of Hashin.(27) 

Parallel Modes 

For the parallel modes, the yield criterion is composed of two terms involving two of the 

five stress invariants of a transversely isotropic material. These invariants are 1 = σ11 I

and 2

13

2 I4 = σ12 +σ . This criterion predicts that the normal and shear stresses are 

mutually weakening (i.e., the presence of shear stress reduces the strength below that 

measured in the uniaxial stress tests). Yielding occurs when f|| ≥ 0, where: 

⎩

⎨

⎧

<

> − = + = +

0

0

1 ( )

C 11

T 11

2

||

2

13

2

12

2

2

11

|| σ

σ σ σ σ

X for

X for X

X S

f (113)

Perpendicular Modes 

For the perpendicular modes, the yield criterion is also composed of two terms involving 

two of the five stress invariants of a transversely isotropic material. These invariants are 

I2 = σ22 + σ33 and 22 33

2 I3 = σ 23 −σ σ . Yielding occurs when f⊥ ≥ 0, where: 

⎩

⎨

⎧

+ <

+ > − = − +

+ =

⊥

⊥ 0

0

1 ( ) ( )

C 22 33

T 22 33

2

22 33

2

23

2

2

22 33

σ σ

σ σ σ σ σ σ σ

Y for

Y for Y

Y S

f (114)

Each yield criterion is plotted in three dimensions in figure 26 in terms of the parallel and 

perpendicular stresses. Each criterion is a smooth surface (no corners). 

Plastic Flow 

The plasticity algorithms limit the stress components once the yield criteria in equations 

111 and 112 are satisfied. This is done by returning the trial elastic stress state back to 

the yield surface. The stress and strain tensors are partitioned into elastic and plastic 

parts. Partitioning is done with a return mapping algorithm that enforces the plastic 

consistency condition. Separate plasticity algorithms are formulated for the parallel and 

perpendicular modes by enforcing separate consistency conditions. No input 

parameters are required. 
80

s L

s LT

s LR

(a) Parallel modes 

s TR

s T

s R

(b) Perpendicular modes 

Figure 26. Yield criteria for wood produce smooth surfaces in stress space. 

σLR

σLT

σL

σR

σT

σTR
81

Hardening 

Wood exhibits prepeak nonlinearity in compression parallel and perpendicular to the 

grain. Separate translating yield surface formulations are modeled for the parallel and 

perpendicular modes, which simulate gradual changes in moduli. Each initial yield 

surface hardens until it coincides with the ultimate yield surface, as shown in figure 27. 

The initial location of the yield surface determines the onset of plasticity. The rate of 

translation determines the extent of nonlinearity. 

Square Root of Parallel Shear Invariant

XT

-XC

-(1-NII)XC

SII

Tension

Compression

Longitudinal Stress

Ultimate Yield Surface

Initial Yield Surface

XC

(1-NII)XC

NIIXC

Initial Yield Strength

Ultimate Yield Strength

Increasing Rate of Translation

CII small

CII large

Strain

 | Longitudinal Stress |

(a) Initial and ultimate yield surfaces (b) Stress-strain behavior 

Figure 27. Prepeak nonlinearity in compression is modeled with translating yield 

surfaces that allow user to specify hardening response. 

For each mode (parallel and perpendicular), the user inputs two parameters―the initial 

yield surface location in uniaxial compression, N, and the rate of translation, c. If the 

user wants prepeak nonlinearity to initiate at 70 percent of the peak strength, then the 

user would input N = 0.3, so that 1 – N = 0.7. If the user wants rapid hardening, then a 

large value of c is input (e.g., c = 1000). If the user wants gradual hardening, then a 

small value of c is input (e.g., c = 10). 

The state variable that defines the translation of the yield surface is known as 

backstress and is denoted by αij. Hardening is modeled in compression, but not shear, 

so the only backstress required for the parallel modes is α11. The value of the 

backstress is α11 = 0 upon initial yield and α11 = –N|| Xc at ultimate yield (in uniaxial 

compression). The maximum backstress occurs at ultimate yield and is equal to the 

total translation of the yield surface in stress space. The backstress components 

required for the perpendicular modes are α22 and α33. The value of the backstress sum 

is α22 + α33 = 0 upon initial yield and α22 + α33 = –N⊥ Yc at ultimate yield (biaxial 

compression without shear). 

Square Root of Parallel Shear Invariant Strain 

Initial Yield Strength 

Increasing Rate of Translation 

C|| small

C|| large

(1-N||)Xc

Longitudinal Stress 

Longitudinal Stress 

-(1-N||)Xc 

-Xc

Compression 

Tension 

Initial Yield Surface 

N||Xc
82

Damage 

Separate damage formulations are modeled for the parallel and perpendicular modes. 

These formulations are loosely based on the work of Simo and Ju.(21) If failure occurs in 

the parallel modes, then all six stress components are degraded uniformly. This is 

because parallel failure is catastrophic and will render the wood useless. If failure 

occurs in the perpendicular modes, then only three perpendicular stress components 

are degraded. This is because perpendicular failure is not catastrophic; the wood may 

continue to carry load in the parallel direction. Based on these assumptions, the 

following degradation model is implemented: 

max( ( ), ( )) dm = d τ || d τ ⊥ (115)

( ) || || d = d τ (116)

11 || 11 σ = (1− d )σ (117)

22 22 σ = (1− dm )σ (118)

33 33 σ = (1− dm )σ (119)

12 || 12 σ = (1− d )σ (120)

13 || 13 σ = (1− d )σ (121)

23 23 σ = (1− dm )σ (122)

Here, each scalar damage parameter, d, transforms the stress tensor associated with 

the undamaged state, σ ij , into the stress tensor associated with the damaged state, σij. 

The stress tensor σ ij is calculated by the plasticity algorithm (including viscoplasticity) 

prior to application of the damage model. Each damage parameter ranges from zero for 

no damage to approaching unity for maximum damage. Thus, 1 – d is a reduction factor 

associated with the amount of damage. Each damage parameter evolves as a function 

of a strain energy-type term. Mesh-size dependency is regulated via a length scale 

based on the element size (cube root of volume). Damage-based softening is brittle in 

tension, less brittle in shear, and ductile (no softening) in compression, as demonstrated 

in figure 28. 

Element erosion occurs when an element fails in the parallel mode and the parallel 

damage parameter exceeds d|| = 0.99. Elements do not automatically erode when an 

element fails in the perpendicular mode. A flag is available that, when set, allows 

elements to erode when the perpendicular damage parameter exceeds d⊥ = 0.99. 

Setting this flag (IFAIL) is not recommended unless excessive perpendicular damage is 

causing computational difficulties. 
83

(a) Tensile softening 

(b) Shear softening 

(c) Compressive yielding 

Figure 28. Softening response modeled for parallel modes 

of southern yellow pine. 
84

Rate Effects 

Data available in the literature for pine indicate that dynamic strength enhancement is 

more pronounced in the perpendicular direction than in the parallel direction.(16)

Therefore, separate rate-effect formulations are modeled for the parallel and 

perpendicular modes. The formulations increase strength with increasing strain rate by 

expanding each yield surface: 

Parallel 

XT 

DYNAMIC

 = XT + EL || ε

( − ) & 1 n|| η|| Tension 

XC 

DYNAMIC

 = XC + EL || ε

( − ) & 1 n|| ηc|| Compression 

S||

DYNAMIC

 = S|| + GLT || ε

( − ) & 1 n|| η|| Shear 

(123)

Perpendicular 

YT 

DYNAMIC

 = YT + ET ⊥

⊥ ε

( − ) & 1 n η⊥ Tension 

YC 

DYNAMIC

 = YC + ET ⊥

⊥ ε

( − ) & 1 n ηc⊥ Compression 

S⊥

DYNAMIC

 = S⊥ + GTR ⊥

⊥ ε

( − ) & 1 n η⊥ Shear 

(124)

where: 

X and Y = static strengths 

X DYNAMIC and Y DYNAMIC = dynamic strengths 

Cε& (1-n)

η = excess stress components 

The excess stress components depend on the value of the fluidity parameter, η; the 

stiffness, C; and the effective strain rate, ε& . When rate effects are requested via the flag 

IRATE = 1, the dynamic strengths are used in place of the static strengths in the yield 

surface formulations. 
85

Setting n > 0 allows the user to model a nonlinear variation in dynamic strength with 

strain rate. Setting n = 0 allows the user to model a linear variation in dynamic strength 

with strain rate. 

2.5 WOOD MODEL OUTPUT 

One example output file is given here for completeness. It gives the output format for 

the PINE and FIR options. If the <BLANK> option is selected, then the lines for wood 

type through units are omitted. The example values printed are for saturated grade 1 

pine at room temperature in units of MPa, cm, and ms. 

Printed for PINE and FIR options only. 

Default Material Property Selection: 

 WOOD Wood type = 0.0000e+00 

 EQ. 0: Southern yellow pine (default) 

 EQ. 1: Douglas fir 

 MC Moisture content (%) = 2.3000e+01 

 TEMP Temperature = 2.0000e+01 

QT Tension/shear quality factor = 0.4700e+00 

QC Compression quality factor = 0.6300e+00 

 Quality factor option = 0.0000e+00 

 EQ. 0: Grades 1, 1D, 2, 2D (default) 

 EQ. -1: DS-65 or SEL STR 

 EQ. -2: Clear wood 

 GT. 0: User-defined 

 Apply quality factors perpendicular = 1.0000e+00 

 EQ. 0: Yes (default) 

 EQ. 1: No (use clear wood strengths) 

 UNITS Units = 1.0000e+00 

 EQ. 0: GPa, mm, ms, kg/mm3

, kN (default) 

 EQ. 1: MPa, mm, ms, g/mm3, N 

 EQ. 2: MPa, mm, s, Mg/mm3, N 

 EQ. 3: lbf/inch2, inch, s, lb-s2/inch4, lbf
86

Printed for all options. 

Control Parameters: 

 NPLOT Plot parameter = 0.0000e+00 

 EQ. 0: Maximum of parallel and perpendicular damage (default) 

 EQ. 1: Perpendicular damage 

 ITERS Number of plasticity algorithm iterations = 1.0000e+00 

 IRATE Rate-effect option = 1.0000e+00

EQ. 0: Rate effects off (default) 

 EQ. 1: Rate effects on 

 GHARD Perfect plasticity override = 0.0000e+00 

 EQ. 0: Perfect plasticity (default) 

 GT. 0: Continuous hardening 

 IFAIL Erode with perpendicular damage = 0.0000e+00 

 EQ. 0: No (default) 

 EQ. 1: Yes (not recommended)

Stiffness: 

EL Parallel normal modulus = 1.1350e+04 

ET Perpendicular normal modulus = 2.4680e+02 

GLT Parallel shear modulus = 7.1520e+02 

GTR Perpendicular shear modulus = 8.7510e+01 

νLT Parallel major Poisson’s ratio = 0.1568e+00 

Strength: 

XT Parallel tensile strength = 4.0250e+01 

XC Parallel compressive strength = 1.3330e+01 

YT Perpendicular tensile strength = 2.0500e+00 

YC Perpendicular compressive strength = 4.0820e+00 

S|| Parallel shear strength = 6.2410e+00 

S⊥ Perpendicular shear strength = 1.2730e+01 

 

Damage: 

Gf I || Parallel fracture energy in tension = 4.2660e+01 

Gf II || Parallel fracture energy in shear = 8.8260e+01 

B Parallel softening parameter = 3.0000e+01 

dmax|| Parallel maximum damage = 0.9999e+00 

Gf I ⊥ Perpendicular fracture energy in tension = 4.0090e−01 

Gf II ⊥ Perpendicular fracture energy in shear = 8.2950e−01 
87

D Perpendicular softening parameter = 3.0000e+01 

dmax⊥ Perpendicular maximum damage = 0.9900e+00 

Rate Effects: 

η|| Parallel fluidity parameter in tension/shear = 0.0045e+00 

ηc|| Parallel fluidity parameter in compression = 0.0045e+00 

n|| Parallel power = 0.1070e+00 

η⊥ Perpendicular fluidity parameter in tension/shear = 0.0962e+00 

ηc⊥ Perpendicular fluidity parameter in compression = 0.0962e+00 

n⊥ Perpendicular power = 0.1040e+00 

Hardening: 

 N|| Parallel hardening initiation = 0.5000e+00 

 c|| Parallel hardening rate = 4.0000e+02 

 N⊥ Perpendicular hardening initiation = 0.4000e+00 

c⊥ Perpendicular hardening rate = 1.0000e+02 

Orientation Options: 

AOPT = 2.0000e+00 

 Orientation parameter #1 = 0.0000e+00 

 Orientation parameter #2 = 0.0000e+00 

 Orientation parameter #3 = 1.0000e+00 

 Orientation parameter #4 = 1.0000e+00 

 Orientation parameter #5 = 0.0000e+00 

 Orientation parameter #6 = 0.0000e+00

89

3. EXAMPLES MANUAL 

This section contains two sets of example problems that the user can review to become 

familiar with the use of the wood material model. The first set of problems are singleelement simulations in tension and compression, both parallel and perpendicular to the 

grain. The second problem is a bogie impact into a wooden post. 

3.1 SINGLE-ELEMENT SIMULATIONS 

Wood material model input is given in figure 29 for default pine parameters, and in 

figure 30 for user-specified input parameters. A complete input file, with nodes and 

elements, is given in appendix G. This file is for tensile loading of a single element 

parallel to the grain. To convert to compressive loading, change the sign of the ordinate 

under *DEFINE_CURVE at the bottom of the file. To convert to loading perpendicular to 

the grain, change the orientation vectors to [A1, A2, A3] = [1, 0, 0] and [D1, D2, D3] 

= [0, 0, 1]. 

Figure 29. Example wood model input for selection of default input 

parameter (option MAT_WOOD_PINE). 

Single-element stress-strain results are shown in figure 31 for saturated, roomtemperature clear wood pine. These results can be achieved using either the default 

input shown in figure 29 or the user-specified input shown in figure 30. Note that the 

peak strengths attained match the input strength parameters listed in figure 30. These 

are XT = 85.2 MPa, XC = 21.2 MPa, YT = 2.05 MPa, and YC = 4.08 MPa. The results are 

plotted with LS-POST as cross-plots of element z stress versus node 8 displacement. 

As additional exercises, the user can vary the moisture content (MC), temperature 

(TEMP), grade (QT, QC, and IQUAL), units, hardening (GHARD), and rate effects 

(IRATE with rate of loading) to examine how the default behavior of pine varies with 

these quantities. 
90

Figure 30. Example wood model input for user specification 

of input parameters (option MAT_WOOD). 
91

A Tension

B Compression

−0.05 0 0.05 0.1

−20 

0 

20 

40 

60 

80 

Z−displacement (cm) 

Z−stress (MPa) 

A

(a) Parallel 

 

A Tension

B Compression

−0.1 −0.05 0 0.05

−4 

−3 

−2 

−1 

0 

1 

2 

Z−displacement (cm) 

Z−stress (MPa)

A

A

A

A

B

 (b) Perpendicular 

Figure 31. Example single-element stress-strain results for clear wood pine. 
92

3.2 BOGIE IMPACT SIMULATION 

The wood model input is similar to that shown in figure 29, although units of GPa 

(UNITS = 0 with RO = 6.73-e07) and rate effects (IRATE = 1) are requested. In addition, 

the post is grade 1, rather than clear wood. The post is saturated along the entire 

length. The post is modeled with simplified boundary conditions to attain a fast run time. 

Bogie-post output results are shown in figures 32 and 33. The post breaks just below 

ground level. Damage occurs in the impact regime and in the region of post breakage. 

The post displacement is linear with time. The post cross-sectional force (between the 

impact point and ground level) reaches a peak of approximately 47 MPa within the first 

5 ms before decaying, in an oscillatory manner, back down to zero. The results are 

plotted with LS-POST. 

Figure 32. Deformed configuration of post at 40 ms, including fringes of damage. 
93

Node ID

A 22

0 10 20 30

0 

100 

200 

300 

Time (msec)

X−displacement (mm)

(a) Deflection 

 

Section ID

A 1

0 10 20 30

0 

10 

20 

30 

40 

Time (msec)

X−Force (kN)

A

(b) Force 

Figure 33. Post deflection and cross-sectional force histories.

Time (ms) 

Time (ms) 

95

APPENDIX A. MEASURED VARIABILITY OF SOUTHERN 

YELLOW PINE 

 (a) 4-percent moisture content (b) 8-percent moisture content 

 

(c) 12-percent moisture content (d) Saturated 

Figure 34. Measured load displacement curves of southern yellow pine 

exhibit variability in tension parallel to the grain. 

Source: Forest Products Laboratory (14). 
96

 

(a) 4-percent moisture content (b) 8-percent moisture content 

 

(c) 12-percent moisture content (d) Saturated 

Figure 35. Measured load displacement curves of southern yellow pine 

exhibit variability in compression perpendicular to the grain. 

Source: Forest Products Laboratory (14).
97

APPENDIX B. QUADRATIC EQUATIONS FIT TO MOISTURE 

CONTENT DATA 

(a) Data (b) Fit 

Figure 36. Effect of moisture content on tensile modulus parallel to the grain. 

Source: Forest Products Laboratory.(14)

 

(a) Data (b) Fit 

Figure 37. Effect of moisture content on tensile modulus perpendicular 

to the grain. Source: Forest Products Laboratory.(14)
98

 

(a) Data (b) Fit 

Figure 38. Effect of moisture content on tensile strength parallel 

to the grain. Source: Forest Products Laboratory.(14) 

(a) Data (b) Fit 

Figure 39. Effect of moisture content on tensile strength perpendicular 

to the grain. Source: Forest Products Laboratory.(14)
99

 

(a) Data (b) Fit 

Figure 40. Effect of moisture content on compressive strength parallel 

to the grain. Source: Forest Products Laboratory.(14)

 

(a) Data (b) Fit 

Figure 41. Effect of moisture content on compressive strength perpendicular 

to the grain. Source: Forest Products Laboratory.(14)
100

(a) Data (b) Fit 

Figure 42. Effect of moisture content on shear strength parallel 

to the grain. Source: Forest Products Laboratory.(14)

No measurements were made for shear strength perpendicular to the grain. 
101

(a) Data (b) Fit 

Figure 43. Effect of moisture content on mode I fracture intensity. 

Source: Forest Products Laboratory.(14)

(a) Data (b) Fit 

Figure 44. Effect of moisture content on mode II fracture intensity. 

Source: Forest Products Laboratory.(14)

103

APPENDIX C. ANALYTICAL FORM OF CANDIDATE 

FAILURE CRITERIA 

C.1 REVIEW OF CRITERIA 

Numerous failure criteria available in the literature were reviewed and evaluated for 

modeling the yield strength of wood and composite materials. To our knowledge, a 

validated three-dimensional theory for modeling wood is not documented in the 

literature. However, numerous failure criteria have been documented for fiber-reinforced 

plastic (FRP) composites (e.g., see the survey by Nahas(28)). Composites are similar to 

wood because they are transversely isotropic materials with distinct failure modes in the 

parallel (fiber) and perpendicular (transverse fiber) directions. Therefore, many of the 

criteria originally developed for modeling composites were evaluated as candidates for 

modeling wood. 

The functional form of each failure criterion that was evaluated is given in section C.2. 

These include one limit and six interactive criteria. Both orthotropic and transversely 

isotropic criteria are reported. All criteria are stress-based criteria. The stresses are 

transformed to the principal material axes (L-T-R axes) before application of the failure 

criteria. Strain-based criteria were not evaluated because failure strains are not reported 

in the literature for wood. One cannot derive failure strains from stresses if the stressstrain behavior is nonlinear, as it is for wood in compression. A brief summary of each 

criterion is given here: 

Maximum Stress (commonly applied limit theory): Failure occurs when any 

component of stress exceeds its corresponding strength. 

Tsai-Wu (tensor polynomial theory that was originally developed for anisotropic 

materials): It contains linear and quadratic stress terms. Seven coefficients must be 

defined for transversely isotropic applications. The noninteraction coefficients are 

determined from measured uniaxial and pure shear strengths. By noninteraction, one 

means terms that contain one component of stress (e.g., F1σ11). The interaction terms 

are determined from measured biaxial strengths. By interaction, one means terms that 

have two or more components of stress multiplied together (e.g., F12σ11σ22). 

Hoffman: Hoffman extended Hill’s distortional energy criterion for orthotropic materials 

to account for different strengths in tension and compression. The criterion contains 

linear and quadratic stress terms. Six coefficients are determined from uniaxial stress 

and pure shear tests. Biaxial strengths are not needed. 

Norris: Norris developed three yield criteria for mutually orthogonal planes. Each 

criterion contains quadratic stress terms (no linear terms). Nine coefficients are 

determined from uniaxial and pure shear tests. Tensile strengths are used when the 

corresponding stresses are tensile. Compressive strengths are used when the 

corresponding stresses are compressive. 
104

Extended Yamada-Sun: Three yield criteria are reported for mutually orthogonal 

planes. Each criterion predicts that the normal and shear stresses are mutually 

weakening (the presence of shear stress reduces the strength below that measured in 

uniaxial stress tests). Nine coefficients are determined from uniaxial and pure shear 

tests. 

Hashin: Hashin formulated a quadratic stress polynomial in terms of the invariants of a 

transversely isotropic material. Separate formulations are identified for parallel and 

perpendicular modes by assuming that failure is produced by the normal and shear 

stresses acting on the failure plane. In addition, the parallel and perpendicular modes 

are subdivided into tensile and compressive modes. Assumptions include: (1) biaxial 

compressive strength perpendicular to the grain is much greater than the uniaxial 

compressive strength and (2) shear stress does not contribute to compressive failure 

parallel to the grain. All coefficients are determined from six uniaxial and shear 

strengths. 

Modified Hashin (extended form of Hashin’s criteria): More terms are retained in 

this modified form than in the original form because fewer assumptions are made 

regarding material behavior. All coefficients are determined from six uniaxial and shear 

strengths. 

C.2 FORM OF CRITERIA 

Here, we give the functional form of the various failure criteria that were evaluated in 

section C.1 for modeling the strength of wood. Both orthotropic and transversely 

isotropic criteria are reported. 

The following notation is used for orthotropic criteria: 

 XT Tensile strength in longitudinal direction 

XC Compressive strength in longitudinal direction 

YT Tensile strength in tangential direction 

YC Compressive strength in tangential direction 

ZT Tensile strength in radial direction 

ZC Compressive strength in radial direction 

Sxy Shear strength parallel to the grain in L-T plane 

Sxz Shear strength parallel to the grain in L-R plane 

Syz Shear strength perpendicular to the grain in T-R plane 

The following notation is used for transversely isotropic criteria: 

 XT Tensile strength parallel to the grain 

XC Compressive strength parallel to the grain 

YT Tensile strength perpendicular to the grain 

YC Compressive strength perpendicular to the grain 

S|| Shear strength parallel to the grain 

S⊥ Shear strength perpendicular to the grain 
105

Here, X, Y, and Z are the strengths in the longitudinal, tangential, and radial directions, 

respectively, and S is the shear strength. The subscripts T and C refer to the tensile and 

compressive components, respectively. 

Maximum Stress: This is one of the most common limit theories. Failure occurs when 

any component of stress in the principal material directions exceeds its corresponding 

strength. Its application to wood as an orthotropic material is: 

( 0) σ 11 ≥ XT σ 11 > (125)

( 0) σ 11 ≥ XC σ11 < (126)

( 0) σ 22 ≥YT σ 22 > (127)

( 0) σ 22 ≥YC σ 22 < (128)

( 0) σ 33 ≥ZT σ 33 < (129)

( 0) σ 33 ≥ZC σ 33 < (130)

XY ≥S σ 12 (131)

XZ ≥S σ 13 (132)

YZ ≥S σ 23 (133)

Nine independent modes of failure are predicted: tensile, compressive, and shear failure 

parallel to the grain; tensile, compressive, and shear failure in the tangential direction; 

and tensile, compressive, and shear failure in the radial direction. The number of failure 

modes reduces to six for a transversely isotropic material. 

Tsai-Wu: Tsai and Wu developed a stress tensor component polynomial theory as a 

failure criterion for anisotropic materials.(29) A reduced form of their criterion is applicable 

to transversely isotropic materials. The criterion contains both linear and quadratic 

stress terms. Failure occurs when the following equation is satisfied: 

2 ( ) 2 ( ) 1

( ) ( 2 ) ( )

22 33

2

12 11 22 11 33 23 23

2

13

2

66 12

2

23

2

33

2

22 22

2

1 11 2 22 33 11 11

+ + − ≥

+ + + + + + + + +

σ σ σ σ σ σ σ

σ σ σ σ σ σ σ σ σ

F F

F F F F F

(134)

Seven coefficients must be defined for wood modeled as a three-dimensional 

transversely isotropic material. Six coefficients (F1, F2, F11, F22, F23, and F66) are 

determined from uniaxial and shear tests on unidirectional specimens. Each of the 

coefficients F1, F2, F11, F22, and F23 includes contributions from both tensile and 

compressive strengths: 
106

XT XC

F 1 1

1 = − (135)

YT YC

F 1 1

2 = − (136)

XT XC

F 1

11 = (137)

YTYc

F 1

22 = (138)

2

||

66

1

S

F = (139)

S YtYc

F 1 2 2 2 23 = −

⊥

(140)

One coefficient, F12, must be determined from biaxial tests, a variety of which are 

available. Different biaxial tests produce different values of F12. The choice was to fit F12

to off-axis compression test data at 45 degrees: 

( ) ( ) 1 1 2 1 2 11 22 66 2 F12 = − F + F − F + F + F

σ σ

(141)

Here, σ is the biaxial strength measured in the principal material directions. It is equal to 

half the ultimate strength measured in off-axis tests at 45 degrees. If σ = XcYc / (Xc + Yc), 

then the Tsai-Wu model will be in agreement with Hankinson’s two-dimensional formula 

plotted in ultimate stress versus grain angle space. The failure envelope is a smooth 

surface in stress space. Only the onset of failure is predicted, not the mode of failure. 

Hoffman: Hill generalized von Mises’ distortional energy criterion for isotropic materials 

to include orthotropic materials.(30) Failure occurs when the following equation is 

satisfied: 

( ) ( ) ( ) 1 2

13

2

23

2

12

2

33 11

2

22 33

2 A σ 11 −σ 22 + B σ −σ + C σ −σ + Dσ + Eσ + Fσ = (142)

The six coefficients (A, B, D, E, F, and G) are determined from uniaxial stress and pure 

shear tests. However, Hill’s orthotropic criterion is not directly applicable to wood 

materials because it does not model different strengths in tension and compression. 

Hoffman modified Hill’s quadratic criterion by adding linear stress terms that take into 

account different strengths in tension and compression.(31) Failure occurs when the 

following equation is satisfied: 
107

2

13

2

23

2

12

2

33 11

2

22 33

2

11 22 A(σ −σ ) + B(σ −σ ) + C(σ −σ ) + Dσ + Eσ + Fσ

+Gσ 11 + Hσ 22 + Iσ 33 ≥ 1

(143)

XT XC YTYC ZT ZC

A

2

1

2

1

2

1 = + −

\

(144)

XT XC YTYC ZT ZC

B

2

1

2

1

2

1 =− + + (145)

XT XC YTYC ZT ZC

C

2

1

2

1

2

1 = − + (146)

2

1

Sxy

D = (147)

2

1

Sxz

E = (148)

2

1

Syz

F = (149)

XT XC

G 1 1 = − (150)

YT YC

H 1 1 = − (151)

ZT ZC

I 1 1 = − (152)

This criterion predicts a parabolic increase in strength with confining pressure. The nine 

coefficients are determined from uniaxial stress and pure shear tests. The criterion is 

readily simplified for materials with transversely isotropic strength values. One 

advantage of this criterion is that the interaction terms are not based on biaxial data, so 

it is easier to fit than the Tsai-Wu criterion. One disadvantage of this criterion (and the 

Tsai-Wu criterion) is that the onset of failure is predicted, but not the mode of failure. 

Norris: Tsai and Azzi simplified the Hill criterion to account for transverse isotropy and 

plane stress conditions of composite materials:(32) 

1 2

2

12

2

2

22

2

11 22

2

2

11 − + + ≥

X X Y S

σ σ σ σ σ (153)
108

Tsai showed that the criterion is applicable to composites with different properties in 

tension and compression. Tensile strengths are used when the corresponding stresses 

are tensile; compressive strengths are used when the corresponding stresses are 

compressive. Tsai also developed two additional equations for mutually orthogonal 

planes (similar to equation 153) for failure analysis of three-dimensional materials.(27) 

Similarly, Norris reports three yield criteria for mutually orthogonal planes. His criteria 

are similar to the Tsai-Azzi criteria except that the interaction terms are not biased 

toward one particular strength. In addition, he applied his criteria to wood materials, not 

composites. Norris’s criterion for the 1-2 (L-R) plane is as follows: 

1 2

2

12

2

2

11 22 22

2

2

11 − + + ≥

X XY Y S

σ σ σ σ σ (154)

Similar equations can be written for the 1-3 and 2-3 planes by proper interchange of 

subscripts. This criterion for modeling wood was evaluated by using tensile strengths 

when the corresponding stresses are tensile, and compressive strengths when the 

corresponding stresses are compressive. For each of the three equations, three 

combinations of tensile and compressive stresses are possible. Therefore, nine modes 

of failure are modeled. 

Extended Yamada-Sun: Yamada and Sun developed a plane stress criterion for 

composites that is a degenerative form of the Tsai-Azzi and Norris criteria.(33) Failure 

occurs when the following equation is satisfied: 

1 2

2

12

2

2

11 + ≥

X S

σ σ (155)

This criterion predicts that the normal and shear stresses are mutually weakening (the 

presence of shear stress reduces the strength below that measured in uniaxial stress 

tests). We extended this concept to three dimensions for application to wood as either 

an orthotropic or transversely isotropic material: 

1 2

2

13

2

2

12

2

2

11 + + ≥

X Sxy Sxz

σ σ σ

(156)

1 2

2

23

2

2

12

2

2

22 + + ≥

xy yz Y S S

σ σ σ

(157)
109

1 2

2

23

2

2

13

2

2

33 + + ≥

Z Sxz S yz

σ σ σ (158)

Tensile strengths are used when the corresponding stresses are tensile; compressive 

strengths are used when the corresponding stresses are compressive. Six independent 

modes of failure are predicted: tensile and compressive failure in the longitudinal 

direction, tensile and compressive failure in the tangential direction, and tensile and 

compressive failure in the radial direction. The application of shear stress contributes to 

failure in each of these modes and is mutually weakening. 

Hashin: The Tsai-Wu and Hoffman interactive failure criteria predict when a given set of 

stresses will produce failure, but they do not predict the mode of failure. Hashin 

developed a set of interactive failure criteria in which distinct failure modes are modeled. 

He applied his failure criteria to fiber composite materials. Since most fiber composites 

are transversely isotropic (e.g., wood), Hashin defined a general failure criterion in 

terms of the stress invariants of a transversely isotropic material (27). The five stress 

invariants (I1, I2, I3, I4, and I5) are: 

1 = σ 11 I (159)

2 = σ 22 +σ 33 I (160)

22 33

2 I3 = σ 33 −σ σ (161)

2

13

2 I4 = σ12 +σ (162)

2

33 12

2 I5 = 2σ12σ 23σ13 −σ 22σ13 −σ σ (163)

Hashin’s three-dimensional failure criterion is a quadratic stress polynomial of the 

general form: 

1 12 1 2 3 3 4 4

2

2 2 2 2

2

1 1 1 1 A I + B I + A I + B I + C I I + A I + A I ≥ (164)

The I5 invariant does not appear in the criterion because only linear and quadratic terms 

are retained in this polynomial. 

To identify distinct failure modes, Hashin argued that failure is produced by the normal 

and shear stresses acting on the failure plane. For failure parallel to the grain, the failure 

plane is the 2-3 plane, acted on by stresses σ11, σ12, and σ13. The perpendicular 

stresses (σ22, σ23, and σ33) do not contribute to parallel failure. The implicit assumption 

here is that the perpendicular stresses do not impede compression bucking; thus, an 

interaction mechanism is not required (the term C12I1I2 is neglected). 
110

In perpendicular-to-the-grain failure, failure occurs in any plane with axes parallel and 

perpendicular to the grain. The failure plane is acted on by stresses σ22, σ33, σ23, σ12, 

and σ13. The implicit assumption here is that the stress parallel to the grain (σ11) does 

not contribute to perpendicular failure because this stress is carried almost entirely by 

the fibers. 

By applying these assumptions to the general criterion in equation 164, Hashin 

developed specific yield criteria for the parallel and perpendicular modes: 

Parallel Mode 

1 ( )

2

||

2

13

2

2 12

11 11 = +

+ +

S

Af B f

σ σ

σ σ (165)

Perpendicular Mode 

1 ( ) ( ) ( ) ( ) 2

||

2

13

2

12

2

22 33

2

2 23

22 33 22 33 = +

+ − + + + +

S⊥ S

Am Bm

σ σ σ σ σ

σ σ σ σ (166)

Failure mechanisms are different for tensile and compressive modes, so Hashin further 

divided each criterion into tensile and compressive modes. 

Tension Parallel: Hashin assumed that tensile and shear stresses are mutually 

weakening; therefore, both contribute to tensile failure. Data from direct pull and biaxial 

tests are needed to solve for both coefficients (Af and Bf). If data from a direct pull test 

are the only data available, then one can solve for either Af or Bf. The failure surface 

remains smooth and elliptical if one solves for Bf and neglects Af. Failure occurs when 

the following equation is satisfied: 

1 ( 0) ( )

11 2

||

2

13

2

12

2

2

11 ≥ >

+

+ σ

σ σ σ

XT S

(167)

This criterion is the same as our extension of the Yamada-Sun criterion. 

Compression Parallel: If compressive and shear stresses are assumed to be mutually 

weakening, then one can develop a compression parallel criterion similar to that for 

tension. However, Hashin argued that there is no physically reasonable method for 

including the effect of shear stress, at least for composites. Therefore, he represents 

parallel compressive failure in simple maximum stress form: 
111

( 0) σ11 ≥ XC σ11 < (168)

Tension Perpendicular: Data from direct pull and biaxial tests are needed to solve for 

both coefficients (Am and Bm). If data from a direct pull test are the only data available, 

then one can solve for either Am or Bm. The failure surface remains smooth and elliptical 

if one solves for Bm and neglects Am. Failure occurs when the following equation is 

satisfied: 

1 ( 0) ( ) ( ) ( )

22 33 2

||

2

13

2

12

2

22 33

2

23

2

2

22 33 ≥ + >

+

+ − +

+

⊥

σ σ

σ σ σ σ σ σ σ

Y S S T

(169)

Compression Perpendicular: Hashin argued that the biaxial compressive strength 

(Ycc) of composites is much greater than the uniaxial compressive strength (Yc). 

Therefore, he solved for both coefficients (Am and Bm) and retained only first-order 

terms in Yc /Ycc. Failure occurs when the following equation is satisfied: 

1 ( ) ( ) ( ) 1

4 2

( )

2

||

2

13

2

12

2

22 33

2

23

2

2

22 33

2

2

||

2

22 33 ≥

+

+ − +

+

⎥

⎥

⎦

⎤

⎢

⎢

⎣

⎡

− ⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛

+

+

S⊥ Y S⊥ S

Y

S c

C σ σ σ σ σ σ σ σ σ

( 0) σ 22 +σ 33 <

(170)

The Modified Hashin criteria predict four independent modes of failure: tensile and 

compressive failure parallel to the grain, and tensile and compressive failure 

perpendicular to the grain. Although Hashin applied these criteria to fiber composites, 

these criteria were also evaluated for wood. A plane stress version is currently 

implemented in model 22 of LS-DYNA, along with an augmentation by Chang that takes 

into account nonlinear shear stress-strain behavior. 

Modified Hashin: In addition to the composite criteria proposed by Hashin, a simple 

modification of the Hashin criteria was evaluated. 

Tension and Compression Parallel: Here, it is assumed that shear stress weakens 

wood in compression as well as in tension. In this case, the tensile and compressive 

yield criteria have the following form: 

⎩

⎨

⎧

<

>

≥ = ++

0

0

1 ( )

C 11

T 11

2

||

2

13

2

12

2

2

11

σ

σ σ σ σ

X for

X for X

X S (171)

Tension and Compression Perpendicular: It is not clear whether the biaxial 

compressive strength of wood is greater or lesser than the uniaxial compressive 

strength. The only available data for southern yellow pine are uniaxial stress data, so no 
112

assumptions are made regarding the relative strengths in biaxial and uniaxial 

compression. In this case, the tensile and compressive yield criteria have the following 

form: 

⎩

⎨

⎧

+ <

+ >

≥ = ++ − +

+

⊥ ⊥ 0

0

1 ( ) ( ) ( )

C 22 33

T 22 33

2

2

13

2

12

2

22 33

2

23

2

2

22 33

σ σ

σ σ σ σ σ σ σ σ σ

Y for

Y for X

Y S S

(172)

The Modified Hashin criteria predict four independent modes of failure: tensile and 

compressive failure parallel to the grain, and tensile and compressive failure 

perpendicular to the grain. 

Hankinson: Hankinson developed an empirical formula that is frequently applied to offaxis wood tests in two dimensions.(16,34) Goodman and Bodig extended the formula to 

three dimensions.(7) Their three-dimensional formula predicts the ultimate compressive 

strength (σult) of wood relative to the both the grain (θ) and ring angles (φ). The threedimensional formula has the following form: 

θ φ

σ

φ

φ

2 2 sin F cos

F

ult X

+ = (173)

where Fφ varies sinusoidally with the ring angle between the relative strength in the 

tangential direction (FT = Y/X) and the relative strength in the radial direction (FR = Z/X). 

Fφ has the following form: 

2

( ) ( sin 2 ) 90

( R T ) R T

φ T

F F K F F F F +

+ − − = + φ θ (174)

Here, K is an empirical constant that is typically 0.4 for softwoods. The last term on the 

right of 172 is a sinusoidal correction to the straight line interpolation between the 

tangential (Y) and radial (Z) strengths. This correction is illustrated in figure 45. It can be 

compared with the data shown later in figure 48. 

Comparisons of Hankinson’s three-dimensional formula with Douglas fir test data were 

previously shown in figure 2. Hankinson’s formula predicts large reductions in strength 

when the load is inclined at a small angle to the grain or ring, in agreement with test 

data. 
113

Hankinson model 

Figure 45. Compressive strength variation of clear wood is readily modeled by a 

sinusoidal correction in the R-T plane. Source: Krieger Publishing Company.(16) 

For a transversely isotropic material with no variations in strength with the ring angle 

(FT = FR = Fφ), the three-dimensional formula reduces to its two-dimensional form: 

θ θ

σ 2 2 X sin Y cos

XY

ult

+ = (175)

Although Hankinson’s formula has been shown to provide good fits to off-axis test data, 

it is not a general-purpose formulation that can be applied to other types of tests. Thus, 

it is not suitable for use in finite element codes. However, it is reported here so that it 

can be compared with other criteria to help in their evaluation. 

115

APPENDIX D. GRAPHICAL COMPARISON OF CANDIDATE 

FAILURE CRITERIA 

D.1 COMPARISON RESULTS 

Each failure criterion was evaluated by comparing it with available test data and with 

other criteria. Comparisons include: 

• Off-axis strength predictions. 

• Parallel-to-the-grain strength predictions. 

• Perpendicular-to-the-grain strength predictions. 

Details of the comparisons are given in section D.2. The results of the comparisons are 

briefly given here: 

• The failure stresses for all candidate criteria agree for states of uniaxial 

stress in the principal material directions. The criteria disagree on what 

constitutes failure for biaxial and triaxial stress states. 

• Nearly sufficient uniaxial data are available for evaluating the candidate 

criteria. The exception is the lack of a perpendicular shear strength 

measurement for southern yellow pine. Insufficient biaxial and triaxial data 

are available to thoroughly evaluate each criterion. Biaxial and triaxial data 

are limited to off-axis test data for Douglas fir. 

• The Modified Hashin and Tsai-Wu criteria fit off-axis test data the best. All 

candidate failure criteria were compared with off-axis test data in two 

planes: the L-T and T-R planes. In the L-T plane, all criteria except 

Maximum Stress and Hashin accurately predict the trend observed in offaxis test data. In the R-T plane, no criterion accurately predicts the trend 

observed in off-axis test data. However, the transversely isotropic criteria 

(Hashin, Modified Hashin, and Tsai-Wu) provide the best fit. These 

comparisons suggest that the Modified Hashin and Tsai-Wu criteria are 

the best candidates. 

• The Modified Hashin criterion is more practical than the Tsai-Wu criterion 

because it identifies four distinct failure modes. The Hashin, Yamada-Sun, 

Norris, and Maximum Stress criteria also identify the mode of failure. On the 

other hand, the Tsai-Wu and Hoffman criteria predict the onset of failure, but 

not the mode of failure. This makes it difficult to model postpeak softening, 

because the extent of softening depends on the failure mode (tensile, 

compression, or shear) and failure direction (parallel or perpendicular). 
116

• Simultaneous and/or sequential failure in perpendicular and parallel 

modes is modeled by those criteria that predict distinct failure modes 

(Modified Hashin, Hashin, Yamada-Sun, Norris, and Maximum Stress). 

This is because these criteria are formulated with multiple equations and 

more than one failure equation can be satisfied at one time. Those criteria 

formulated with a single equation (Tsai-Wu and Hoffman) cannot make 

such predictions. 

• Most criteria produce smooth failure surfaces in stress or stress invariant 

space. Having a smooth failure surface is important when applying the 

plasticity algorithm because nonsmooth surfaces (with corners) require 

additional coding and add computational intensity. The only nonsmooth 

surface that was evaluated is the Hashin criterion, and it occurs in the 

perpendicular modes. Most criteria, such as the Modified Hashin, produce 

segmented surfaces when cross-plots are made involving both the parallel 

and perpendicular modes; however, this is not a computational issue 

because plasticity will be treated separately for each mode. The Modified 

Hashin criterion produces a smooth surface for the parallel modes and a 

separate smooth surface for the perpendicular modes. 

• Most criteria model little or no interaction between the parallel and 

perpendicular modes. This means that parallel strength is not significantly 

affected by perpendicular confinement. Five of the seven criteria predict 

little or no increase in parallel compressive strength with perpendicular 

confinement. The Modified Hashin criterion is one of the five criteria. The 

two exceptions are the Hoffman and Tsai-Wu criteria. Four of the seven 

criteria predict little or no increase in parallel tensile strength with 

perpendicular confinement. The Modified Hashin criterion is one of the 

four criteria. The three exceptions are the Hoffman, Tsai-Wu, and Norris 

criteria. Interaction terms could readily be added to the Modified Hashin 

criterion if test data become available that suggest that an interaction 

mechanism is appropriate. 

• The transversely isotropic criteria conservatively model parallel failure 

under parallel shear stress. These criteria (Modified Hashin, Hashin, 

Yamada-Sun, and Tsai-Wu) predict shear failure at lower shear stresses 

than the orthotropic criteria. The lower shear strength results from 

simultaneous application of two shear stress components. 

• The strength predictions of the transversely isotropic criteria (Modified 

Hashin, Hashin, Yamada-Sun, and Tsai-Wu) are realistic under 

transformation of stress in the isotropic plane. These criteria recognize 

σT = –σR as a state of pure shear stress; the orthotropic criteria do not. 

• The transversely isotropic criteria are more flexible in modeling failure and 

yielding in the perpendicular modes than the orthotropic criteria. This is 
117

because the transversely isotropic predictions are sensitive to the value of 

the perpendicular shear strength selected by the analyst. 

D.2 GRAPHICAL COMPARISONS 

The criteria are plotted using the transversely isotropic strength values previously 

reported in table 4 for southern yellow pine at 12-percent moisture content. 

In some plots, the various failure criteria are compared with Hankinson’s formula. 

Hankinson’s formula is an empirical equation that is frequently fit to off-axis 

compression tests of wood in two dimensions. Goodman and Bodig extended the 

formula to three dimensions.(7) Although Hankinson’s formula fits off-axis test data fairly 

well, it is not a general-purpose formulation that can be applied to other types of tests. 

This is because it is explicitly formulated in terms of the off-axis grain and the ring 

angles. Thus, it is not suitable for use in finite element codes. However, it is reported 

here so that it can be compared with the various criteria to help in their evaluation. 

D.3 OFF-AXIS STRENGTH COMPARISONS 

Off-axis strength is characterized by performing uniaxial tests with the symmetry axis 

(L-T-R axis) oriented at-angle to the loading axis, as schematically shown in figure 46. 

The measured strength depends on two off-axis angles―the angle θ between the grain 

and loading axis, and the angle φ between the rings and loading axis. A grain angle of 

θ = 0 degrees means that the load is being applied in the longitudinal direction. A grain 

angle of θ = 90 degrees means that the load is being applied in the perpendicular 

direction. If φ = 0 degrees when θ = 90 degrees, then the load is being applied in the 

tangential direction. If φ = 90 degrees when θ = 90 degrees, then the load is being 

applied in the radial direction. 

The various failure criteria are compared here with off-axis test data (i.e., Hankinson’s 

formula). Two sets of comparisons are made. In the first, the off-axis strength in the L-T 

plane (parallel and perpendicular to the grain) is examined. In the second, the off-axis 

strength in the T-R plane (perpendicular to the grain only) is examined. 
118

Figure 46. Geometry of an off-axis test specimen. 

Source: Krieger Publishing Company.(16) 

D.3.1 Off-Axis Strength in the L-T Plane 

Data: Off-axis compressive strength data for Douglas fir were previously shown in the 

three-dimensional plot of figure 2.(7) Also shown is the three-dimensional Hankinson 

formula. The comparison is good, although Hankinson’s formula tends to underestimate 

the measured strength, particularly for grain angles between 15 and 45 degrees. 

Hankinson’s formula also tends to underestimate the measured strength of Engelmann 

spruce, oak, and aspen (not shown), as reported by Goodman and Bodig.(7) 

Although Hankinson’s formula is an empirical fit to wood data, the good fit does not 

validate the formula. This is because longitudinal, tangential, radial, and shear stresses 

all act on the wood in the material coordinate system; however, only longitudinal, 

tangential, and radial stresses are assumed to contribute to compressive failure in 

Hankinson’s formula. The relative effects of the normal and shear stresses are not 

established. 

Strength Comparisons: The measured off-axis strengths are not tabulated in the 

paper by Goodman and Bodig and are difficult to extract from figure 2.(7) In addition, no 

off-axis test data are available from FPL for southern yellow pine. Therefore in figure 

47(a), the failure criteria are compared with Hankinson’s formula rather than test data 

(keep in mind that Hankinson’s formula tends to underestimate the measured strength). 
119

The comparisons are for the ultimate compressive strength of southern yellow pine 

calculated as a function of grain angle for φ = 0 degrees. 

All but two criteria, the Maximum Stress and Hashin criteria, are in good agreement with 

Hankinson’s formula. In fact, most criteria predict a slightly greater strength than 

Hankinson’s formula, a trend consistent with the measured data. On the other hand, the 

Maximum Stress criterion significantly overestimates the strength of the wood compared 

with Hankinson’s formula. The Maximum Stress criterion predicts an increase in off-axis 

strength with increasing grain angle between 0 and 23 degrees. This is opposite the 

trend predicted by Hankinson’s formula and observed in off-axis tests of various wood 

species. In addition, the Hashin criterion overestimates the strength of wood compared 

with Hankinson’s formula for small grain angles and underestimates the strength at 

moderate grain angles. These comparisons suggest that the Maximum Stress and 

Hashin criteria are not good candidates for modeling the off-axis strength of wood. 

In addition to compressive strength, tensile strength is also plotted as a function of grain 

angle in figure 47(b). The primary reason these tensile comparisons were made was to 

check the Tsai-Wu criterion. One coefficient of the Tsai-Wu criterion must be fit to 

biaxial strength data. The criterion was fit to the off-axis compressive strength predicted 

by Hankinson’s formula at 45 degrees. Even though the fit was made in compression, 

the Tsai-Wu criterion is in agreement with Hankinson’s formula and the other interactive 

criteria in tension. Although Goodman and Bodig did not report comparisons of 

Hankinson’s formula with wood test data in tension, various authors suggest that 

Hankinson’s formula is a reasonable fit in tension as well. 
120

(a) Compressive strength 

(b) Tensile strength 

Figure 47. Most of the interactive failure criteria are in agreement 

with Hankinson’s formula for the off-axis strength 

of southern yellow pine in the L-T plane. 
121

Failure Mode Comparisons: Although most of the interactive criteria accurately predict 

the off-axis compressive strength, it is not known if the criteria predict the correct failure 

mode. This is because the measured failure mode was not reported by Goodman and 

Bodig.(7) We expect the failure to be in the parallel modes at low grain angles and in the 

perpendicular modes at high grain angles. Five of the failure criteria (Hashin, Modified 

Hashin, Yamada-Sun, Norris, and Maximum Stress) predict distinct failure modes. 

Considering the Modified Hashin criterion, for grain angles less than about 24 degrees, 

the predicted mode is compression failure parallel to the grain. For grain angles greater 

than about 24 degrees, the predicted mode is compression failure perpendicular to the 

grain. This change in mode is evident by the discontinuity in slope in figure 47(a). 

Failure modes are readily identified because Hashin formulated separate criteria 

(equations) for the parallel and perpendicular modes. Two of the criteria, the Tsai-Wu 

and Hoffman criteria, do not distinguish the mode of failure because they are formulated 

with a single equation. 

D.3.2 Off-Axis Strength in R-T Plane 

Data: Off-axis test data are also useful for evaluating the assumption of transverse 

isotropy. Transversely isotropic strength values were previously reported in table 4 for 

pine and in table 5 for fir. This means that no distinction is made between the strengths 

in the tangential and radial directions, so only one measurement is made and labeled as 

perpendicular. The assumption of transverse isotropy is assessed by examining the offaxis strength data measured by Goodman and Bodig in the R-T (isotropic) plane. These 

data are shown in the two-dimensional plot of figure 48. Goodman and Bodig measured 

the ultimate strength at various ring angles, holding the grain angle constant at 

θ = 90 degrees. These data are for four different wood species, one of which is Douglas 

fir. 

If the wood species in figure 48 were transversely isotropic, then the tangential strength 

would be equal to the radial strength and the data would form a straight line between 0 

and 90 degrees. The data at 0 degrees is the tangential strength. The data at 90 

degrees is the radial strength. For Douglas fir, the radial strength is about 85 percent of 

the tangential strength, which is reasonably close. However, all of the data follow a 

similar pattern―the off-axis strength measured at 45 degrees is less than that 

measured at 0 or 90 degrees. For Douglas fir, the strength at 45 degrees is about 

60 percent of the tangential strength. 
122

Figure 48. Effect of ring angle variation at 90-degree grain angle on 

the relative compression strength of four wood species. 

Source: Society of Wood Science and Technology.(7) 

Although strength varies with the ring angle for these wood species, the variation is not 

great when compared with the compressive strength in the longitudinal direction. For 

Douglas fir, the tangential strength at 0 degrees is about 10 percent of the longitudinal 

strength. The off-axis strength at 45 degrees is about 6 percent of the longitudinal 

strength. Therefore, the variation in perpendicular strength is about 4 percent of the 

parallel strength. The assumption of transverse isotropy in strength is probably 

reasonable, especially if wood posts fail catastrophically in the parallel modes rather 

than in the perpendicular modes. 

Strength Comparisons: Hankinson’s formula in the R-T plane is an excellent fit to 

wood data, as previously shown in figure 48. All failure criteria are compared with 

Hankinson’s formula in the R-T plane in figure 49. None of the criteria is in agreement 

with Hankinson’s formula, except in uniaxial compression at 0 and 90 degrees. At 

45 degrees, the ultimate strength predicted by the various failure criteria is greater than 

that predicted by Hankinson’s formula. Hankinson’s formula suggests a 40-percent 

reduction in strength at 45 degrees. The Tsai-Wu, Hashin, and Modified Hashin criteria 

predict no reduction in strength at 45 degrees and no variation with ring angle. This is 

because they were derived from the invariants of a transversely isotropic material and 

the strength of a transversely isotropic material does not vary with ring angle. All other 

orthotropic criteria predict an increase in ultimate strength at 45 degrees (in poor 

Ring Angle (degrees)

Relative Compression Strength 

Legend

 Measured Predicted 

Englemann Spruce o _____ 

Douglas Fir  __ __ 

Oak Δ _ _ _ _ 

Aspen • . 

o



Δ

•
123

Figure 49. Failure criteria comparison for perpendicular modes 

as a function of the ring angle. 

agreement with Hankinson’s formula and the trend observed in the measured data 

(figure 48)). It is interesting to note that the orthotropic failure criteria do not predict 

transversely isotropic behavior even though transversely isotropic strength values are 

used. These comparisons suggest that the best criteria for modeling biaxial 

compressive strength perpendicular to the grain are the transversely isotropic criteria 

(Tsai-Wu, Hashin, and Modified Hashin). 

Transformations: Failure criteria are applied to stresses in the principal material 

directions. The off-axis strength, σult, is transformed into stresses in principal material 

directions: 

⎥

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎢

⎣

⎡

−

=

⎥

⎥

⎥

⎦

⎤

⎢

⎢

⎢

⎣

⎡

mn

n

m

ult

2

2

TR

R

T

σ

σ

σ

σ

(176)

where: 

m = cosφ

n = sinφ

φ = ring angle between the tangential direction and loading axis 

(F) 
124

The state of stress at 45 degrees is biaxial, with σT = σult / 2, σR = σult / 2, and 

σTR = σult / 2. Referring back to the Douglas fir data in figure 48, one sees that the equal 

biaxial compressive stress perpendicular to the grain is 60 percent of the tangential 

strength, which, in turn, is half of the applied stress. Therefore, the biaxial compressive 

strength of Douglas fir is 30 percent of the uniaxial (tangential) compressive strength, at 

least in the presence of shear stress. On the other hand, the transversely isotropic 

failure criteria predict a biaxial compressive strength that is 50 percent of the uniaxial 

compressive strength, also in the presence of shear stress. 

To accurately measure the biaxial compressive strength, one would need data from 

biaxial tests performed with and without the application of shear stress. If one argues 

that shear stress is mutually weakening, as it is for composites, then one would expect 

the biaxial compressive strength measured with shear stress to be less than that 

measured without shear stress. Some failure criteria predict such a trend, as discussed 

in section C.3. 

D.4 PARALLEL-TO-THE-GRAIN STRENGTH COMPARISONS 

Wood posts are observed to fail in tension or shear parallel to the grain. Therefore, it is 

particularly important to accurately model the critical combinations of stresses that 

produce failure in the parallel modes. Although ultimate stress versus grain angle plots 

from off-axis tests include assessment of parallel failure, they reveal few differences 

among the various criteria. A more exacting assessment is attained with biaxial stress 

plots. The biaxial stress plots discussed in subsequent paragraphs indicate that 

significant differences in strength are predicted by the various criteria. These differences 

are revealed in three sets of stress plots: 

• Longitudinal stress versus perpendicular stress. 

• Longitudinal stress versus parallel shear stress. 

• Longitudinal stress versus combined stress (perpendicular plus parallel shear 

stress). 

D.4.1 Biaxial Comparisons of Longitudinal Versus Confining Stress 

Strength Comparisons: The combinations of longitudinal and confining stresses that 

satisfy the various failure criteria are compared in the biaxial strength plot of figure 50 

for southern yellow pine. By confining stress, we mean the sum of the tangential and 

radial stresses. This sum is one of the invariants of a transversely isotropic material. 

The sum was obtained for the specific case of equal tangential and radial stresses. This 

figure is plotted with the stresses positive in tension. 

Stress states that lie on the vertical axes (σT + σR = 0) are uniaxial stress states. When 

the confining stress is zero, each curve intersects the vertical axis twice―once in 

tension (positive) and once in compression (negative). These intersections are the 
125

longitudinal strengths in uniaxial tension and compression. Stress states that lie on the 

horizontal axes (σL = 0) are biaxial stress states. When the longitudinal stress is zero, 

each curve (except Hoffman) intersects the horizontal axis twice―once in tension and 

once in compression. These intersections are twice the perpendicular strengths in equal 

biaxial tension and compression. Stress states that do not lie on the vertical or 

horizontal axes are triaxial stress states. All criteria are in agreement for states of 

uniaxial stress. Uniaxial strength measurements are available for fitting the criteria. The 

criteria disagree on what constitutes failure or yielding for states of biaxial and triaxial 

stress. Biaxial and triaxial strength measurements are not currently available for fitting 

the criteria. Biaxial and triaxial test data are needed to validate one or more of the 

proposed criteria. 

Figure 50. Predicted effect of perpendicular confinement and extension on the 

longitudinal strength of southern yellow pine in tension and compression. 

Failure Mode Comparisons: All criteria plotted in figure 50 are closed, but not 

necessarily smooth surfaces.6

 The failure surfaces that are smooth are the Tsai-Wu and 

Hoffman criteria. This is because failure is modeled by a single formula or equation. The 

drawback of these criteria is that only the onset of failure is predicted, not the mode of 

failure. This means that when a particular combination of stresses indicates that failure 

or yielding has occurred, the criteria do not indicate the type of failure (tensile, 

compressive, or shear) or its direction (parallel or perpendicular). It is important to know 

the type and direction of failure so that the correct softening response (brittle or ductile) 

and fracture energy are modeled, as discussed in section 1.7. 

 

6 Although the Hoffman criterion appears to model infinite biaxial compressive strength perpendicular to 

the grain, it actually models an extremely large biaxial compressive strength (off the scale of the figure). 

(C) 

(F) 

(D) 
126

The nonsmooth criteria indicate both the type and direction of failure. This is because 

failure is modeled by more than one equation. For example, the Modified Hashin 

criterion models four modes of failure because four different failure equations are 

satisfied on different portions of the failure surface, which is shaped like a square. 

These modes are compressive yielding perpendicular to the grain (left side of square), 

tension failure parallel to the grain (top side of square), tension failure perpendicular to 

the grain (right side of square), and compressive yielding parallel to the grain (bottom 

side of square). 

All criteria, except Hoffman and Tsai-Wu, predict no increase, or a small increase in 

compressive strength with compressive confinement (refer to the lower left-hand 

quadrant of figure 50). The Hoffman criterion predicts a large increase in compressive 

strength with confinement. The Tsai-Wu criterion predicts a reduction in strength with 

confinement. The only strength data available for comparison that include confinement 

are the off-axis test data; however, such data also include contributions from shear 

stress. The effect of shear stress on longitudinal strength is discussed in subsequent 

paragraphs. 

D.4.2 Biaxial Comparisons of Longitudinal Versus Combined Stress 

The combinations of longitudinal and tangential stresses that satisfy the various failure 

criteria are compared in the biaxial strength plot of figure 51(a) for southern yellow pine. 

These curves include a contribution from the parallel shear stress and were calculated 

for the specific case of σLT = σT, with σLR = 0 and σR = 0. Also included in this figure is 

one point from Hankinson’s formula that represents off-axis test data measured at 

θ = 45 degrees. Almost all criteria, expect the Maximum Stress criterion, are in good 

agreement with the Hankinson point, as expected from previous comparisons with offaxis test data. 

The failure criteria differ most in the upper left-hand quadrant of figure 50. This is the 

quadrant that predicts the longitudinal tensile strength as a function of perpendicular 

confinement and parallel shear. The two criteria that differ most are the Tsai-Wu and 

Norris criteria. The Tsai-Wu criterion predicts an increase in longitudinal tensile strength 

with compressive confinement and shear, while the Norris criterion predicts a large 

decrease in longitudinal tensile strength with compressive confinement and shear. The 

remaining criteria predict little or no reduction in strength with combined confinement 

and shear. The small reduction in strength predicted by the Hashin, Modified Hashin, 

and Yamada-Sun models is a result of the application of shear. This is evident by 

comparing the failure criteria in figure 51(a), which were calculated with shear, with the 

failure criteria in figure 51(b), which were calculated without shear. No data are 

available in this quadrant for evaluating each model. More clear wood data are needed 

to understand the relative contributions of shear stress and perpendicular confining 

stress on the longitudinal strength in tension. 
127

D.4.3 Biaxial Comparisons of Longitudinal Versus Shear Stress 

All clear wood strength data available for southern yellow pine or Douglas fir are a 

measurement of either σLT or σLR, or, more generally, of parallel shear strength. The 

combined effect of applying two parallel shear stress components simultaneously has 

not been measured. Wood posts have been observed to fail in a shear mode, so it is 

important to establish the correct shear strength, particularly if failure is caused by 

simultaneous application of two shear components. The effects of unilateral and 

bilateral shear stresses on longitudinal strength are examined here. 

The combinations of longitudinal and shear stresses that satisfy the various failure 

criteria are compared in the biaxial strength plot of figure 52 for southern yellow pine. By 

shear stress, we mean the sum of the squares of the L-T and L-R parallel shear stress 

terms ( ) 2 2 σ LT +σ LR . This sum is one of the invariants of a transversely isotropic 

material. The square root is taken to retain units of stress. Two plots are shown. In one 

plot, the sum was obtained for the specific case of σLT = σLR. In the other plot, the sum 

was obtained for the specific case of σLR = 0. 
128

(a) Calculated with parallel shear stress equal to tangential stress 

(b) Calculated without parallel shear stress 

Figure 51. Predicted effect of parallel shear and tangential stresses on the 

longitudinal strength of southern yellow pine in tension and compression. 

(A) 

(B) 

(C) 

(F) 

(E) 

(D) 

(G) 

(A) 

(B) 

(C) 

(D) 

(E) 

(F) 

(G) 

(H) 
129

(a) Calculated with parallel shear stress equal to tangential stress 

(b) Calculated without parallel shear stress 

Figure 52. Predicted effect of parallel shear invariant on the longitudinal strength 

of southern yellow pine in tension and compression. 

(B) 

(F) 

(G) 

(A) 

(D) 

(C) 

(E) 

(E) 

(F) 

(G) 

(A) 

(B) 

(D) 

(C) 
130

In figure 52(a), the transversely isotropic and orthotropic criteria predict different shear 

strengths. The orthotropic criteria (Maximum Stress and Norris) predict a combined 

shear strength that is 40 percent greater than that predicted by the transversely 

isotropic criteria. 

The transversely isotropic criteria predict a value of 16.8 MPa, which is equal to the 

shear strength measured parallel to the grain for either the L-T or L-R shear stress 

component. The individual values of σLT and σLR do not affect the parallel shear strength 

of a transversely isotropic material, only the sum of their squares. This is because this 

sum is an invariant of a transversely isotropic material and the parallel shear stress 

components are not included in the other invariants used in the transversely isotropic 

models. However, this sum is not an invariant if the criterion is orthotropic, so the 

specific values of each stress component affect the predicted strength. For the specific 

case of equal shear stress in the L-T and L-R planes, this sum is equal to 23.7 MPa. 

These comparisons demonstrate that the combined shear strength predicted by 

transversely isotropic criteria is lower than that predicted by orthotropic criteria. Use of 

transversely isotropic criteria in guardrail post calculations could result in a lower shear 

strength than that predicted by orthotropic criteria. Planned guardrail post calculations 

will examine the value of each shear stress component and check the sensitivity of the 

results to shear strength. 

D.5 PERPENDICULAR-TO-THE-GRAIN STRENGTH COMPARISONS 

Wood guardrail posts fail catastrophically in the parallel modes (tension and shear). 

However, failure and yielding in the perpendicular modes could precede parallel failure 

without causing catastrophic failure of the guardrail post. Failure in the perpendicular 

modes occurs at much lower strengths than failure in the parallel modes. It is important 

to accurately model the perpendicular failure criteria in order to limit the perpendicular 

stresses. 

This section demonstrates that significant differences exist in the perpendicular strength 

predicted by the various failure criteria. In addition, the transversely isotropic criteria 

have an advantage over the orthotropic criteria in that they are more flexible in fitting 

data. This is demonstrated by comparing the candidate failure criteria with each other 

and with Hankinson’s formula. Two sets of biaxial stress plots are evaluated: 

• Radial versus tangential stress. 

• Shear versus tangential stress. 

The biaxial stress plots provide a more exacting assessment of the perpendicular failure 

criteria than the off-axis strength predictions previously discussed. 
131

D.5.1 Biaxial Comparisons of Radial Versus Tangential Stress 

The biaxial strength of southern yellow pine calculated without the application of shear 

stress is examined here. The combinations of radial and tangential stresses that satisfy 

the various failure criteria are compared in the biaxial strength plot of figure 53 for 

southern yellow pine. Stress states that lie on the horizontal (σR = 0) and vertical 

(σT = 0) axes are uniaxial stress states. Stress states that do not lie on the horizontal 

and vertical axes are biaxial stress states. When the radial stress is zero, each curve 

intersects the horizontal axis twice―once in tension (positive) and once in compression 

(negative). These intersections are the tangential strengths in uniaxial tension and 

compression. When the tangential stress is zero, each curve intersects the vertical axis 

twice―once in tension and once in compression. These intersections are the radial 

strengths in uniaxial tension and compression. The failure stresses for all criteria are in 

agreement for states of uniaxial stress. The criteria disagree on what constitutes failure 

for states of biaxial stress. Here, three biaxial states are examined: equal biaxial 

compression, equal biaxial tension, and pure shear. 

Figure 53. Predicted strength of southern yellow pine perpendicular 

to the grain (no perpendicular shear stress applied). 

Equal Biaxial Compression: States of equal biaxial compression and tension (σT = σR) 

are visualized as a diagonal line extending from the lower left-hand to the upper righthand corners of the plot. There are three main clusters of curves in biaxial compression 
132

(lower left-hand quadrant). First, the Maximum Stress, Norris, and Yamada-Sun criteria 

predict a biaxial compressive strength equal to the uniaxial compressive strength. 

Second, the Hoffman criterion predicts a biaxial compressive strength that is much 

greater (off the scale of the plot) than the uniaxial compressive strength. Third, the 

Hashin, Modified Hashin, and Tsai-Wu criteria predict biaxial compressive strengths that 

are less than the uniaxial compressive strength. No strength data are available in biaxial 

compression for comparison with the failure criteria that do not include contributions 

from shear stress. 

Equal Biaxial Tension: There are two main clusters of curves in biaxial tension (upper 

right-hand quadrant). First, the Maximum Stress, Norris, and Yamada-Sun criteria 

predict a biaxial tensile strength that is equal to the uniaxial tensile strength. In addition, 

the Hoffman criterion predicts a biaxial tensile strength that is approximately equal to 

the uniaxial tensile strength. The biaxial tensile strength predicted by the Hoffman 

criterion is sensitive to the strengths used to fit the model. Second, the Hashin, Modified 

Hashin, and Tsai-Wu criteria predict a biaxial tensile strength that is less than the 

uniaxial tensile strength. No data are available in biaxial tension for comparison with the 

failure criteria. 

Pure Shear: If wood is assumed to be transversely isotropic, then the R-T plane is the 

isotropic plane. Mohr’s circle indicates that equal normal stresses of opposite sign 

(σ22 = –σ33), calculated in the material principal directions, can be transformed into a 

state of pure shear stress. States of shear can be visualized as a diagonal line 

extending from the upper left-hand to the lower right-hand corners of the plot. There are 

three main clusters of curves in shear. First, the Maximum Stress and Yamada-Sun 

criteria predict a shear strength that is equal to the uniaxial tensile strength 

perpendicular to the grain. Second, the Hashin, Modified Hashin, and Tsai-Wu criteria 

predict a shear strength that is equal to the shear strength measured perpendicular to 

the grain (S⊥). A value for S⊥ is not reported in table 4; therefore, it is assumed that S⊥ is 

equal to 1.4 times the parallel shear strength (S||) and is thus greater than the uniaxial 

tensile strength. These transversely isotropic criteria explicitly include σT = –σR as a 

state of shear with strength S⊥. Third, the Norris and Hoffman criteria predict shear 

strengths that are less than the uniaxial tensile strength and less than S⊥. These criteria 

do not explicitly recognize σT = –σR as a state of shear stress upon transformation. No 

strength data are available for σT = –σR for comparison with the failure criteria. 

Assessment: One advantage of the transversely isotropic criteria is that the shape of 

the failure surface is readily modified as a function of S⊥. Figure 54 displays the Hashin, 

Modified Hashin, and Tsai-Wu criteria for three different values of S⊥. Although S⊥ is 

hard to measure, its inclusion in the failure criteria is not only realistic, but it provides us 

with flexibility for modeling failure and yielding in the perpendicular modes. The 

orthotropic criteria, plotted in this plane, do not vary with S⊥ and, therefore, lack the 

flexibility of the transversely isotropic criteria. 
133

Figure 54. Shape of the failure surface is sensitive to perpendicular shear 

strength if the criteria are transversely isotropic. 

D.5.2 Biaxial Comparisons of Shear Versus Tangential Stress 

The purpose of this section is to determine if the biaxial strength is affected by the 

application of shear stress. The comparisons indicate that most criteria predict a 

reduction in biaxial strength with increasing shear stress. No test data are available to 

validate this trend. 

The combinations of shear and tangential stresses that satisfy the various failure criteria 

are compared in the biaxial strength plots of figure 55. One plot was calculated with 

S⊥ = 1.4 MPa and S|| = 23.7 MPa, while the other was calculated with S⊥ = 0.33 MPa 

and S|| = 5.58 MPa. This figure allows us to examine the effect of shear stress on 

perpendicular compressive strength. The only data available to the author for evaluating 

the failure criteria under states of biaxial stress with shear are off-axis test data, such as 

that previously shown in figure 2 for Douglas fir. As previously discussed, off-axis 

strength data include a contribution from shear stress. 

The stress combinations are plotted for the specific case of equal biaxial compression 

(σR = σT). Stress states that lie on the horizontal (σRT = 0) and vertical (σT = 0) axes are 

biaxial and pure shear stress states, respectively. When the shear stress is zero, each 

curve intersects the horizontal axis twice―once in tension (positive) and once in 

compression (negative). These intersections are the biaxial strengths in tension and 
134

compression. When the tangential stress is zero, each curve intersects the vertical axis 

twice―once in tension and once in compression. 

These intersections are the pure shear strengths. All of the criteria, except Maximum 

Stress and Hoffman, predict a reduction in biaxial strength with the application of shear 

stress. 

Also plotted is a solid black dot in the lower left-hand quadrant of each plot where the 

tangential stress is 30 percent of the uniaxial compressive strength and σTR = σR = σT. 

This is the biaxial stress state attained in off-axis tests at 45 degrees and predicted by 

Hankinson’s formula. Recall that Hankinson’s formula is in excellent agreement with test 

data in the T-R plane for a variety of woods. As expected, all failure criteria predict 

biaxial strengths that are greater than that predicted by Hankinson’s formula. 
135

(a) Calculated with S⊥ = 23.7 MPa 

(b) Calculated with S⊥ = 5.58 MPa 

Figure 55. Combinations of perpendicular and shear stresses 

that satisfy the failure criteria in the isotropic plane. 

(G) 

(A) 

(F) 

(D) 

(E) 

(C) 

(B) 

(H) 

(H) 

(H) 

(A) 

(F) 

(C) 

(B) 

(G) 

(D) 

(E) 

137

APPENDIX E. DERIVATION OF CONSISTENCY PARAMETER 

FOR PLASTICITY ALGORITHM 

Plasticity is modeled by enforcing separate consistency conditions for the parallel and 

perpendicular modes. The goal of each consistency condition is to partition the total 

strain increments into elastic and plastic components: 

p

ij

e

ij ij Δε = Δε + Δε (177)

The superscripts e and p indicate the elastic and plastic components, respectively. The 

total strain increments (Δεij) are calculated by the finite element code from the dynamic 

equations of motion and the time step. Once this partition is known, then the stress 

increments are updated from the elastic strain increments: 

( ) 1 p kl kl ijkl n

ij n σ ij =σ + C Δε − Δε + (178)

Here, n denotes the nth time step in the finite element analysis. 

E.1 CHECK FOR YIELDING 

The partition into elastic and plastic components requires two steps. The first step is to 

check for yielding. This is done by temporarily updating the stress components from the 

incremental strains by assuming that the entire strain increment is elastic: 

ijkl kl nij n σ ij =σ + C Δε * +1 (179)

These updated stresses are called the trial elastic stresses ( * σ ij ). The trial elastic stress 

invariants ( , , , ) *

4 *3 *2 *1I I I I are updated from the trial elastic stresses. The value of the 

yield function is evaluated from the trial elastic invariants and is denoted as f *. The new 

stress state is elastic if f * ≤ 0 and plastic if f * > 0. 

E.2 CALCULATE CONSISTENCY PARAMETER 

The second step is to enforce the consistency condition if f *> 0. Enforcement of the 

consistency condition requires an assumption about the direction of plastic flow. It is 

assumed that the plastic strain increments are normal to the yield surface: 
138

n ij

p

ij

f

σ

ε λ

∂

∂ Δ = Δ (180)

where Δλ is a proportionality constant known as the consistency parameter. 

This assumption is known as an associated flow rule, or normality condition. Use of a 

potential function other than the yield function in equation 180 results in a 

nonassociated flow rule. Recent studies reported by Pucik(35) suggest that rateindependent models with nonassociated flow lead to spurious (non-unique) dynamic 

solutions, so only associated flow is proposed for the present model. 

The plasticity algorithm calculates Δλ by enforcing the plastic consistency condition. 

This condition is expressed as: 

0 1 Δ = − = n+ n f f f (181)

where f is the yield surface function at time increments n to n + 1. The stress state at the 

beginning of the time step lies on the yield surface, thus, f

 n = 0. The stress state at the 

end of the time step is returned to the yield surface by the plasticity algorithm, thus, 

f n+1 = 0. Therefore, Δf = 0. 

The solution of the consistency condition in equation 181 determines Δλ, which, in turn, 

determines the partitioning of the total strain rate into elastic and plastic components. 

The stresses are updated from the elastic strain components. Separate Δλ solutions are 

proposed for the parallel and perpendicular modes. 

Parallel Modes 

The consistency condition is derived in terms of the invariants rather than the stresses, 

because the parallel failure criterion in equation 13 is formulated in terms of two 

transversely isotropic stress invariants. For the purposes of this derivation, the parallel 

failure criterion from equation 13 is defined as f||(I1,I4) ≥ 0, with: 

( , ) 1 2

||

4

2

2

1

|| 1 4 = + − S

I

X

I f I I (182)

A first-order Taylor series expansion of the consistency condition 0 ||

1

|| Δ || = − = n+ n f f f

from time increment n to n + 1 gives: 

0 4

4

||

1

1

|| Δ = ∂

∂

Δ +

∂

∂

I

I

f I

I

f

n n

(183)
139

Expansion of the stress invariant increments gives: 

||

||

* 1

1 1 λ

λ

Δ

∂

∂ Δ = Δ +

n

I I I (184)

||

||

* 4

4 4 λ

λ

Δ

∂

∂ Δ = Δ +

n

I I I (185)

where *

1 ΔI and *4 ΔI are the trial elastic increments calculated with Δλ|| = 0. Substitution 

of the updates from equations 184 and 185 into the consistency condition in equation 

183 results in the following expression for Δλ||: 

n n n n

n n

I

I

f I

I

f

I

I

f I

I

f

||

4

4

||

||

1

1

||

*

4

4

* ||

1

1

||

||

λ λ

λ

∂

∂

∂

∂

+

∂

∂

∂

∂

Δ

∂

∂

Δ − ∂

∂ −

Δ = (186)

The expression in the numerator is recognized as the first-order Taylor series expansion 

of *

|| f , where ( , ) *

4 * || 1 *|| f = f I I is the value of the failure criterion calculated from the trial 

elastic invariants. Therefore, the expression for Δλ|| reduces to: 

n n n n

I

I

f I

I

f

f

||

4

4

||

||

1

1

||

*

||

||

λ λ

λ

∂

∂

∂

∂

+

∂

∂

∂

∂

− Δ =

(187)

Perpendicular Modes 

The perpendicular failure criterion in equation 14 is derived in terms of the invariants 

rather than the stresses. The perpendicular failure criterion is defined as f⊥(I2,I3) = 0, 

with: 

( , ) 1 2

3

2

2

2

2 3 = + −

⊥

⊥

S

I

Y

I

f I I (188)
140

Expansion of the consistency condition Δf⊥ = 0 gives: 

3 0

3

2

2

Δ = ∂

∂ Δ +

∂

∂ ⊥ ⊥ I

I

f I

I

f

n n

(189)

The stress invariant updates are: 

⊥

⊥

Δ

∂

∂ Δ = Δ + λ

λ

n

I I I * 2

2 2 (190)

⊥

⊥

Δ

∂

∂ Δ = Δ + λ

λ

n

I

I I * 3

3 3 (191)

Substitution of these updates into the consistency condition gives the following express 

for Δλ⊥. 

n n n n

I

I

I f

I

f

f

⊥

⊥

⊥

⊥

⊥

⊥

∂

∂

∂

∂

+

∂

∂

∂

∂

− Δ =

λ λ

λ

3

3

2

2

*

(192)

E.3 UPDATE STRESSES 

The third step is to update the stresses. There are two options: (1) a purely elastic 

update and (2) an elastoplastic update. If the trial elastic stress state lies inside the 

failure surface (f *< 0), then Δλ = 0 and the stress state is purely elastic. In this case, the 

stress update from equations 178 and 179 is trivial: 

+1 * +1 = n

ij n σ ij σ (193)

If the trial elastic stress state lies outside the failure surface (f* > 0), then Δλ ≠ 0 and the 

stress state is elastoplastic. In this case, equations 178, 179, and 180 combine to give: 

kl n

ijkl

n

ij

n

ij

f C

σ

σ σ λ

∂

∂ = − Δ +1 * +1 (194)
141

APPENDIX F. DERIVATION OF LIMITING FUNCTION FOR 

HARDENING MODEL 

The functions that restrict the motion of each yield surface so that they cannot translate 

outside the ultimate surfaces are labeled G|| and G⊥ for the parallel and perpendicular 

modes, respectively. Each function is derived from the yield surface definition and 

hardening stress update. The limiting function for the parallel modes is derived here. 

The desired attributes of the limiting function are G|| = 1 at initial yield and G|| = 0 at 

ultimate yield. Hardening is modeled in compression, but not shear or tension, so the 

only stress component with hardening is σ11. The initial yield strength in compression is 

defined as σ11, and the ultimate strength in compression is defined as F

11 σ . The 

relationship between these strengths is: 

(1 ) 11 11 N||

F σ = σ − (195)

where 1 – N|| is the user-supplied reduction factor. For combined stress states, the 

ultimate yield strength from equation 13 is: 

2

||

4

11 1

S

I Xc

F σ = − (196)

For the case of uniaxial compressive stress, the ultimate yield strength is c

F σ11 = X . 

The longitudinal stress update with hardening is: 

σ 11 =σ 11 +α11 (197)

where α11 is the hardening stress (backstress). At ultimate yield, this relationship 

becomes: 

max

11 11 11 σ =σ +α F (198)

where max

11 α is the maximum backstress that can be attained. As previously defined, 

F

11 σ is the total stress with hardening (at ultimate yield) and σ11 is the stress without 

hardening (at initial yield). 
142

Substitution of equation 198 into equation 195 and rearranging gives: 

1 0

|| 11

max

11 − = F N σ

α

(199)

The above function has the desired attribute in that it equals zero when the stress state 

lies on the ultimate yield surface. Thus, one defines: 

F N

G

|| 11

11

|| 1

σ

α = − (200)

The value of the limiting function is G|| = 1 at initial yield because α11 = 0 at initial yield. 

The value of the limiting function is G|| = 0 at ultimate yield because max

11 11 α = α from 

equation 199. Thus, G|| limits the growth of the backstress as the ultimate yield surface 

is approached. 
143

APPENDIX G. SINGLE-ELEMENT INPUT FILE 

Figure 56. Single-element input file. 
144

Figure 57. First continuation of single-element input file. 
145

Figure 58. Second continuation of single-element input file. 

147 

REFERENCES 

1. LS-DYNA Keyword User’s Manual, Version 960, Volumes 1 and 2, Livermore 

Software Technology Corporation, Livermore, CA, 2001. 

2. Murray, Y.D., and J.D. Reid, et al., Evaluation of LS-DYNA Wood Material Model 

143, Report No. FHWA-HRT-04-096, Federal Highway Administration, 2004. 

3. Lewis, B.A., Manual for LS-DYNA Soil Material Model 147, Report No. FHWAHRT-04-095, Federal Highway Administration, 2004. 

4. Reid, J.D., B.A. Coon, B.A. Lewis, S.H. Sutherland, and Y.D. Murray, Evaluation 

of LS-DYNA Soil Material Model 147, Report No. FHWA-HRT-04-094, Federal 

Highway Administration, 2004. 

5. LS-DYNA User’s Manual, Livermore Software Technology Corporation, 

Livermore, CA, 1997. 

6. Cramer, S.M., and W.B. Fohrell, “Method for Simulating Tension Performance of 

Lumber Members,” Journal of Structural Engineering, Vol. 116, No. 10, October 

1990, pp. 2729–2747, figure 1, p. 2731. 

7. Goodman, James R., and Jozsef Bodig, “Orthotropic Strength of Wood in 

Compression,” Wood Science, Vol. 4, No. 2, October 1971, pp. 83–94, figure 6, 

p. 89; figure 9, p. 92. 

8. Stanzl-Tschegg, S.E., D.M. Tan, and E.K. Tschegg, “Fracture Resistance to the 

Crack Propagation in Wood,” International Journal of Fracture, Vol. 75, 1996, pp. 

347-356, figure 1, p. 349, figure 3, p. 350. 

9. Kretschmann, D.E., General Engineers, Forest Products Laboratory, Madison, 

WI. The Results from Trial of the CEN Draft Standard to Determine Fracture 

Energy for Tension Perpendicular to the Grain on Southern Pine. Investigates the 

proposed CIB W18A standards. Testing results presented at the 1990 meeting of 

CIB-W18A, figure 1, p. 11, figure 3, p. 15. 

10. Rohde, John, and John Reid, “Evaluation of the Performance Criteria for Wood 

Posts in Strong-Post W-Beam Guardrail,” TRB Paper 97-1206, Transportation 

Research Board, January 1997. 

11. Reid, S.R., and C. Peng, “Dynamic Uniaxial Crushing of Wood,” International 

Journal of Impact Engineering, Vol. 19, Nos. 5-6, 1997, pp. 531–570, figure 7, 

p. 545, figure 14, p. 560. 

12. Cramer, S.M., and J.R. Goodman, “Failure Modeling: A Basis for Strength 

Prediction of Lumber,” Wood and Fiber Science, Vol. 18, No. 3, 1986, pp. 446–

459, figure 4, p. 452. 

13. Green, D.W., and D.E. Kretschmann, “Moisture Content and the Properties of 

Clear Southern Pine,” Research Paper FPL-RP-531, Forest Products Laboratory, 

U.S. Department of Agriculture, 1994. 

14. Kretschmann, D.E., and David W. Green, “Modeling Moisture ContentMechanical Property Relationships for Clear Southern Pine,” Wood and Fiber 

Science, 28(3), 1996, pp. 320–337. 
148 

15. McBurney, R.S., and J.T. Drow, “The Elastic Properties of Wood―Young’s 

Moduli and Poisson’s Ratios of Douglas Fir and Their Relations to Moisture 

Content,” No. 1528-D, Forest Products Laboratory, U.S. Department of 

Agriculture, March 1956. 

16. Bodig, J., and B.A. Jayne, Mechanics of Wood and Wood Composites, Krieger 

Publishing Company, Malabar, FL, 1993, figure 313, p. 119, figure 7.32, p. 312, 

figure 11.28, p. 576. 

17. Woodward, Clinton, and John Minor, “Failure Theories for Douglas Fir in 

Tension,” Journal of Structural Engineering, Vol. 114, No. 12, December 1988, 

pp. 2808-2813. 

18. Forest Products Laboratory, U.S. Department of Agriculture, Wood Handbook: 

Wood as an Engineering Material, Agriculture Handbook 72, figure 4-14, p. 4–36. 

19. Patton-Mallory, Marcia, Steven M. Cramer, Frederick W. Smith, and Patrick J. 

Pellicane, “Nonlinear Material Models for Analysis of Bolted Wood Connections,” 

Journal of Structural Engineering, Vol. 123, No. 8, August 1997, pp. 1063–1070. 

20. Sandler, I., F.L. DiMaggio, and M.L. Barron, “An Extension of the Cap Model 

Inclusion of Pore Pressure Effects and Kinematic Hardening to Represent an 

Anisotropic Wet Clay,” in C.S. Desai and R.H. Gallagher, Mechanics of 

Engineering Materials, Chapter 28, Wiley, New York, NY, 1984. 

21. Simo, J.C., and J.W. Ju, “Stress- and Strain-Based Continuum Damage Models, 

Parts I and II,” International Journal of Solids and Structures, Vol. 23, No. 7, 

1987. 

22. Sih, G.C., P.C. Paris, and G.R. Irwin, “On Cracks in Rectilinearly Anisotropic 

Bodies,” International Journal of Fracture Mechanics, Vol. 1, No. 3, September 

1965, pp. 189–203. 

23. Bragov, A., and A.K. Lomunov, “Dynamic Properties of Some Wood Species,” 

Journal de Physique IV France 7, Colloque C3, Supplement au Journal de 

Physique III d’août [Journal of Physics IV France 7, C3 Conference, Supplement 

to the Journal of Physiques III for August], 1997, figure 1, p. C3-490. 

24. Murray, Y.D., “Modeling Rate Effects in Rock and Concrete,” Proceedings of the 

8th International Symposium on the Interaction of the Effects of Munitions With 

Structures, Defense Special Weapons Agency, McLean, VA, April 1997. 

25. Duvaut, G., and J.L. Lions, “Les Inequations en Mechanique et en Physique” 

[Equations in Mechanics and in Physics], Dunod, Paris, France, 1972. 

26. Standard Grading Rules for Southern Pine Lumber, Southern Pine Inspection 

Bureau, Pensacola, FL, 1991. 

27. Hashin, Z., “Failure Criteria for Unidirectional Fiber Composites,” Transactions of 

the American Society of Mechanical Engineers, Journal of Applied Mechanics, 

Vol. 47, June 1980, pp. 329–334. 

28. Nahas, Mahmoud N., “Survey of Failure and Post-Failure Theories of Laminated 

Fiber-Reinforced Composites,” Journal of Composites Technology and 

Research, Vol. 8, No. 4, Winter 1986, pp. 138-153. 

29. Tsai, Stephen W., and Edward M. Wu, “A General Theory of Strength for 

Anisotropic Materials,” Journal of Composite Materials, Vol. 5, January 1971, 

pp. 58–80. 
149 

30. Hill, R. The Mathematical Theory of Plasticity, Oxford University Press, Oxford, 

United Kingdom, 1950. 

31. Hoffman, O., “The Brittle Strength of Orthotropic Materials,” Journal of Composite 

Materials, Vol. 1, 1967, pp. 200–206. 

32. Tsai, Stephen W., and Victor D. Azzi, “Strength of Laminated Composite 

Materials,” AIAA Journal, Vol. 4, No. 2, February 1966, pp. 296–301. 

33. Yamada, S.E., and C.T. Sun, “Analysis of Laminate Strength and Its Distribution,” 

Journal of Composite Materials, Vol. 12, July 1978, pp. 275–284. 

34. Hankinson, R.L., “Investigation of Crushing Strength at Varying Angles of Grain,” 

Air Service Information Circular 3(259), Material Section Paper No. 130, 1921. 

35. Sandler, I.S., and T.A. Pucik, “Non-Uniqueness in Dynamic Rate-Independent 

Non-Associated Plasticity,” Mechanics of Materials and Structures, Elsevier 

Science, 1994, pp. 221–239. 
HRDS-04/08-07(1M)E

Recycled

Recyclable