An investigation of the effects of deformation-induced anisotropy on isotropic classical elastic-plastic materials

The effects of recoverable deformation induced anisotropy on the inelastic response of elastic materials are described. Starting by quantifying the degree of deformation induced anisotropy in isotropic materials, it is proved that the resultant anisotropy is significant only in materials capable of realizing large elastic deviatoric strains. For common engineering materials, this condition requires that the material strength increase strongly with pressure. For materials whose strength does not vary strongly with pressure, such as metals, recoverable deformation induced anisotropy is shown to be negligible.

AN INVESTIGATION OF THE EFFECTS OF DEFORMATION-INDUCED ANISOTROPY ON ISOTROPIC CLASSICAL ELASTIC-PLASTIC MATERIALS by Timothy J. Fuller A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy Department of Mechanical Engineering The University of Utah December 2010 Copyright c Timothy J. Fuller 2010 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of has been approved by the following supervisory committee members: , Chair Date Approved , Member Date Approved , Member Date Approved , Member Date Approved , Member Date Approved and by , Chair of the Department of and by Charles A. Wight, Dean of The Graduate School. Timothy J. Fuller Rebecca M. Brannon September 1, 2010 Daniel O. Adams September 1, 2010 Timothy Ameel September 1, 2010 James E. Guilkey September 1, 2010 O. Erik Strack September 1, 2010 Timothy Ameel Mechanical Engineering ABSTRACT The effects of recoverable deformation induced anisotropy on the inelastic re-sponse of elastic materials are described. Starting by quantifying the degree of deformation induced anisotropy in isotropic materials, it is proved that the resultant anisotropy is significant only in materials capable of realizing large elastic deviatoric strains. For common engineering materials, this condition requires that the material strength increase strongly with pressure. For materials whose strength does not vary strongly with pressure, such as metals, recoverable deformation induced anisotropy is shown to be negligible. In those materials that are capable of realizing large elastic strains, the effects of recoverable deformation induced anisotropy are revealed through the predicted coupling of hydrostatic and deviatoric responses in isotropic materials. It is shown that the coupling of the two responses is more significant than previously recognized in the literature. Properly accounting for the coupling of hydrostatic and deviatoric responses requires re-evaluating elastic materials characterization data, allowing for the coupled response. The result is an apparent decrease in the pressure sensitivity of the elastic shear modulus. The decrease in the pressure sensitivity of the shear mod-ulus leads to stress paths that are more tangential to the yield surface in stress space, resulting in an increase in predicted elastic strain at each step of an elastic-plastic stress update. Consequently, predicted plastic strains and, in particular, volumetric plastic strains, are smaller than if recoverable deformation induced anisotropy had been neglected, giving the appearance of a nonassociated plastic model. It is shown that this behavior agrees with what is experimentally observed. Numerical algorithms for the incorporation of recoverable deformation-induced anisotropy in existing classical elastic-plastic constitutive models are given. Up-grading existing code base to include recoverable deformation-induced anisotropy involves very few lines of extra coding and very little computational cost. Using the provided algorithms, model results for problems of interest to the geomechanics, defense, or any other engineering community where large pressures and deformations are characterized, are expected to be more predictive than if recoverable deformation induced anisotropy is neglected. iv To Kirsten, Miles, and Graham CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii ACKNOWLEDGMENTS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv PART I BACKGROUND AND MOTIVATION . . . . . . . . . . . . . . . . . 1 CHAPTER 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1 Motivation and Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Qualitative View of RDIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.1.1 Thought Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Overview of Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Importance of Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2. LITERATURE SURVEY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 RDIA in Geomechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 RDIA in Other Engineering Disciplines . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4 Research Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3. PRELIMINARIES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1 Notation and Tensor Representation . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2 Isotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 Tensor Valued Function Isotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4 Independent and Dependent Variables . . . . . . . . . . . . . . . . . . . . . . . . . 14 4. ELEMENTS OF THERMOMECHANICS . . . . . . . . . . . . . . . . . . . . 15 4.1 Conservation Laws in Thermomechanics . . . . . . . . . . . . . . . . . . . . . . . 15 4.2 Second Law of Thermodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.2.1 The Clausius-Duhem Inequality . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.3 Thermoelasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.4 Thermoplasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.4.1 Conservation Laws for Thermoplastic Materials . . . . . . . . . . . . . 20 4.4.2 Plastic Yield . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.4.3 Plastic Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4.4 Temperature Evolution Due to Dissipation . . . . . . . . . . . . . . . . . 26 PART II THERMODYNAMIC REQUIREMENT OF RDIA . . . . . 28 CHAPTER 5. THERMODYNAMIC REQUIREMENT OF RDIA . . . . . . . . . . . . 29 5.1 Elasticity in Isotropic Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5.2 Elastic Stiffness Isotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.3 Elastic Stiffness Isotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.4 First Order Approximation of C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.5 Elastic Compliance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 6. NUMERICAL IMPLEMENTATION OF RDIA . . . . . . . . . . . . . . . 34 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 6.2 Elastic Stress Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 6.3 Elastic-Plastic Stress Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 6.3.1 Elastic-Plastic Stress Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 6.3.2 Elastic Strain Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 6.3.3 Numerical Solution Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 6.3.4 Stress Update Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 PART III EFFECTS OF RDIA ON ISOTROPIC MATERIALS . . . . 40 CHAPTER 7. QUANTIFICATION OF ANISOTROPY . . . . . . . . . . . . . . . . . . . . . 41 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 7.2 Scalar Measure of RDIA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 7.2.1 Properties of Anisotropy Measure . . . . . . . . . . . . . . . . . . . . . . . . 44 7.3 Evaluation of Degree of Anisotropy in Isotropic Materials . . . . . . . . . . 44 7.4 Degree of Anisotropy in Anisotropic Materials . . . . . . . . . . . . . . . . . . 45 7.5 Degree of RDIA in Isotropic Materials . . . . . . . . . . . . . . . . . . . . . . . . 45 7.5.1 Constitutive Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.5.1.1 Yield Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.5.1.2 Elastic Moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 7.5.2 Degree of RDIA in Elastic Materials . . . . . . . . . . . . . . . . . . . . . . 47 7.5.3 Degree of RDIA in Elastic-Plastic Materials . . . . . . . . . . . . . . . . 47 7.5.3.1 Factors Affecting the Degree of RDIA . . . . . . . . . . . . . . . . . 47 7.5.3.1.1 Pressure-dependence of the shear modulus. . . . . . . . . . 48 7.5.3.1.2 Pressure-dependence of the yield. . . . . . . . . . . . . . . . . 48 7.6 Magnitude of RDIA in Engineering Materials . . . . . . . . . . . . . . . . . . . 48 7.7 Discussion and Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 vii 8. EFFECTS OF RDIA ON ISOTROPIC MEDIA . . . . . . . . . . . . . . . 63 8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 8.2 Elastic Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 8.2.1 Coupled Hydrostatic-Deviatoric Response . . . . . . . . . . . . . . . . . . 66 8.2.2 Elastic Moduli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 8.2.3 Elastic Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 8.3 Inelastic Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 8.3.1 Trial Stress Increment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8.3.2 Return Direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 8.4 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 8.4.1 Uniaxial Strain Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 8.4.1.1 Elastic Limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 8.4.1.2 Trial Stress Increments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 8.4.1.3 Permanent Volume Change . . . . . . . . . . . . . . . . . . . . . . . . . 73 8.4.2 Triaxial Strain Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 8.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 9. CASE STUDY: EFFECTS OF RDIA . . . . . . . . . . . . . . . . . . . . . . . . 90 9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 9.2 Constitutive Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 9.3 Constitutive Model Parametrization . . . . . . . . . . . . . . . . . . . . . . . . . . 90 9.4 Computational Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 9.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 9.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 PART IV DISCUSSION AND CONCLUSION . . . . . . . . . . . . . . . . . . . 98 CHAPTER 10. DISCUSSION AND CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . 99 10.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 10.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 APPENDIX: NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 viii LIST OF FIGURES 1.1 In the compressed state, the stiffness of the isotropic collection of foam spheres is anisotropic. Deformation induces the anisotropy. . . . . . . . . . . 7 4.1 Oblique return projection of the trial stress state on to the yield function isosurface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.1 Elastic Kirchhoff stress-update algorithm. New terms involving RDIA, not expected to be included in an existing stress-update algorithm, are set apart explicitly by • in Step 1. . . . . . . . . . . . . . . . . . . . . . . . . . . 38 6.2 Elastic-plastic Kirchhoff stress-update algorithm. New terms involving RDIA, not expected to be included in an existing stress-update algo-rithm, are set apart explicitly by • in Steps 1, 3, and 5. . . . . . . . . . . 39 7.1 Geometric interpretation of decomposition of fourth-order tensor into isotropic and nonisotropic parts, showing the angle θ between the tensor A and the isotropic hyper-plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.2 Measure of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic symmetry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (mono-clinic symmetry), and Copper (cubic symmetry). Material data from B¨ohlke and Br¨uggemann (2001). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 7.3 Variation of the measure of anisotropy for elastic uniaxial strain de-formation ×+. Also shown are the measures of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic symmetry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (monoclinic symmetry), and Copper (cubic symmetry). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 7.4 Variation of the measure of anisotropy for uniaxial strain deformation for the same material as in Figure 7.3 but allowing for plastic yielding ×+. Also shown are the measures of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic symmetry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (monoclinic symmetry), and Copper (cubic sym-metry). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 7.5 Variation of the measure of anisotropy with μ for uniaxial strain defor-mation plotted against axial strain for dμ/dp = 50 −·−·−·−·−, dμ/dp = 20 ············, dμ/dp = 10 −−−−−−, dμ/dp = 2 ----, and dμ/dp = 0.0 ----. For large values of μ , the measure of anisotropy decreases due to the limiting value of shear strain. . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 7.6 Variation of the measure of anisotropy with a4 for uniaxial strain de-formation plotted against axial strain for a4 = 0.399 −·−·−·−·−, a4 = 0.2 ············, a4 = 0.1 −−−−−−, a4 = 0.05 ----, and a4 = 0.0 ----. Except for the largest values of a4, the measure of anisotropy never exceeds 5%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 7.7 Degree of anisotropy in elastic-plastic materials in uniaxial strain com-pression for linearly pressure-dependent elastic moduli for 99.5 % alu-mina ············, TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. Note that, with the exception of the powdered metal (alu-mina), the degree of recoverable anisotropy never exceeds 6%. In the small to medium strain regime, the degree of recoverable anisotropy does not exceed 2%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.8 Degree of recoverable anisotropy in elastic-plastic materials in uniaxial strain compression for elastic moduli given by Mie Gr¨uneisen equation of state for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. The degree of recoverable anisotropy never exceeds 4% in the small to medium strain regime for any of the materials. . . . . . . . . . . . . . . . . . . 60 7.9 Degree of recoverable anisotropy in elastic-plastic materials in triaxial strain compression for linearly pressure-dependent elastic moduli for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. With the exception of alumina, the degree of recoverable anisotropy never exceeds 5% even in the finite strain regime. . . . . . . . . . . . . . . . . . . . . . 61 7.10 Degree of recoverable anisotropy in elastic-plastic materials in triaxial strain compression for elastic moduli given by Mie Gr¨uneisen equation of state for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. Interestingly, the degree of recoverable anisotropy in limestone exceeds 8% when the bulk modulus is computed from the Mie Gr¨uneisen equa-tion of state while the degree of recoverable anisotropy does not exceed 6% and 1% for TiB2 and 2024 aluminum, respectively. . . . . . . . . . . . . 62 8.1 Difference in the z and r components of the stress at the elastic limit for titanium diboride. The yield zero isosurface of the yield function is denoted by ----, the admissible response by ---, and inadmis-sible response by −··−··−··−. Allowing for RDIA predicts an increase in pressure and, thus, an increase in shear stress at the elastic limit. . . . . . 79 8.2 Comparison of the trial stress increment for the admissible model--- and the inadmissible model −·−·−·−·− in uniaxial compression, showing that the increment in the z component of the admissible trial stress is greater than the corresponding increment in the inadmissible trial stress while the increment in the r component of the trial stress in the admissible model is smaller than the corresponding component in the inadmissible model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 x 8.3 Detail of the comparison of the trial stress increment for the admissible model --- and the inadmissible model −·−·−·−·− in uniaxial compres-son. In this more detailed figure, it is obvious that the increment in the z component of the admissible trial stress is greater than the corre-sponding increment in the inadmissible trial stress while the increment in the r component of the trial stress in the admissible model is smaller than the corresponding component in the inadmissible model. . . . . . . . 81 8.4 Predicted permanent volume strain in the admissible and inadmissible models in uniaxial compression. At 200 MPa, the admissible model predicts nearly 24% less plastic volume strain than the inadmissible model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 8.5 Predicted permanent volume strain in the admissible and inadmissible models using associative and nonassociative flow rules in uniaxial com-pression. At 200 MPa, the admissible model predicts 24% less plastic volume strain than the inadmissible model for associative flow. For nonassociative flow, the flow potential parameter β is 22% larger in the admissible model than in the inadmissible model. . . . . . . . . . . . . . . . . . 83 8.6 Ratio of plastic work in the admissible Πp and inadmissible model ˜Π p from uniaxial compression. The plastic work is less in the admissible model than the inadmissible model in the pressures investigated. . . . . . 84 8.7 Comparison of the trial stress increment for the admissible model--- and the inadmissible model −·−·−·−·− showing that the increment in the z component of the admissible trial stress is greater than the correspond-ing increment in the inadmissible trial stress while the increment in the r component of the trial stress in the admissible model is smaller than the corresponding component in the inadmissible model. . . . . . . . . . . . . 85 8.8 Detail of the comparison of the trial stress increment for the admissible model --- and the inadmissible model −·−·−·−·−. In this more de-tailed figure, it is obvious that the increment in the z component of the admissible trial stress is greater than the corresponding increment in the inadmissible trial stress while the increment in the r component of the trial stress in the admissible model is smaller than the corresponding component in the inadmissible model. . . . . . . . . . . . . . . . . . . . . . . . . . . 86 8.9 Predicted permanent volume strain in the admissible and inadmissible models. At 200 MPa, the admissible model predicts nearly 50% less plastic volume strain than the inadmissible model. . . . . . . . . . . . . . . . . 87 8.10 Predicted permanent volume strain in the admissible and inadmissible models using associative and nonassociative flow rules. At 200 MPa, the admissible model predicts 50% less plastic volume strain than the inadmissible model for associative flow. For nonassociative flow in the inadmissible model, the predicted permanent volume strain is still larger than in the associated admissible model. . . . . . . . . . . . . . . . . . . . . . . . . 88 xi 8.11 Ratio of plastic work in the admissible Πp and inadmissible model ˜Π p. The plastic work is less in the admissible model than the inadmissible model in the pressures investigated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 9.1 Schematic for shock compression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 9.2 Degree of anisotropy in shock loaded 6061 T6 aluminum. The degree of anisotropy is smaller than the noise in the computational simulation, never exceeding 0.4%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 9.3 Axial and lateral stress predictions at the back surface of the target. In shock compression, the stress state is nearly isotropic, thus the two stress responses are nearly indistinguishable. The difference between the admissible and inadmissible model is negligible. . . . . . . . . . . . . . . . 96 9.4 Longitudinal velocities at the back surface of the target in the admissible model ----, inadmissible model −·−·−·−·−, and experimental data ---. The simulated velocities are not in good agreement with the experimental data after the Hugoniot state. However, the results show no perceivable difference between the admissible and inadmissible models. 97 xii LIST OF TABLES 7.1 Anisotropic elastic stiffness coefficients of several inherently anisotropic materials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 7.2 Material properties representative of ceramic, geologic, and metallic materials used in simulations shown in Figures 7.7 - 7.10. . . . . . . . . . . . 58 8.1 Material properties adopted in this study. Material properties are sim-ilar to those in Salari et al. (2004). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 9.1 Material properties for the 6061 T6 aluminum used in this study. . . . . . 94 ACKNOWLEDGMENTS I would like to thank the members of my committee, who have each contributed in some way to me getting where I am today. I would especially like to thank those professors at the University of Utah: Professor Timothy Ameel, who served as the graduate school adviser when I first enrolled and assisted me in finding an adviser and area of study, Professor Dan Adams who served as my adviser for the first year and a half of my graduate studies , and Professor Jim Guilkey for the weekly discussions in the last several months that have helped me keep things in proper perspective. I am especially grateful to Dr. O. Erik Strack from Sandia National Laboratory for taking a chance and providing the resources necessary for me to complete graduate studies. His counsel, advice, and technical knowledge have helped considerably over the last three and a half years. Of course, this list is incomplete without mentioning Professor Rebecca Brannon who has served as my adviser and committee chair since her first semester at the University of Utah in 2007. From the first time I heard her keys jangling from her ever present fanny pack, and attended her introductory course on continuum mechanics, I have become passionate about computational mechanics in all of its aspects. Her passion and expertise in seemingly every field of mechanics have been invaluable to me over the last three and a half years. Lastly, I would like to thank my beautiful wife, Kirsten. When we first met, she had no way of knowing that the hopeful professional alpine ski racer she met and married would decide to begin university studies years after most of his friends had finished, and then proceed on to doctoral studies. Her patience and love are immeasurable, and I would not be half of the person I am today without her. PART I BACKGROUND AND MOTIVATION CHAPTER 1 INTRODUCTION 1.1 Motivation and Problem Statement In recent years, considerable attention has been paid to advanced constitutive models that attempt to account for the dependence of material strength and elastic properties on pressure. Classes of pressure sensitive materials of interest include geologic media, protective armor ceramics in the defense community, and compliant polymers for explosives detonation. In these materials, the constitutive response is sensitive to pressure, and accurately predicting material response to large deforma-tions at high pressures is of great research interest. Often overlooked in model development is the elastic component of the gener-ally elastic-plastic constitutive response, with most attention paid to the inelastic response. Relatively little attention has been paid to the mathematical and physical well posedness of the elastic component. The elastic component, however, of a general elastic-plastic model is an important element in developing a model that is predictive under general loading and, in fact, serves as the foundation from which the plastic component is developed. Accordingly, the assumptions used when interpreting the elastic component of materials characterization data must be thermodynamically admissible to ensure greater accuracy of the general model. However, in the literature, it is common for thermodynamic admissibility to be inadvertently abandoned in favor of the adoption of certain simplifying assumptions and interpretations of experimental data. For adiabatic large distortions of an elastic material, commonly adopted assump-tions used in the modeling of shock data are: 1. Stress is derivable from an elastic energy potential 3 2. Material is isotropic 3. Elastic stiffness is isotropic 4. Shear modulus varies with pressure (or is nonconstant) Historically, for small strain analysis, the first three assumptions have been adopted. For larger strains, material nonlinearity is significant, and the last assumption is necessarily adopted. However, the simultaneous adoption of a nonconstant shear modulus and an isotropic elastic stiffness conflicts with the first assumption, resulting in a thermodynamically inadmissible model (Zytynski et al., 1978). The inadmissibility arises from the erroneous assumption that material isotropy implies isotropy of the elastic stiffness tensor. However, it is well known that ther-modynamic considerations require that isotropic materials develop recoverable defor-mation induced anisotropy (RDIA) in their elastic stiffnesses, except in the special cases of isotropic deformation, infinitesimal distortional deformation, or indifference of elastic moduli to loading (constant elastic moduli) (Zytynski et al., 1978; Simo and Pister, 1984; Marsden and Hughes, 1994). Thus, the traditional approach to nonlinear elasticity of isotropic materials, to place material nonlinearity on the isotropic elastic moduli, is thermodynamically inadmissible. For nonlinear deformations, thermody-namic admissibility is recovered by allowing the stiffness tensor to develop RDIA, thereby altering its integrity basis. The commonly accepted traditional approach to nonlinear elasticity notwithstand-ing, the requirement of RDIA in the elastic stiffness of isotropic materials is known and documented in well respected texts on applied mechanics (cf. Marsden and Hughes, 1994; Simo and Hughes, 2000). However, RDIA has received very little attention in the research literature, particularly in the context of elastic-plastic deformation where it is often assumed that the elastic contribution to the constitutive response is negligible in comparison to the plastic portion during plastic loading cycles. To the contrary, this dissertation, which is the first investigation into the effects of RDIA on isotropic elastic-plastic materials, shows that RDIA in the elastic response of 4 pressure-dependent isotropic elastic-plastic materials can have as much as first-order implications in constitutive predictions such as plastic strain and work. 1.1.1 Qualitative View of RDIA Not to be confused with inherent anisotropy, such as in fiber reinforced composite materials or single crystals, RDIA in isotropic materials is induced by loading and recovered upon unloading. Though seemingly counter-intuitive, the following thought experiment demonstrates RDIA. 1.1.1.1 Thought Experiment Imagine a collection of uniformly distributed foam spheres, as shown in Fig-ure 1.1(a). In this undeformed state, the collection generally possesses cubic symme-try. Imagine that the axial and lateral stiffness, Kx and Ky, respectively, are equal. When loaded axially, elastic deformation causes the spheres to deform into flattened ellipsoids, as in Figure 1.1(b). In response to this deformation, the stiffness increases in the direction of loading as the material is compressed. In the lateral direction, there would be relatively little change in stiffness. Thus, in the deformed state Kx > Ky, and the collection of spheres transforms from one of cubic to orthotropic symmetry. For a body with isotropic symmetry, deformation generally induces anisotropy in the elastic stiffness. Anisotropy induced by deformation is consistent with isotropy be-cause the body can be arbitrarily rotated prior to loading and the resulting anisotropy would be the same. Since, upon unloading, the body returns to their original spherical shapes, the anisotropy is recoverable, or elastic. 1.2 Overview of Study As motivated by the previous thought experiment, and to be proved in the follow-ing chapters, thermodynamics requires that isotropic materials develop deformation-induced anisotropy in their elastic stiffnesses in response to deformation. Though documented in well respected texts on applied mechanics (Malvern, 1969; Marsden and Hughes, 1994; Simo and Hughes, 2000) and recognized in fields such as rheology and geomechanics (Scheidler, 1996; Hueckel et al., 1992; Johnson and Rasolofosaon, 5 1996) where large material distortions are commonly observed, RDIA has seen little to no consideration in other fields, such as shock physics, where stiff materials undergo large distortions thru shock processes and elastic moduli are inferred from measured wavespeeds. Additionally, the extent to which deformation-induced anisotropy affects isotropic materials in general has not been recognized in the literature. This disser-tation is the first systematic analytical investigation of the effects of RDIA in the elastic portion of an elastic-plastic model on the overall inelastic material response in general loading. The topics addressed in this dissertation are • The conditions necessary for an isotropic material to have an isotropic stiffness • Thermodynamically admissible inclusion of RDIA in existing classical elastic-plastic stress-update frameworks. • The magnitude of RDIA in isotropic materials. • The conditions and materials for which RDIA is significant. • The effects of RDIA on the classical elastic-plastic stress-update. In the process of addressing the above topics, guidance for the type of materials most affected by RDIA is given along with computational algorithms for the inclusion of RDIA in a classical elastic-plastic stress-update framework. 1.3 Importance of Study Thermodynamic admissibility in the modeling of continuous media has been the subject of literally thousands of publications. In the study of materials loaded to the elastic limit, thermodynamic considerations have given rise to Green elasticity, hyperelasticity, hypoelasticity, equation of state development, among others, where the stress is directly related to the overall deformation. Beyond the elastic limit, path dependence of material response requires that the direct stress-deformation relationships be replaced by incremental relationships, requiring knowledge of the 6 fourth-order elastic tangent1 stiffness tensor. For isotropic materials, the elastic stiffness is almost always presumed to be itself isotropic, and material nonlinearity is modeled through nonlinearity of elastic moduli. However, this assumption is in conflict with the basic tenants of thermodynamics, which require RDIA to develop in the elastic stiffness of an isotropic materials. In recent years, the development of hyperplasticity has implicitly included RDIA in some models, though its effects have not been described in the literature beyond passing mention (Houlsby et al., 2005; Collins and Houlsby, 1997; Hueckel et al., 1992). In this dissertation, the effects of RDIA on isotropic elastic-plastic materials are described for the first time. An important result of the study is that the degree of RDIA in materials whose shear strength does not vary significantly with pressure, such as metals, is negligible. However, in materials whose shear strength and elastic moduli do vary significantly with pressure, such as rocks, the degree of RDIA is not negligible and model predictions accounting for RDIA differ in the first order from those that do not. Algorithms for the inclusion of RDIA in an elastic-plastic stress-update are provided. The necessary source code revisions add little computational complexity to existing codes. 1From here forward, "tangent" will be understood. 7 x y Kx = Ky (a) Undeformed state Kx > Ky (b) Deformed state Figure 1.1. In the compressed state, the stiffness of the isotropic collection of foam spheres is anisotropic. Deformation induces the anisotropy. CHAPTER 2 LITERATURE SURVEY 2.1 Introduction When subjected to large distortions, initially isotropic materials respond by devel-oping anisotropy in their response to deformation. Well-known examples of inelastic deformation-induced anisotropy (IDIA) include directional anisotropy in cold worked metals (Hill, 1948), plastic flow induced anisotropy (Stouffer and Bodner, 1979), and nucleation and growth of oriented micro-cracks in brittle media (Horii and Nemat- Nasser, 1983). Not to be confused with inherent anisotropy, such as in fiber reinforced composite materials or single crystals, deformation-induced anisotropy is caused by loading and is manifest in changes in the isotropy of the fourth-order tangent stiffness tensor of the material. In addition to the above mentioned examples of nonrecoverable deformation-induced anisotropy, the first and second laws of thermodynamics imply that the fourth-order tangent stiffness tensor of initially isotropic media develop generally recoverable deformation-induced anisotropy (RDIA). This consequence of thermody-namics is known in the applied mechanics community. For example, in his text The Mathematical Foundations of Elasticity, Marsden and Hughes give the most general form of the elastic stiffness tensor in terms of the right Cauchy-Green stretch tensor (Marsden and Hughes, 1994, p. 23) which is, in general, an anisotropic fourth order tensor. Though Marsden and Hughes do not explicitly state this subtle fact, it is left as an exercise to readers to derive the conditions under which the stiffness is isotropic. Those conditions, also given by Hueckel for the secant stiffness of isotropic materials (Hueckel et al., 1992) and Simo and Pister (1984) for generally hyperelastic materials, are that the bulk modulus vary at most with the first invariant of the strain tensor and 9 the shear modulus be constant. Willam (2002), also describes the anisotropy in the elastic stiffness of isotropic materials resulting from thermodynamic considerations. The violation of thermodynamic principles resulting from the adoption of a non-constant shear modulus in an isotropic elastic stiffness was first described by Zytynski et al. (1978) who showed that if the shear modulus of elastic materials is allowed to vary with pressure, as commonly assumed in geomaterials, it is possible to construct elastic cycles closed in stress that are not closed in strain. If deformed in a particular manner, a net energy increase is possible - a clear violation of the first law of thermodynamics for elastic adiabatic loading. However, Zytynski asserted that his analysis lead to one of two alternatives: (1) adopt a constant shear modulus and allow Poisson's ratio to vary with deformation, allowing for Poisson's ratio to become negative, or (2) abandon the notion of an elastic region of material behavior. Zytynski failed to consider a third alternative: allowing the elastic stiffness to develop reversible anisotropy in response to deformation. 2.2 RDIA in Geomechanics Houlsby (1985) describes the theoretical difficulties of modeling pressure-dependence of the shear modulus in clays in a thermodynamically admissible manner. Particular attention is paid to the fact that coupling of isotropic and deviatoric responses through the pressure-dependence of the shear modulus also requires dependence of the bulk modulus on the stress deviator. In particular, Houlsby notes that pressure-dependence of the shear modulus cannot be simplistically achieved through pressure-dependence of the bulk modulus and constancy of Poisson's ratio. Instead, thermodynamics requires additional terms, in the form of deformation-induced anisotropy. Isotropic-deviatoric coupling in isotropic materials in the context of RDIA has also been described by Hueckel et al. (1992) and Scheidler (1996). Using the form of pressure-dependent shear modulus proposed by Wroth (1971), Hueckel uses reci-procity relations to derive the shear stress dependence of the bulk modulus. Triaxial experimental records on consolidated clays imply the need for pressure-dependence of the shear modulus, though the shear dependence of the bulk modulus is considered 10 negligible in comparison and is omitted in their analysis. However, using analytic methods, Scheidler showed that neglecting the shear dependence of the bulk response of isotropic materials could lead to overestimation of the apparent Hugoniot elastic limit in uniaxial strain loading. Scheidler claims that for Titanium Diboride, TiB2, the difference could be as much as 18%, though the claim is not yet validated. More recently, the theory of hyperplasticity (Borja et al., 1997; Collins and Houlsby, 1997; Einav and Puzrin, 2004; Houlsby et al., 2005) has been developed in which a yield function is not explicitly defined but, instead, a dissipation potential is given from which dissipative inelastic behavior is derived, similar to how elastic behavior is derived from an elastic energy potential. In this theory, RDIA plays an implicit role through cross anisotropic terms in the elastic stiffness and compliance (Borja et al., 1997; Collins and Houlsby, 1997; Einav and Puzrin, 2004). Houlsby et al. (2005) has proposed a hyperelastic framework for modeling the stress-dependent nonlinear elastic stiffness of soils using this technique. Results from hyperplasticity have been shown to match experimental loading/unloading data suggestive of pressure dependence of the shear modulus in undrained clays more closely than the classical approach to plasticity (Borja et al., 1997). 2.3 RDIA in Other Engineering Disciplines Outside of the geomechanics community, evidently RDIA has seen little, if any, consideration. In metals plasticity, for example, material nonlinearity is observed (Richmond and Spitzig, 1980; Burakovsky et al., 2003) and based on the observation that the ratio of flow strength to shear modulus is approximately constant, pressure dependence of the shear modulus is inferred from experimental data (Guinan and Steinberg, 1975; Hua et al., 2002). Pressure-dependence of the shear modulus is also inferred from recordings of ultrasonic shear wave speeds taken at varying pressures. However, material data are still interpreted in the context of elastic stiffness isotropy (Hayes et al., 1999). Evidently, RDIA has never been considered in the context of metals plasticity. In the study of metal powder compaction, large pressures and material distortions 11 are needed to form powdered metals into desired shapes. Achieving constant density gradients through the formed materials is largely a process of trial and error, though finite-element analysis and engineering principles help reduce the number of trials. However, material response is still widely assumed to be governed by an isotropic stiffness (Zeuch et al., 2001), despite extremes of pressure and deformation. In other words, thermodynamics has been implicitly abandoned by not accounting for RDIA. 2.4 Research Scope Despite being a known consequence of thermodynamics for the last three decades, and being included implicitly in many geomechanics analysis, the effects of RDIA on classical elastic-plasticity theory have not been considered in the engineering litera-ture. In this dissertation, the effects of RDIA on components of classical plasticity theory are quantified. In materials whose strength and elastic moduli vary with pressure, it is shown that the degree of anisotropy can obtain magnitudes that are significant. An immediate effect of RDIA is that isotropic and deviatoric components of the constitutive response of isotropic materials to deformation are coupled. Though this coupling has previously been mentioned in the research literature, its effects beyond the elastic limit have not. Properly accounting for the coupling of isotropic and deviatoric responses alters the interpretation of experimental data suggestive of pressure varying shear moduli so that higher order shear modulus terms appear to be smaller than if the coupling is neglected. As a consequence, any term acted on by the shear modulus must by adjusted to compensate. Accounting for the coupling of isotropic and deviatoric responses in the elastic response of isotropic materials is described. Guidelines for the for the computational adoption of RDIA in a classical elastic-plastic stress update framework are given. CHAPTER 3 PRELIMINARIES 3.1 Notation and Tensor Representation Focusing on applications in mechanics, first, second, and fourth-order tensors are R3. Components of tensors in R3 are defined relative to an orthonormal basis (e1, e2, e3) so that a = aiei, A = Aijeiej , and A = Aijkleiejekel, where implied summation is assumed from 1 - 3. Additionally, second-order and fourth-order tensors are presumed symmetric and minor symmetric, respectively, and thus are also equivalently cast as first and symmetric second-order tensors in R6 (Helnwein, 2001). Components of tensors in R6 are defined relative to an orthonormal basis (E1,E2, . . . ,E6), thus A = AIEI and A = AIJEIEJ , with implied summation from 1 - 6. The ei and EI are related by EI = eiei, i= 1, 2, 3, E4 = 1 √2 (e1e2 + e2e1) , E5 = 1 √2 (e2e3 + e3e2) , E6 = 1 √2 (e3e1 + e1e3) . (3.1) In the above, adjacent basis tensors represent dyads and a raised dot · between tensor arguments represents the inner product of a pair of basis tensors, e.g., ei·ej = δij , where δij is the Kronecker delta, a·b = aibjei·ej = aibi, A:B = AijBij, etc. Linear mapping, composition, and inner product of tensors of arbitrary order are thus constructed by the appropriate number of raised dots between arguments. The L2, or Frobenius, norm of a tensor of arbitrary order A is denoted A and is defined relative to the orthonormal Euclidean components by A 2 = i,j,...,m=1 nAij...mAij...m, (3.2) 13 where m is the order of A ∈ Rn. In other words, this norm is the square root of the sum of the squares of components with respect to any underlying orthonormal basis, whether R3, R6, or any other dimension. 3.2 Isotropy Definition 3.2.1. An isotropic second-order tensor A is expressible as A = Q∗A = √3αAEz ∀Q ∈ SO(3), αA = 1 3 trA, (3.3) where SO(3) is the set of all proper orthogonal1 second-order tensors in R3, Ez = δ/√3, δ is the second-order identity tensor, and ∗ is the Rayleigh product, defined for second-order tensors as Q∗A = Q·A·QT (3.4) For a tensor of arbitrary order, the Rayleigh product is defined analogously as (Q∗A)ij···kl = i=1 3Qiq · · ·QjnAn···qo···pQkp · · ·Qlo. (3.5) Definition 3.2.2. A minor-symmetric isotropic fourth-order tensor A is expressible as A = Q∗A = αAEiso + √5βAEsd, ∀Q ∈ SO(3), (3.6) where α and β are the eigenvalues of A and the orthonormal eigenprojectors of A, Eiso and √5Esd, are given in indicial form by Eiso = 1 3 δijδkl, √5Esd = 1 2 (δikδjl + δilδjk) − 1 3 δijδkl. (3.7) When acting on a second-order tensor A, Eiso and √5Esd return the isotropic part of A and the symmetric deviatoric part of A, respectively. The factor of √5 is introduced 1In general, the restriction Q ∈ SO(3) is overly restrictive and can be relaxed to Q ∈ O(3), where O(3) is the set of all orthogonal tensors. However, physical arguments lead to the restriction Q ∈ SO(3) for physical tensors common to engineering mechanics, thus the restriction. 14 so that both Eiso and Esd are unit tensors, thus making them an orthonormal basis for the linear manifold of all minor-symmetric fourth-order tensors. This fact is later exploited to quantify the degree of anisotropy. 3.3 Tensor Valued Function Isotropy Definition 3.3.1. A tensor valued function h(A) is isotropic if h(Q∗A) = Q∗h(A) ∀Q ∈ SO(3). (3.8) By the representation theorem (Smith, 1971; Truesdell et al., 2004), h(A) is express-ible as h(A) = ω1δ + ω2A + ω3 dev (A ·A ), (3.9) where A is the deviatoric2 part of A , the dev operator returns the deviatoric part of its argument, and the ωi are scalar functions of the "mechanics" invariants of A defined as JA 1 = trA, JA 2 = 1 2 trA ·A , JA 3 = 1 3 trA ·A ·A . (3.10) Alternatively, the ωi can be chosen as functions of any other set of independent invariants of A, in which case the tensor basis can also change to linearly independent triplet that span tensors commuting with A. 3.4 Independent and Dependent Variables Throughout the entirety of this work, strain is regarded as the independent vari-able. Thus, constraints are developed on the elastic stiffness and moduli on volume change and distortion. By the chain rule, these constraints can be readily cast in terms of stress. Inversion of the elastic stiffness gives similar constraints on the elastic compliance. 2Superscripted primes are used to both denote the deviatoric part of a tensor argument and differentiation with respect to the primary variable for a function of a single variable. In cases where context alone does not distinguish between the two, explicit differentiation notation will be adopted. CHAPTER 4 ELEMENTS OF THERMOMECHANICS The conservation laws developed in classical thermomechanics form the basis on which all other physical laws of continuum mechanics are derived. As such, a review of classical thermomechanics is an appropriate starting point for the topics considered in later chapters. Rather than provide a comprehensive overview of the discipline, only the notation and key concepts pertinent to later material will be considered. For the reader comfortable with the discipline of classical thermomechanics, in particular, its relationship to thermoelasticity and thermoplasticity, this chapter can be skipped without loss of continuity. A more detailed treatise on classical thermomechanics may be found in seminal works by Malvern (1969), Truesdell et al. (2004), and Gurtin (2003). 4.1 Conservation Laws in Thermomechanics Given an extensive quantity contained in an enclosed domain Ω, the rate of change of that quantity must be equal to the sum of the production of the quantity within the domain and the flux of the quantity across the boundary of the domain ∂Ω. Mathematically, conservation laws can be expressed in the following general form d dt Ω f(x, t) dV = ∂Ω f(x, t) (vn(x, t) − u˙ (x, t)·n(x, t)) dA + ∂Ω g(x, t) dA + Ω h(x, t) dV, (4.1) where f is a scalar, vector, or tensor valued conserved quantity, vn is the normal velocity of the boundary ∂Ω, x is the spatial position vector, u˙ is the material velocity, n is the outward unit normal to ∂Ω, g is the surface source of f, h is the volume source of f, and t the time, respectively. 16 Using Eq. (4.1), the conservation of mass, momentum, and energy can be written in local form as ρ˙ + ρ ∇ ·u˙ = 0, (4.2) ∇ ·σ + ρb = ρa, (4.3) u˙ − Jσ:d + J ∇ ·q − ρ0r = 0, (4.4) where ρ is the material density, σ is the Cauchy stress, b is the body force per unit mass, a is the material acceleration, u is the internal energy per unit reference volume, d is the symmetric part of the velocity gradient, J is determinant of the deformation gradient F, q is the heat flux vector, and r is the energy production per unit mass. For shock loading, Eq. (4.1) continues to apply and leads to additional Rankine-Hugoniot jump conditions that supplement the above local differential equations (Drumheller, 1998). 4.2 Second Law of Thermodynamics The internal energy is commonly regarded as a state function of the other state variables; the general form of this function is restricted by the second law of ther-modynamics. For any reversible process, the integral of the heat divided by the temperature is zero over any cycle closed in the thermodynamic state variables. ¯ dQ θ = 0. (4.5) Thus, even though the heat increment ¯ dQ is an inexact differential, dividing it by the temperature θ produces an exact differential in reversible loading. Accordingly, there must exist a state variable, S, such that, for reversible processes, dS = ¯ dQ θ , (4.6) where S is the total entropy of the system. It is established experimentally that the entropy change for isolated (Q = 0) systems is never negative and reaches its maximum at equilibrium. This experimental fact is known as the second law of 17 thermodynamics. For nonisolated systems, the entropy change for the system and its surroundings is always non-negative. Mathematically, we can write the second law as ˙S = ˙S sys + ˙S env ≥ 0 (4.7) where ˙S sys and ˙S env are the change in entropy of the system and the surrounding environment, respectively. Moreover, entropy is an extensive quantity, thus implying existence of a specific entropy, η. 4.2.1 The Clausius-Duhem Inequality Revising the balance law in Eq. (4.1) to allow for imbalance, the second law of thermodynamics can be written d dt Ω ρη dV ≥ ∂Ω ρη (vn(x, t) − u˙ (x, t)·n(x, t)) dA − ∂Ω q·n θ dA + ∂Ω ρr θ dV. (4.8) In local form, Eq. (4.8) becomes ρ0η˙ ≥ −J ∇ · q θ + ρ0 r θ , (4.9) which is known as the Clausius-Duhem inequality or the entropy inequality. 4.3 Thermoelasticity Expanding the first term on the right hand side of Eq. (4.9), the Clausius-Duhem inequality becomes ρ0η˙ ≥ − J θ ∇ ·q + J θ2 q· ∇ θ + ρ0 r θ (4.10) Assuming that θ > 0, substituting the balance of energy, Eq. (4.4), into Eq. (4.10), the Clausius-Duhem inequality may be written ρ0θη˙ ≥ u˙ − Jσ:d + J θ q· ∇ θ. (4.11) The Clausius-Duhem inequality may also be written as a sum of internally dissipative and heat conductive parts D +F ≥0. (4.12) 18 where the internally dissipative part D is given by D = Jσ:d + ρ0θη˙ − u˙ , (4.13) and the heat conductive part by F = − 1 ρθ q· ∇ θ. (4.14) For large deformations, where d is not a true rate of strain, it is customary to introduce alternative reference conjugate stress and deformation measures, σ, , such that the stress power is given by Jσ:d = σ:˙ , (4.15) where is a reference elastic strain, ˙ is its true rate, and σ is its work conjugate stress. If, for example, the strain is the deformation gradient, then σ is the first Piola-Kirchhoff stress , or, if the strain is the Green-Lagrange then σ is the second Piola-Kirchhoff stress. It is essential that the stress-strain measures be reference tensors (i.e., unaffected by a change in frame) to ensure that the rate ˙ can be treated as a true rate. Using Eq. (4.15) and the generalized stress and strain measures, Eq. (4.13) be-comes D = σ:˙ + ρ0θη˙ − u˙ . (4.16) In thermoelasticity it is presumed that the internal energy is a function of the deformation through the strain tensor and the entropy. In this case, using the chain rule of differentiation, u˙ = ∂u ∂ :˙ + ∂u ∂η η˙, so that Eq. (4.16) can be written D = σ − ∂u ∂ :˙ + ρ0θ − ∂u ∂η η˙. (4.17) 19 In thermoelastic cases where the stress and temperature are regarded as functions of the deformation and entropy, requiring non-negative dissipation in Eq. (4.17) implies that σ = ∂u ∂ , θ = 1 ρ0 ∂u ∂η . (4.18) Thus, for thermoelasticity, D = 0 and, from Eq. (4.12) and Eq. (4.13), − 1 θ q· ∇ θ ≥ 0. (4.19) Eq. (4.19) implies that, because θ > 0, heat flows in a direction of decreasing temperature. For elastic materials in which the stress is directly derivable from a strain energy potential as described, the material is said to be "hyperelastic". 1 Using Eq. (4.18), the balance of energy in Eq. (4.4) for thermoelastic materials can be expressed as ρ0θη˙ = −J ∇ ·q + ρ0r. (4.20) In other words, the only source of entropy production is from an internal heat source or flow of heat through conduction. This situation is distinguished from plasticity which allows for entropy production through dissipation. Using the conjugate relations in Eq. (4.18), the rates of stress and temperature are given by ˙σ = ∂2u ∂ ∂ :˙ + ∂2u ∂ ∂η η˙, θ˙ = 1 ρ0 ∂2u ∂η2 η˙ + 1 ρ0 ∂2u ∂η∂ :˙ . (4.21) Using the following Maxwell and Gibbs relations, 1 ρ0 ∂2u ∂η2 = ∂θ ∂η = θ cv , ∂2u ∂ ∂η = ∂σ ∂η = −ρ0θΓ, ∂2u ∂ ∂ = C, (4.22) 1The term "hyperelasticity", sometimes also referred to as "Green" elasticity, refers to the form of elasticity pioneered by Green in which the stress is derivable as the derivative of an elastic energy potential with respect to strain. In contrast, "hypoelastic" elastic models are those in which the stress is not directly derivable from a strain energy potential, see (Truesdell et al., 2004). 20 where C is the isentropic elastic stiffness, Γ is the Gr¨uneisen tensor, and cv is the specific heat at constant volume, Eq. (4.21) can be written in terms of measurable quantities ˙σ = C:˙ − ρ0θΓη˙, θ˙ = θ cv ˙ η − θΓ:˙ . (4.23) For an adiabatically loaded thermoelastic material, since η˙ = 0, the stress and temperature rates can be found from Eq. (4.23) ˙σ = C:˙ , ˙ θ = −θΓ:˙ . (4.24) For isothermal loading, since θ˙ = 0, Eq. (4.23) reduces to ˙σ = (C − ρ0θcvΓΓ) :˙ = Cθ:˙ , ˙ η = cvΓ:˙ , (4.25) where Cθ is the isothermal elastic stiffness. 4.4 Thermoplasticity 4.4.1 Conservation Laws for Thermoplastic Materials For thermoplastic materials, the internal energy depends on elastic strain2, en-tropy, and a set of internal variables that evolve with plastic loading. The stress is allowed to reach a limiting value and entropy is produced not only through heat sources and heat conduction but also through dissipation. The rate of change of internal energy can be expressed as u˙ = ∂u ∂ e :˙ e + ∂u ∂η η˙ + n k=1 ∂u ∂qk q˙k (4.26) 2Alternatively, it can be argued that the elastic energy depends on deformation through total strain, and not only the elastic part. In that case, the stress increment depends on the strain increment holding all other variables, including those that change only with dissipation, constant. In other words, the stress increment depends on the elastic strain increment, implying that the energy potential depends on the deformation through the elastic part of the strain. Throughout this document, these arguments are avoided by assuming a priori that the energy depends on the elastic strain. 21 where ˙ e is the rate of elastic strain3,qk are the internal variables that change only with dissipation, and n is the number of internal state variables for the thermoplastic material. Using Eq. (4.26), the dissipation inequality for thermoplastic materials may now be written as D = σ − ∂u ∂ e :˙ e + ρ0θ − ∂u ∂η η˙ − n k=1 ψk ˙ qk + σ:˙ p . (4.27) This form of the dissipation equality is similar to that in Wright (2002) and Rosakis et al. (2000), except that Wright and Rosakis worked with the Gibbs energy. Note that the dissipation inequality for the thermoplastic material is identical to that of the thermoelastic material in Eq. (4.17) with the addition of the last two terms on the right hand side associated with dissipation. The quantities ψk that appear in Eq. (4.27) are work conjugate to the internal state variables and are defined as ψk = ∂u ∂qk . (4.28) We now assume that the temperature gradient does not depend on the rates of strain, entropy, or internal state variables. Then the entropy inequality for a thermoplastic material can be expressed as σ = ∂u ∂ e, θ= 1 ρ0 ∂u ∂η , D = σ:˙ p − n k=1 ψkq˙k ≥ 0, F = − 1 θ q· ∇ θ ≥ 0, (4.29) and the balance of energy as ρ0θη˙ = H + D (4.30) where H = −J ∇·q+ρ0r. For codes which do not solve the heat conduction equation, such as in shock physics, H = 0. 3Here, and throughout this dissertation, superscripts e represent the elastic part of the argument and p the plastic part. 22 4.4.2 Plastic Yield In plasticity theory, for a given entropy and state of internal variables, the stress in the material is allowed to reach a limiting value, above which plastic deformation will occur. This threshold is defined by a scalar-valued function of stress, entropy, and internal state variables, known as the yield function, f. The yield criterion is expressed mathematically as f (σ, η, qk) = 0. (4.31) The yield surface is the set of all stress states satisfying this yield criterion. Elastic states correspond to f (σ, η, qk) < 0, (4.32) and plastic states, according to classical rate-independent theories of plasticity, cor-respond to f (σ, η, qk) = 0. (4.33) Viscous, or rate-dependent, theories of plasticity, which allow the stress state to lie outside of the yield surface are usually formulated as extensions of the type considered in this dissertation (cf. Duvaut and Lions, 1976). The stress and temperature evolve according to the first and second equations in Eq. (4.29), expressed in rate form as ˙σ = ∂2u ∂ e∂ e :˙ e + ∂2u ∂ e∂η η˙ + n k=1 ∂2u ∂ e∂qk q˙k, ρ0 θ˙ = ∂2u ∂η2 η˙ + ∂2u ∂η∂ e :˙ e + n k=1 ∂2u ∂η∂qk q˙k. (4.34) Using the Maxwell relations in Eq. (4.22), these relationships can be expressed as ˙σ = C:˙ e − ρ0θΓη˙ − n k=1 ∂ψk ∂ e q˙k, θ˙ = θ cv ˙ η − θΓ:˙ e − 1 ρ0 n k=1 ∂ψk ∂η q˙k. (4.35) 23 4.4.3 Plastic Flow For strain-controlled loading, the solution to the plasticity problem begins by computing the "trial" stress, in which it is assumed that the entire strain increment is elastic. For adiabatic conditions, the trial stress rate is given by ˙σ trial = C:˙ (4.36) and the trial stress can be found by first order integration of Eq. (4.36) σtrial n+1 = σn + ˙σ trial n+1Δt. (4.37) If f σtrial n+1, η, qk > 0, plastic flow occurs and corrections to the trial state state are needed to satisfy the yield criterion f σtrial n+1, η, qk ≤ 0. The assumption of additive decomposition of the strain rate into elastic and plastic parts then lead to rates of stress and entropy in Eq. (4.35) for the plastic state given by ˙σ = C: ˙ − ˙ p − ρ0θΓη˙ − n k=1 ∂ψk ∂ e q˙k, θ˙ = θ cv η˙ − θΓ: ˙ − ˙ p − 1 ρ0 n k=1 ∂ψk ∂η q˙k. (4.38) The rate of plastic strain is typically expressed in terms of its magnitude and direction ˙ p = ˙λ m, (4.39) where ˙λ is the magnitude of the rate of plastic deformation and m is its direction, typically presumed to be given by a "flow rule" (Brannon, 2007) m = ∂ϕ/∂σ ∂ϕ/∂σ , (4.40) where ϕ is the flow "potential."4The flow rule is said to be associative if m = n, where n is the yield surface normal, defined as n = ∂f/∂σ ∂f/∂σ . (4.41) 4The existence of a flow potential is questioned by Brannon and Leelavanichkul, who argue in favor of setting m via a state-dependent transformation of n (Brannon and Leelavanichkul, 2009). In either case, a key assumption is that m depends only on the material state, not the rate of change of the state. 24 Assuming each internal state variable qk changes only in response to plastic loading, they are usually expressed in terms of ˙λ as qk = hkλ˙ , (4.42) where hk is an evolution modulus corresponding to each internal state variable (cf. Brannon, 2007). Substituting Eq. (4.39) and Eq. (4.42) into Eq. (4.38) gives the nonlinear coupled evolution of the stress and temperature ˙σ = C: ˙ − ˙λ m − ρ0θΓη˙ − ˙λ n k=1 ∂ψk ∂ e hk, θ˙ = θ cv η˙ − θΓ: ˙ − ˙λ m − ˙λ n k=1 ∂ψk ∂η hk. (4.43) Using the balance of energy in Eq. (4.30), the rate of stress and temperature (assuming adiabatic conditions5 with no heat sources) in Eq. (4.43) can be expressed as ˙σ = ˙σ trial − ˙λ (C:m + Λ) , = ˙σ trial − ˙λ p, (4.44) θ˙ = θ˙trial − λ˙Pθ, (4.45) where ˙σ trial is the "trial" stress rate found by presuming the entire strain increment is elastic, p, the "return direction," is given by p = C:m+Λ, and the elastic-plastic coupling tensor Λ and Pθ are given by Λ = (σ:m)Γ + n k=1 ψkΓ + ∂ψk ∂ e hk, (4.46) Pθ = − 1 ρ0cv σ:m − θΓ:m − 1 ρ0 n k=1 ψk cv − ∂ψk ∂η hk. (4.47) Specific forms of the hk depend on the evolution equations for qk (cf. Brannon, 2007). The value of ˙λ is found by requiring that, after the onset of yield, the stress remain on the yield surface. This requirement, known as the consistency condition, 5This includes locally adiabatic conditions for which heat flow is neglected, as is the case in any code (such as a shock physics code) that does not solve the heat flow equations. 25 is represented mathematically by forcing the yield functions rate to be zero. Namely, recalling Eq. (4.33), f˙ = ∂f ∂σ :˙σ + ∂f ∂η η˙ + ∂f ∂qk q˙k = 0. (4.48) Dividing by ∂f/∂σ to normalize, gives n:˙σ = G˙λ +Ξ˙ η, (4.49) where G, and Ξ are given by G = − ∂f/∂qk ∂f/∂σ hk, Ξ = − ∂f/∂η ∂f/∂σ . (4.50) Substituting Eq. (4.47) and Eq. (4.46) into Eq. (4.49) and solving for ˙λ gives ˙λ = n:C:˙ H + n:C:m + n:Λ , (4.51) where H is the ensemble hardening modulus given by H = G − 1 ρ0θ Ξ σ:m − n k=1 ψkq˙k (4.52) It can be shown that first order integration of Eq. (4.44) leads to an updated stress of the following form (Brannon, 2007), whether or not the strain increment was partially or fully plastic σnew = σtrial − Λp, (4.53) where the scalar Λ is the magnitude of the plastic strain increment over the timestep and is determined by requiring that f σtrial − Λp = 0, which can be solved by standard numerical methods. Equation (4.53) has a convenient physical interpretation: the updated stress σnew is the oblique projection of σtrial on to the yield surface defined by f (σ, η, qk) = 0, p is the direction of the projection, and Λ is its magnitude, as depicted in Figure 4.1. Thus completes the system of equations for the stress and temperature evolution in a thermoplastic material. For further details consult (Brannon, 2007) and (Bran-non and Leelavanichkul, 2009). Discussion of how to solve the resulting system of equations will be postponed until later in this dissertation. 26 4.4.4 Temperature Evolution Due to Dissipation For the case of adiabatic loading of a thermoplastic material, common in applica-tions which involve high rates of deformation, and in the absence of any heat sources, it was shown that the temperature evolves according to ˙ θ = −θΓ:˙ e + 1 ρ0cv σ:˙ p + 1 ρ0 n k=0 1 cv ψk − ∂ψk ∂η q˙k. (4.54) In the literature, the temperature evolution in a thermoplastic material is often expressed as θ˙ = χ ρ0cv σ:˙ p , (4.55) where χ, known as the Taylor-Quinney coefficient, is commonly assigned a constant value between 0.7 and 1.0, a reflection of the experimental evidence that not all plastic work is converted to heat (Taylor and Quinney, 1934). Writing Eq. (4.54) in the above form we get the following for the Taylor-Quinney coefficient χ = 1− 1 σ:˙ p θΓ:˙ e + 1 ρ0 n k=0 ∂ψk ∂η − 1 cv ψk q˙k . (4.56) Clearly, χ evolves with both θ and qk and, thus, should not be assumed to be constant. Doing so will undoubtedly lead to errors in model predictions which can lead to the assignment of erroneous values of model parameters to match experimental data. Rosakis et al. (2000) derived a similar expression for χ and developed a series of experiments to measure χ indirectly. Depending on the specific form of the functions ψk and q˙k, Eq. (4.54) has the desired quality that only a portion of the plastic work is converted to heat without resorting to introducing empirical parameters such as χ. A specific example is presented later in Section 8.2. 27 Elastic Domain Λp σtrial σnew ˙σ trial Δt f (σ, ξ) = 0 σold ˙σ Δt Figure 4.1. Oblique return projection of the trial stress state on to the yield function isosurface. PART II THERMODYNAMIC REQUIREMENT OF RDIA CHAPTER 5 THERMODYNAMIC REQUIREMENT OF RDIA 5.1 Elasticity in Isotropic Media If a material is capable of elastic behavior, a necessary condition of thermodynamic admissibility is that the stress be derivable from a strain energy potential (cf. Truesdell et al., 2004; Malvern, 1969) σ = ∂u ( , η, qk) ∂ η,qk , (5.1) where u is the internal energy per unit reference volume, is the strain measure1 work conjugate to the stress tensor σ, η is the specific entropy, and qk are internal state variables that change only with dissipation. For an isotropic material, the stress remains invariant under rotation, requiring that σ = w( , η) = w(Q∗ , η) = Q∗w( , η), (5.2) where w = ∂u/∂ . By Definition 3.3.1, the tensor valued function w may be expressed as w( , η) = w1δ + w2γ + w3h, (5.3) where the wi are scalar functions of the mechanics invariants of , γ is the deviatoric part of , and h is the deviatoric part of γ·γ. As can be readily verified by back substitution, Eq. (5.3) implies that the internal strain energy of an isotropic material 1The definition of the strain tensor is intentionally left ambiguous so that the following results can be considered general in nature. In specific applications, of course, the strain tensor would take on a specific meaning and the corresponding stress tensor would necessarily by its work conjugate. For example, if is the Lagrange strain, then its work conjugate stress would be the 2nd Piola-Kirchhoff stress tensor. 30 is a function of the mechanics invariants2 of . Mathematically, u = u(J 1 , J 2 , J 3 , η) and the fourth-order elastic stiffness tensor of an isotropic material is found by C = ∂2u J 1 , J 2 , J 3 , η ∂ ∂ η . (5.4) Performing the indicated differentiation, the most general form of the elastic stiffness tensor of an isotropic material is C = ∂2u ∂ ∂ = 3u11 Eiso + u2√5Esd + u12 (γδ + δγ) + u22γγ + u13 hδ + δh + u23 hγ + γh + u33hh + u3 L hδ + δh , (5.5) where ui and uij denote ui = ∂u/∂J i and uij = ∂2u/∂J i ∂J j respectively, and not first-order and second-order tensors. The operator L is given in indicial form by L(Γijkl) = 1 4 (Γikjl +Γkijl +Γiklj +Γkilj) . (5.6) Since γδ, δγ, γγ, hδ, δh, hγ, γh, and hh are not expressible as linear combinations of the basis tensors Eiso and Esd in Eq. (3.7), C is not expressible as in Eq. (3.6) and is anisotropic except in the following cases: u3 = uij = 0 (except i = j = 1), = αδ, and = γ = h = 0. Thus, thermodynamics requires that the elastic stiffness of isotropic materials develop RDIA in response to deformation. 5.2 Elastic Stiffness Isotropy If instead the elastic stiffness is presumed isotropic, by Definition 3.2.2, it is expressible as Ciso = 3κEiso + 2μ√5Esd, (5.7) where κ and μ are the tangent bulk and shear moduli, respectively, modulo a factor of J = detF (where F is the deformation gradient) depending on the definition of . 2Alternatively, any other triplet of independent invariants of the deformation measure can be used. 31 5.3 Elastic Stiffness Isotropy It is natural to consider the conditions under which the elastic stiffness of an isotropic material is itself isotropic. Comparing tensor coefficients of Eq. (5.5) and Eq. (5.7) leads to the following necessary conditions for elastic stiffness isotropy u11 = κ, (5.8a) u2 = 2μ, (5.8b) u12 = u22 = u13 = u23 = u33 = u3 = 0. (5.8c) Substituting Eq. (5.8b) into Eq. (5.8c) implies that ∂μ ∂J 1 = ∂μ ∂J 2 = ∂μ ∂J 3 = 0. (5.9) A necessary and sufficient condition for Eq. (5.9) to be satisfied is that μ be indepen-dent of the deformation measure. Eq. (5.8a) and Eq. (5.8c) also imply that κ = κ(J 1 ). Thus, an isotropic elastic material in a distorted state will have an isotropic elastic stiffness if, and only if, its shear modulus is constant and bulk modulus varies at most with J 1 . Hueckel et al. (1992) arrived at a similar conclusion with respect to the secant bulk and shear moduli and Zytynski et al. (1978) showed that the use of a pressure-dependent shear modulus in conjunction with elastic stiffness isotropy leads to a model which is nonconservative in closed cycles of elastic stress. Zytynski et al., however, wrongly asserted that the alternatives were to either use a constant shear modulus or abandon the notion of elastic behavior, neglecting the possibility of satisfying energy conservation through RDIA. For elastic-plastic materials, similar analysis implies that μ is independent of the deformation measure but can still depend on η and qk. This may be a means of accounting for apparent dependence of μ on pressure through elastic-plastic coupling, though this is not considered further in this dissertation. 32 5.4 First Order Approximation of C In this and subsequent sections, the first-order approximation of C in Eq. (5.5) C = 3u11 Eiso + u2√5Esd + u12 (γδ + δγ) (5.10) is adopted. Comparing tensor coefficients of Eq. (5.10) to Eq. (5.5), the derivatives of the strain energy can be given in terms of the familiar bulk and shear modulus u11 = κ J 1 (5.11a) u2 = 2μ J 1 (5.11b) u12 = 2 dμ J 1 dJ 1 = 2μ . (5.11c) where the μ indicates differentiation with respect to J 1 . With respect to the orthonormal basis in R6, as described in Eq. (3.1), C is C = H + 2 3α1 2α2 2α2 2μ (5.12) where H = ⎛ ⎝ κ + 4 3μ κ− 2 3μ κ− 2 3μ κ + 4 3μ κ− 2 3μ sym. κ + 4 3μ ⎞ ⎠, μ = ⎛ ⎝ μ 0 0 μ 0 sym. μ ⎞ ⎠, (5.13) α1 = μ ⎛ ⎝ 2 (2 1 − 2 − 3) 1 + 2 − 2 3 1 − 2 2 + 3 −2 ( 1 − 2 2 + 3) −2 1 + 2 + 3 sym. −2 ( 1 + 2 − 2 3) ⎞ ⎠, (5.14) and α2 = μ ⎛ ⎝ 2 4 2 5 2 6 2 5 2 6 sym. 2 6 ⎞ ⎠. (5.15) where the i are the components of = ( 1, 2, 3,√2 4,√2 5,√2 6)T. In matrix form, the anisotropy of C, induced by deformation, is evident in the α1 and α2 sub-matrices, i.e., C is isotropic only if μ = 0 and/or α1 = α2 = 0. 5.5 Elastic Compliance In numerical constitutive models, the fourth-order elastic compliance tensor S = C−1 is needed to evaluate the elastic strain increment when the stress increment 33 is known. Successive application of the Sherman-Morrison formula (cf. Golub and Van Loan, 1996) for the inversion of a rank-one modification of a tensor allows the rank-two modification in Eq. (5.10) to be inverted to give S = 1 3κ Eiso + √5 2μ Esd − 1 φ 3κμμ (γδ + δγ) − (μ )2(9κγγ + 4μJ 2 δδ) , (5.16) where φ = 9κμ κμ − 4J 2 (μ )2 . (5.17) 5.6 Summary In this chapter the distinction between isotropic functions and isotropic tensors was emphasized. An isotropic material is one in which its internal energy is an isotropic function of the invariants of the strain tensor. Twice differentiating this isotropic energy function results in thermodynamically required deformation-induced anisotropy in the fourth-order elastic stiffness. The stiffness will be anisotropic except in the case that the bulk modulus varies at most with J 1 and the shear modulus is constant. Considering laboratory evidence that μ is not constant, even in common engineering materials (Guinan and Steinberg, 1975; Hayes et al., 1999; Duffy et al., 1999), the assumption of isotropy of the elastic stiffness should be abandoned to allow revised analysis of material characterization data in a thermodynamically consistent manner. In other words, if using an isotropic stiffness leads to a pressure-dependent shear modulus, then the data must be reanalyzed allowing for RDIA, or irreversible changes in pressure associated with dissipation rather than elastic volume change. CHAPTER 6 NUMERICAL IMPLEMENTATION OF RDIA 6.1 Introduction In this chapter, algorithms for implementing RDIA in the stress response of pressure-dependent elastic and elastic-plastic materials are given. In each algorithm, because of common use in the engineering community, the unrotated Cauchy stress σ and unrotated symmetric part of the velocity gradient d are taken as input1. Accordingly, σ is first transformed to the work conjugate of d, the Kirchhoff stress τ 2. Output is the unrotated Cauchy stress at the end of the step. Any term of the form • are those that account for RDIA and the other terms are those which are probably already coded in a thermodynamically inadmissible model. 6.2 Elastic Stress Update Implementing RDIA in an existing elastic stress-update routine involves the ad-dition of two easily computed terms, shown set apart from terms in a typical elastic stress-update by • in Algorithm 6.1. 1Here, "unrotated" refers to the now common practice of defining a reference tensor A from a spatial tensor A by A = R∗A where R is the rotation from the polar decomposition of the deformation gradient. 2Generally, the Kirchhoff stress τ is work conjugate to d only if the principal directions of the reference stretch are stationary. For problems in which the principal directions of the reference stretch do change, the difference between the d and the true rate of strain is small (Brannon, 2009); thus, it will be assumed that τ is work conjugate to d, as is commonly done in numerical constitutive modeling. 35 6.3 Elastic-Plastic Stress Update For isotropic elastic-plastic materials the strain in the previous sections might be only the elastic strain. More generally, additional terms due to RDIA, such as those terms set apart in Algorithm 6.1 by • , would be present in any term operated on by the fourth order elastic stiffness. 6.3.1 Elastic-Plastic Stress Update Referring to Section 4.4, given a strain increment d , the updated stress is found by presuming that the entire strain increment is elastic and computing the corresponding "trial" stress σtrial = σold + C:d . (6.1) For elastic loading f σtrial ≤ 0 and σnew = σtrial. For inelastic loading f σtrial > 0 and the actual stress increment is found by projecting the trial stress increment onto the yield surface. Assuming d = d e + d p, first order integration of the trial stress increment leads to an updated stress of the form σnew = σtrial − Λp, (6.2) where the second-order tensor p is given by p = C:m − (σ:m)Γ, (6.3) The magnitude of the plastic strain increment Λ is found by requiring that f σtrial − Λp = 0, (6.4) Referring to Eq. (6.1) - Eq. (6.3), only two terms, σtrial and p, are found from direct operation of the elastic stiffness on other tensors and, therefore, are the only terms which require explicit inclusion of the additional RDIA terms. Specifically, the trial stress and return direction, accounting for RDIA are given by σtrial = σold + κ˙ vδ + 2μ˙γ + 2μ ˙ vγe + γe: ˙ δ (6.5) p = κ trmδ + 2μm + 2μ (trmγe + γe:mδ) (6.6) 36 6.3.2 Elastic Strain Update With the stress increment known, the updated elastic strain e is found by enew = eold + S:˙σ (6.7) where the elastic compliance S is given by Eq. (5.16). In the coding, S is not evaluated explicitly. Instead, the action of S on ˙σ is applied in the same manner as the action of C on ˙ . 6.3.3 Numerical Solution Method Rather than explicitly solving for Λ in Eq. (6.4) numerically and substituting its value in to Eq. (6.2) to obtain the updated stress, the following modification to Newton's method is adopted which works by modifying the components of σtrial directly, iterating until σ − σtrial is smaller than an acceptable tolerance. σ0 = σtrial σn+1 = σn − f (σn) (∂f/∂σ):p p (6.8) The coefficient of p in Eq. (6.8) is found by taking the first order expansion of f σtrial − Λp about Λ = 0 f σtrial − Λp ≈ f σtrial − Λ ∂f ∂σ :p = 0 (6.9) solving for Λ leads to Λ = f σtrial (∂f/∂σ):p (6.10) which is the coefficient of the p in Eq. (6.8). The algorithm described above is classical in the sense that it returns the stress to the yield surface along a path perpendicular to nonzero level sets of the yield function. Accordingly, as pointed out by Brannon and Leelavanichkul (2009) who investigated consequences of nonuniqueness of yield functions, the above algorithm requires small time steps to remain accurate. Larger time steps can be accommodated by using a second stage correction to the result of the above algorithm (Brannon and Leelavanichkul, 2009). 37 6.3.4 Stress Update Algorithm For strain increment driven models, Algorithm 6.2 implements RDIA in the elastic stiffness and compliance for isotropic elastic-plastic materials. As can be seen in Steps 1, 3, and 5 of Algorithm 6.2, endowing the effects of RDIA on an existing elastic-plastic model is a matter of adding very few easily evaluated terms. The only terms which are not readily available in most finite element host codes are the elastic strains, which can be easily carried as internal state variables to the model for a moderate additional computational cost. 38 Algorithm 6.1. Elastic Kirchhoff stress-update algorithm. New terms involving RDIA, not expected to be included in an existing stress-update algorithm, are set apart explicitly by • in Step 1. Input Unrotated Cauchy stress σold, unrotated (approximate) strain increment dε = dΔt, timestep Δt, Jacobian Jold, elastic strain εold. Output Cauchy stress σnew Step 1 Compute updated Kirchhoff stress τ = Joldσold + κ(εv)dεvδ + 2μ(εv)dγ + 2μ (dεvγold + γold:dεδ) Step 2 Update energy, temperature, Jacobian, and all other state variables to end of step Step 3 Convert updated Kirchhoff stress to Cauchy stress σ = 1 J τ 39 Algorithm 6.2. Elastic-plastic Kirchhoff stress-update algorithm. New terms involving RDIA, not expected to be included in an existing stress-update algorithm, are set apart explicitly by • in Steps 1, 3, and 5. Input Unrotated Cauchy stress σold, unrotated (approximate) strain increment dε = dΔt, timestep Δt, Jacobian Jold, elastic strain εe old. Output Cauchy stress σnew Step 1 Compute trial Kirchhoff stress τ trial = Joldσold + κ(εv)dεvδ + 2μ(εv)dγ + 2μ (dεvγe old + γe old:dεδ) Step 2 Check yield condition if f τ trial ≤ 0 then τ = τ trial, dεe = dε, dεp = 0 go to Step 6 endif Step 3 Compute yield normal, flow direction, and return direction n = ∂f/∂τ ∂f/∂τ , m = ∂ϕ/∂τ ∂ϕ/∂τ , p = κ(εv) trmδ + 2μ(εv) devm− Joldσold:mΓ + 2μ (trmγe old + γe old:mδ) Step 4 Improve estimates for stress and check convergence τ = τ trial − f τ trial (∂f/∂τ ):p p if τ − τ trial > TOL then τ trial = τ go to Step 3 endif Step 5 Update elastic and plastic strain κ = 4 Jε 2μ κμ , ζ = 3 2 κμ(1 − κ) μ α1 = 1 ζ 2 3 Jε 2 tr dτ κ − γe old:dτ 2μ , α2 = 1 ζ 3 2 Jε 2 γe old:dτ μ − tr dτ 2μ dεe = 1 3κ tr dτδ + 1 2μ dev dτ + α1δ + α2γe old dεp = dε − dεe Step 6 Update energy, temperature, Jacobian, and all other state variables to end of step Step 7 Convert updated Kirchhoff stress to Cauchy stress σ = 1 J τ PART III EFFECTS OF RDIA ON ISOTROPIC MATERIALS CHAPTER 7 QUANTIFICATION OF ANISOTROPY 7.1 Introduction It is clear from the previous chapters, Eq. (5.5) in particular, that the elastic stiffness of isotropic elastic and elastic-plastic materials must develop RDIA in its integrity basis in response to distortion. To objectively assess the effects of the resultant anisotropy on classical elastic-plastic materials, it is necessary to quantify the degree of the resultant anisotropy. The notion of quantifying the degree of anisotropy of a tensor by a single scalar measure has previously been considered in a number of contexts and disciplines (Nye, 1957; Backus, 1970; Rychlewski, 1984; Pierpaoli and Basser, 1996). Though the formulations differ in their details, they share the common objective of attempting to quantify, in a meaningful way, the degree to which the tensor is not isotropic. One of the first measures of tensor anisotropy ψ was the ratio of the magnitude of a tensor and its deviator (Nye, 1957) ψ (A) = A A , (7.1) where A is the deviatoric part of a tensor of arbitrary order A. This form of ψ has been used, for example, in the study of second-order tensors in crystal optics where the degree of anisotropy has important implications on the behavior of light in anisotropic crystals (Fedorov, 1968). Though attractive for its intuitiveness, this method of quantifying the anisotropy of a tensor has been criticized for being too simplistic (Rychlewski, 1984). For example, as formulated, ψ is incapable of differentiating between tensors whose deviators are equal in magnitude but differ in Lode angle θ, an important distinction when considering the behavior of pressure-sensitive materials. 42 Recognizing the fundamental property of anisotropic tensors that their compo-nents change with rotation, Rychlewski (1984) proposed as the measure of anisotropy of a second-order tensor the ratio of its orbit to its magnitude ψ∗(A) = d(A) 2 A , (7.2) where the orbit of the tensor d(A) is the maximum distance between tensors attain-able by rotation of A d(A) = max R1,R2∈R R1∗A − R2∗A . (7.3) Here, R is the set of rotation tensors and ∗ represents the Rayleigh product. Using this formulation, ψ is reducible to ψ∗(A) = A A sin π 3 + θ = ψ sin π 3 + θ , (7.4) where the factor sin(π/3 + θ) is referred to as the degree of variability of the tensor A and θ is the "modular angle" of A, defined as cos 3θ = 3√6JA 3 2 A 2 . (7.5) Zhang (1990) applied Eq. (7.4) to the study of the second-order stress tensor in linear and nonlinear elastic and inelastic materials, developing bounds on the degree of anisotropy of the stress in these materials. More interesting is the quantification of anisotropy in fourth-order tensors, such as the tangent elastic stiffness and compliance tensors of elastic and inelastic materials. Due to multiple symmetry planes and a large number of independent components in generally anisotropic fourth-order tensors, developing compact closed form expres-sions, such as in Eq. (7.4) is generally not possible, except for limiting cases such as cubic or transversely isotropic symmetry (Zhang and Rychlewski, 1990; Rychlewski, 1995). Other methods of visualizing the degree of anisotropy in fourth-order and higher-order tensors have been established and used in the study of the fourth-order elastic stiffness tensor of linear elasticity for uniaxial loading (B¨ohlke and Br¨uggemann, 2001), though a single scalar measure of the degree of anisotropy was not used. 43 Here, a scalar measure of tensor anisotropy is defined for tensors of arbitrary order. An important objective of quantifying the degree of tensor anisotropy is to determine if the degree of RDIA in isotropic materials is significant and, if RDIA is significant to determine the necessary conditions for RDIA to become significant. 7.2 Scalar Measure of RDIA Commonly, ψ is defined by the ratio of the norm of the difference between a tensor A and its isotropic part and the norm of A (Fedorov, 1968) ψ (A) = A−Aiso A , (7.6) where Aiso, the isotropic part of A, is given by Aiso = P•A, (7.7) and P is the tensor projector which projects A on to the hyper-surface spanned by the isotropic basis for nth-order tensors and • represents the appropriate order tensor contraction. In this form, ψ has previously been used in a variety of applications, for example, in the investigation of crystal optics where the degree of anisotropy has important implications on the behavior of light in anisotropic crystals (Fedorov, 1968), second-order tensors in mechanics (Rychlewski, 1984; Zhang, 1990), and finding the isotropic tensor closest to an anisotropic symmetry (Cavallini, 1999; Norris, 2006b; Moakher and Norris, 2006). Rather than adopt Eq. (7.6), Brannon (2009) makes the following observation: given a tensor A, Aiso is the projection of A onto the hyper-surface spanned by the unit basis tensors of the isotropic part of A. For fourth-order tensors, the basis tensors are given by Eiso and √5Esd, and a geometric description of the projection is shown in Figure 7.1. The angle θ between A and Aiso, defined analogously to the angle between two vectors, is given by cos θ = A::Aiso A Aiso = Aiso A . (7.8) 44 Rearranging and normalizing, the scalar measure of anisotropy is defined as ψ (A) = 2 π cos −1 Aiso A , (7.9) 7.2.1 Properties of Anisotropy Measure As given in Eq. (7.9), ψ has the following properties: ψ Aiso = 0, ψ(A ) = 1, ψ (A) ≥ 0, ψ (αA) = αψ (A), α∈ R, ψ (Q∗A) = ψ (A) . Because ψ is based on the Euclidean norm, ψ is not invariant under inversion (Norris, 2006a), i.e., ψ A−1 = ψ (A) . (7.10) However, several of the previous alternative measures suffer from the same problem. If invariance under inversion is required, the subsequently reported work can be re-run using ψ†(A) = 1 2 ψ (A) + ψ A−1 (7.11) which is invariant under inversion. Despite this drawback, Eq. (7.9) is adopted due to its intuitive geometric interpretation and computational ease. 7.3 Evaluation of Degree of Anisotropy in Isotropic Materials For fourth-order tensors, P = P, where P = Eiso ijkl Eiso pqrs + Esd ijkl Esd pqrs eiejekelepeqeres. (7.12) Fortunately, because of the symmetries of P, one need not actually assemble P and perform the projection operation as indicated in Eq. (7.7). Instead, given the R6 components of A A = A1 A3 A3 A2 , (7.13) 45 the isotropic part is Aiso = λ1 1 0 0 0 + λ2 δ 0 0 δ (7.14) Here, 0 and 1 are the 3 × 3 submatrices with components 0 and 1, respectively, and λ1 = 2a − b 15 , λ2 = 3b − a 15 , (7.15) where a = i=1 3 j=1 3Aij 1, b= i=1 3 Aii 1 + Aii 2 . (7.16) Applying Eq. (7.14) - Eq. (7.16) to the reduced elastic stiffness of Eq. (5.10), the isotropic part C is given by Ciso in Eq. (5.7), and ψ (C) is reducible to ψ (C) = 2 π cos −1 9κ2 + 5(2μ)2 9κ2 + 5(2μ)2 + 12(2μ )2J 2 (7.17) 7.4 Degree of Anisotropy in Anisotropic Materials For perspective in the subsequent sections in which the degree of RDIA in isotropic elastic-plastic materials will be characterized, the degree of anisotropy ψ of the elastic stiffness of several inherently anisotropic materials is shown plotted against Ciso / C in Figure 7.2 for the materials summarized in Table 7.1. In Table 7.1, "Sym. Type" is the symmetry type of the material, N is the number of symmetry axes, and the Cij refer to the coefficients of the elastic stiffness tensor in R6. For further explanation and sources of material data see B¨ohlke and Br¨uggemann (2001) and the references therein. 7.5 Degree of RDIA in Isotropic Materials In contrast to the materials depicted in Figure 7.2, whose degree of anisotropy is inherent and constant, the degree of anisotropy in isotropic elastic materials varies with deformation. As will be shown, the degree of RDIA in classical elastic-plastic materials generally depends on two competing constitutive factors: pressure depen-dence of the shear modulus and pressure-dependence of the yield function. These two features will be explored in this section. 46 7.5.1 Constitutive Model Isotropic engineering materials will be analyzed using Algorithm 6.2 and the following constitutive features/assumptions. 7.5.1.1 Yield Function For each material, the yield function in Algorithm 6.2 is given by the pressure-dependent form f (σ) = Jσ 2 − a1 + a3ea2Jσ 1 + a4Jσ 1 + a5 γp eq m, (7.18) where a1 - a5 and m are material parameters, and γp eq is the equivalent plastic strain. As defined, f is capable of reducing to several common plasticity models, e.g., choosing a2 - a4 equal to zero, J2 plasticity with isotropic strain hardening is recovered, choosing a2, a3, and a5 equal to zero, linear Drucker-Prager plasticity is recovered, or, if a5 is set equal to zero, the yield function employed by the Sandia Geomodel (Fossum and Brannon, 2004, 2006; Brannon et al., 2009) is recovered. 7.5.1.2 Elastic Moduli Two bulk modulus functions will be used in subsequent analysis: Mie-Gr¨uneisen equation of state, common in high pressure applications such as shock loading, κ = 1 ρ 1 − γ (ρ0 − ρ) 2ρ0 − dp d(1/ρ) H + γ 2 pH, (7.19) where ρ is the density, pH is the pressure on the Hugoniot, and γ is the Gr¨uneisen parameter, and the first order expansion of κ about J 1 = 0 κ = κ0 − κ0 dκ dp 0 J 1 . (7.20) The shear modulus is taken as its first order expansion about J 1 = 0 μ = μ0 − κ0 dμ dp 0 J 1 . (7.21) In Eq. (7.20) and Eq. (7.21), since elastic moduli are typically measured as functions of pressure, use has been made of the chain rule to cast first order bulk and shear moduli coefficients in terms of pressure derivatives. 47 7.5.2 Degree of RDIA in Elastic Materials Before quantifying the degree of anisotropy in elastic-plastic materials, the degree of anisotropy for an elastic material subjected to uniaxial strain compression is shown in Figure 7.3. Parameters for the constitutive model are for Titanium Diboride (TiB2) κ0 = 237 GPa, μ0 = 246 GPa, dμ/dp|0 = 9, dκ/dp|0 = 2, ρ0 = 4520, γ = 1.71. The Mie-Gr¨uneisen equation of state was used in the evaluation of the bulk modulus and yield was suppressed by taking a4 to be an arbitrarily large number. The data in Figure 7.3 demonstrate that, for elastic loading, the magnitude of anisotropy induced by deformation in isotropic materials is capable of attaining and even exceeding the magnitude of anisotropy in intrinsically anisotropic materials at finite strains. 7.5.3 Degree of RDIA in Elastic-Plastic Materials Of course, in elastic-plastic materials, due to the limit on the magnitude of realizable elastic deviatoric strain imposed by the yield criterion, the magnitude of RDIA will also be limited. As an illustration, the degree of RDIA for the same material as Section 7.5.2 subjected to the same uniaxial strain compression, is shown in Figure 7.4 allowing for yield by taking a1 = 1800 MPa and a4 = 0.095. As illustrated in Figure 7.4, allowing for yield limits the degree of RDIA in this material quite severely. Only at very large strains does the degree of anisotropy even exceed 5%. 7.5.3.1 Factors Affecting the Degree of RDIA Considering the term involving RDIA in the elastic stiffness in Eq. (5.10), CRDIA = 2μ (γeδ + δγe) , (7.22) two key factors in how large the magnitude of RDIA becomes during loading are: pressure-dependence of the shear modulus, and, implicitly, the pressure dependence of yield which allows the realizable elastic deviatoric strain to increase with pressure. It turns out that under some circumstances these two are competing factors, as will be shown. In the following paragraphs, each factor will be considered separately using the same model parameters for TiB2 as used in Figure 7.4. 48 7.5.3.1.1 Pressure-dependence of the shear modulus. In Figure 7.5, the magnitude of RDIA is shown plotted versus axial strain for uniaxial strain deformation for varying dμ/dp. All other model parameters are the same as for the material in Figure 7.4. Referring to Figure 7.5, increasing dependence of the shear modulus on pressure through increasing values of dμ/dp does not necessarily lead to an increase in the magnitude of RDIA at all strains. This behavior is easily explained by observing that stronger dependence of the shear modulus on pressure results in an increase in deviatoric strain at a given pressure. Thus, unless the pressure-dependence of strength is also large, the elastic-deviatoric strain will decrease with increasing pressure due to the yield criterion, causing a decrease in the magnitude of CRDIA. Depending on relative competing effects of pressure-dependence of strength and shear modulus, even materials with moderate pressure-dependence of shear modulus can exhibit large degrees of RDIA at large strains. 7.5.3.1.2 Pressure-dependence of the yield. In Figure 7.6, the magnitude of RDIA is shown plotted versus axial strain for uniaxial strain deformation for increasing levels of pressure dependence of strength achieved by varying a4. All other model parameters are the same as for the material in Figure 7.4. Based on the observation of Lee et al. (2004) that the slope of the yield function in stress space never exceeded the slope of the stress trajectory for uniaxial strain loading, the upper limit of a4 = 0.399 was chosen as this limiting value. Not surprisingly, since the size of the elastic domain increases monotonically with a4, the degree of RDIA also increases at all strains monotonically with a4. 7.6 Magnitude of RDIA in Engineering Materials It has been shown that, for pressure-sensitive materials, the RDIA of the elastic stiffness can be comparable in magnitude to the degree of anisotropy in inherently anisotropic materials in the finite strain regime for elastic loading, though its mag-nitude was severely restricted when allowing for yield. The degree of RDIA ψ (C) will now be shown in Figures 7.7-7.10 for limestone, 2024 aluminum, TiB2, and 49 99.5% alumina whose properties are outlined in Table 7.2 for uniaxial and triaxial = ( , /10, /10, 0, 0, 0)T strain compression. Strains are taken as positive in com-pression. Though the magnitude of RDIA is characterized for the specific materials outlined, they were chosen to be representative of the larger classes of geologic, ceramic, metallic, and powdered metal materials, respectively. As seen in Figures 7.7 - 7.10, the degree of anisotropy is largely independent of the deformation paths and bulk modulus functions investigated. The most important factors are the dependence of the yield strength on pressure. 7.7 Discussion and Summary In initially isotropic materials whose shear modulus and strength vary with pres-sure, the degree of RDIA increases with axial strain in uniaxial and triaxial strain deformation. The degree of RDIA is not sensitive to the loading paths investigated, being nearly the same in uniaxial and triaxial compression. For the metallic, geologic, ceramic, and powdered metal materials investigated, the following general observa-tions can be made: • For metallic materials, even with the inclusion of slight pressure-dependence of yield as reported by Richmond and Spitzig (1980), the degree of RDIA is, even at finite strains, negligible. This observation makes sense in light of the observation that the ratio of strength and shear modulus is constant (Hua et al., 2002) in metals, which limits the elastic-deviatoric strain from achieving magnitudes which allow for the development of RDIA in the material. • For geologic materials, in which the strength and shear modulus are both moderately dependent on pressure, the magnitude of RDIA exceeds 3% only in large deformation, finite strain, regimes. • For ceramic materials, in which the strength and shear modulus can both be strongly dependent on pressure, the magnitude of RDIA exceeds 2% at moderate strains. 50 • For the powdered metals in which the strength and shear modulus are extremely pressure-dependent, the degree of RDIA is significant even at small strains. RDIA is only one source of deformation induced anisotropy. For ceramics and rocks, which fail through growth of oriented cracks, inelastic deformation-induced anisotropy is expected to be severe. Nevertheless, the inclusion of RDIA in a com-putational framework is only a matter of adding a few extra lines of easily evaluated code to existing models, as shown in Algorithms 6.1 and 6.2 and should, therefore, be included if thermodynamic admissibility is a desired attribute of the computational model. 51 θ A A Esd αAEiso Eiso √5βAEsd Aiso Figure 7.1. Geometric interpretation of decomposition of fourth-order tensor into isotropic and nonisotropic parts, showing the angle θ between the tensor A and the isotropic hyper-plane. 52 Table 7.1. Anisotropic elastic stiffness coefficients of several inherently anisotropic materials. Material Sym. Type N C11 C22 C33 C44 C55 C66 C12 C13 C14 C15 C23 C25 C35 C46 Hornblende Monoclinic 13 130.1 187.7 198.4 61.1 38.7 45 61.1 59.5 9.5 61.4 -6.9 -40.6 -0.9 Uranium Orthorhombic 9 214.8 198.6 267.1 124.4 73.4 44.3 46.5 21.8 107.6 Quartz Trigonal 6 86.8 105.75 58.2 7.04 11.91 -18.04 Zircon Tetragonal 6 73.5 46 13.8 16 9 -5.4 Titanium Hexagonal 5 123.1 152.9 30.7 99.6 68.8 Copper Cubic 3 168 75.4 121.4 53 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 ||Ciso|| / ||C|| ψ(C) 0.14 0.26 0.92 0.98 Enlargement Zircon Quartz HUorarnnbiulemnde Titanium Copper Figure 7.2. Measure of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic symmetry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (monoclinic symmetry), and Copper (cubic symmetry). Material data from B¨ohlke and Br¨uggemann (2001). 54 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 ||Ciso|| / ||C|| ψ(C) εA=-0 εA=0.05 εA=0.116 εA=0.23 εA=0.38 εA=1 0.12 0.34 0.912 0.982 Enlargement εA=0.23 εA=0.38 εA=0.116 Figure 7.3. Variation of the measure of anisotropy for elastic uniaxial strain defor-mation ×+. Also shown are the measures of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic sym-metry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (monoclinic symmetry), and Copper (cubic symmetry). 55 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 ||Ciso|| / ||C|| ψ(C) εA= -0.00 εA= 0.15 0 0.06 0.9965 1 Enlargement εA= 1.00 εA= 0.05 εA= 0.15 εA= -0.00 Figure 7.4. Variation of the measure of anisotropy for uniaxial strain deformation for the same material as in Figure 7.3 but allowing for plastic yielding ×+. Also shown are the measures of anisotropy for Zircon 22251 (tetragonal symmetry), Quartz 62894 (trigonal symmetry), Uranium ◦ (orthorhombic symmetry), Titanium 52743 • (hexagonal symmetry), Hornblende 42420 (monoclinic symmetry), and Copper (cubic symmetry). 56 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1 εA ψ(C) μ' = 50 μ' = 20 μ' = 10 μ' = 2 μ' = 0 Figure 7.5. Variation of the measure of anisotropy with μ for uniaxial strain deformation plotted against axial strain for dμ/dp = 50 −·−·−·−·−, dμ/dp = 20 ············, dμ/dp = 10 −−−−−−, dμ/dp = 2 ----, and dμ/dp = 0.0 ----. For large values of μ , the measure of anisotropy decreases due to the limiting value of shear strain. 57 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 εA ψ(C) a4 = .399 a4 = .2 a4 = .1 a4 = .05 a4 = 0.0 Figure 7.6. Variation of the measure of anisotropy with a4 for uniaxial strain deformation plotted against axial strain for a4 = 0.399 −·−·−·−·−, a4 = 0.2 ············, a4 = 0.1 −−−−−−, a4 = 0.05 ----, and a4 = 0.0 ----. Except for the largest values of a4, the measure of anisotropy never exceeds 5%. 58 Table 7.2. Material properties representative of ceramic, geologic, and metallic materials used in simulations shown in Figures 7.7 - 7.10. Property Material TiB2 1 Limestone2 2024 Al3 99.5% alumina4 κ0 (GPa) 237.0 24.0 76.0 1.549 μ0 (GPa) 246.0 15.1 28.6 2.722 ∂μ/∂p|0 9.0 2 1.8 208.5 ∂κ/∂p|0 2.0 0 4.75 171.6 a1 (MPa) 1800 71.9 256.0 4.2 a2 (1/MPa) 0.0 3180.0 0.0 0.0 a3 (MPa) 0.0 70.1 0.0 0.0 a4 0.095 0.171 0.0017 0.1813 a5 (MPa) 0.0 0.0 426 0.0 m NA NA 0.34 NA ρ0 (Kg/m3) 4520 2300 2770 3560 γ 1.71 1.0 2.0 NA 1 1Data from Grady (1991) 2Data from Brannon et al. (2009) 3Data from Johnson and Cook (1985); Steinberg (1996) 4Data from Zeuch et al. (2001) 59 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0 0.2 0.4 0.6 0.8 1 εA ψ(C) Figure 7.7. Degree of anisotropy in elastic-plastic materials in uniaxial strain compression for linearly pressure-dependent elastic moduli for 99.5 % alumina ············, TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. Note that, with the exception of the powdered metal (alumina), the degree of recoverable anisotropy never exceeds 6%. In the small to medium strain regime, the degree of recoverable anisotropy does not exceed 2%. 60 0 0.03 0.06 0.09 0 0.2 0.4 0.6 0.8 1 εA ψ(C) Figure 7.8. Degree of recoverable anisotropy in elastic-plastic materials in uniaxial strain compression for elastic moduli given by Mie Gr¨uneisen equation of state for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. The degree of recoverable anisotropy never exceeds 4% in the small to medium strain regime for any of the materials. 61 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0 0.2 0.4 0.6 0.8 1 εA ψ(C) Figure 7.9. Degree of recoverable anisotropy in elastic-plastic materials in triaxial strain compression for linearly pressure-dependent elastic moduli for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. With the exception of alumina, the degree of recoverable anisotropy never exceeds 5% even in the finite strain regime. 62 0 0.03 0.06 0.09 0 0.2 0.4 0.6 0.8 1 εA ψ(C) Figure 7.10. Degree of recoverable anisotropy in elastic-plastic materials in triaxial strain compression for elastic moduli given by Mie Gr¨uneisen equation of state for TiB2 ---, limestone −−−−−−, and 2024 aluminum −·−·−·−·−. Interestingly, the degree of recoverable anisotropy in limestone exceeds 8% when the bulk modulus is computed from the Mie Gr¨uneisen equation of state while the degree of recoverable anisotropy does not exceed 6% and 1% for TiB2 and 2024 aluminum, respectively. CHAPTER 8 EFFECTS OF RDIA ON ISOTROPIC MEDIA 8.1 Introduction Thermodynamic admissibility of continuum models in the sense that the stress is derivable from a thermodynamic energy potential has received considerable attention in the modeling of elastic-plastic materials, particularly in the geomechanics com-munity (Collins and Houlsby, 1997; Houlsby and Puzrin, 1999; Puzrin and Houlsby, 2001; Houlsby et al., 2005). In these "hyperplastic" materials, so called because of their similarity with hyperelastic materials, the stress is derivable from a strain energy potential up to the elastic limit and from a combination of the strain energy potential and a dissipation potential beyond. A consequence of the continuum thermomechanics framework is that the elastic stiffness, found by twice differentiating the elastic strain energy potential with respect to the strain tensor, must develop induced anisotropy in response to deformation, even for isotropic materials (Marsden and Hughes, 1994). As pointed out by Houlsby et al. (2005), deformation-induced anisotropy is not inherent to the material, but is a reversible change in the integrity basis of the elastic stiffness tensor in response to deformation. This subtle fact, a consequence of twice differentiating an isotropic function with respect to a tensor, has been documented in the geomechanics literature (Hueckel et al., 1992; Collins and Houlsby, 1997; Houlsby et al., 2005) and texts on applied mechanics (Marsden and Hughes, 1994; Simo and Hughes, 2000), but has seen little consideration in other engineering disciplines. For the secant stiffness of isotropic materials, Hueckel et al. (1992) demonstrated the conditions under which deformation-induced anisotropy is nonexistent: bulk modulus varies at most with pressure and constant shear modulus. These same 64 conditions also hold for the tangent stiffness tensor, as was shown in Chapter 5. Years earlier, Zytynski showed that the adoption of a pressure-dependent shear modulus in an isotropic elastic stiffness formulation leads to energy production in paths closed in elastic stress (Zytynski et al., 1978). In other words, if an isotropic stiffness is used in the modeling of the response of an isotropic material in conjunction with a nonconstant shear modulus, the elastic stiffness is not derivable from an elastic energy potential and the resulting model is thermodynamically inadmissible. The effects of deformation-induced anisotropy are not apparent in all materials. In brittle materials, or materials whose strength is not strongly dependent on pressure, such as most metals, plastic flow and material failure are likely to occur before the magnitude of deformation-induced anisotropy becomes significant. However, in compliant materials or materials whose strength and elastic moduli strongly depend on pressure, the magnitude of deformation-induced anisotropy is capable of attaining significant levels. This can, at least in part, explain why RDIA has seen so little attention outside of the geomechanics community where strongly pressure-dependent materials are characterized. Though the thermodynamic requirement of deformation-induced anisotropy is known, the effects of RDIA on isotropic materials have not been described in the lit-erature. In this chapter, effects and consequences of RDIA on elastic-plastic materials will be described. The effects will be characterized by comparing model predictions for the thermodynamically inadmissible model excluding RDIA and the the admissible model including induced anisotropy. Differences between the two models will be shown by first fitting each model to the same experimental characterization data so that they agree in that loading mode. With the two models aligned in this way, the differences between the two models will be quantified in other loading modes. Throughout this chapter, quantities of the inadmissible model are distinguished by a raised tilde, while those of the admissible model maintain the overbar, heretofore employed. 65 8.2 Elastic Response It is well known that the hydrostatic and deviatoric responses of inherently anisotr-opic materials are coupled. Lesser well known, is the coupling of hydrostatic and deviatoric responses of isotropic materials due to RDIA. In this section, the coupling of hydrostatic and deviatoric behavior in the elastic response of isotropic materials will be shown. Also shown is the effect that the coupling has on the interpretation of elastic characterization data. Up to the elastic limit, the stress response of an elastic material may be computed directly from Eq. (5.1) σ = ∂u ( , η) ∂ η , requiring only that the elastic energy potential u be known. The simplest form of the energy potential capable of satisfying Eq. (5.11), found by direct integration of Eq. (5.11b) with respect to e, is u(J e 1 , J e 2 , η) = ˆu(J e 1 ) + 2μ(J e 1 )J e 2 + ´u(η), (8.1) from which σ = σzEz + σrEσ r (8.2) where σz and σr are given by σz = √3 ∂u ∂J e 1 = √3 ˆu + 2μ J e 2 σr = ∂u ∂γe = 2μ 2J e 2 (8.3) and the unit tensor Eσ r by Eσ r = ∂u/∂γe ∂u/∂γe (8.4) The subscripts "z" and "r" refer to stress measures that are isomorphic to stress space and can be interpreted geometrically as cylindrical coordinates in stress space (cf. Brannon, 2007). 66 For the inadmissible model, σ = ˜σzEz + ˜σrE σ r (8.5) ˜σz = √3˜κJ e 1 ˜σr = 2˜μ 2J 2 (8.6) where E σ r = γe 2J 2 (8.7) 8.2.1 Coupled Hydrostatic-Deviatoric Response Letting ˆσz = √3ˆu and expressing J e 2 in Eq. (8.3) in terms of σr, J 2 = σ2r /8μ2, the coupled thermodynamic pressure is σz = ˆσz + √3μ 4μ2 σ2r (8.8) The presence of both σz, which is a multiple of pressure, and σr, which is the magnitude of the stress deviator, in Eq. (8.8) demonstrates that the hydrostatic and deviatoric responses for an isotropic elastic material are coupled. The coupling arises naturally as a consequence of taking the derivative of the isotropic elastic energy potential. For the case of the inadmissible material, coupling of the hydrostatic and devia-toric responses is not predicted, so that ˜σz = ˆσz, and a relative measure of the change in pressure due to the deviatoric response may be obtained by φz = σz − ˜σz ˜σz ≈ 4√3μ 2˜σ z 0J e 1 +σ˜z 0J e 1 2 J e 2 , (8.9) where ˜σz has been approximated by its second order Maclaurin series expansion about J 1 = 0 (the subscript 0 refers to the argument evaluated at J 1 = 0). In practice, derivatives of elastic moduli are measured with respect to stress and not strain, thus, Eq. (8.9) can be cast in a more useful form through the relation (d˜σz/dJ e 1 )0 = −1/√3κ0 φz = 12(dμ/dˆσz)0 2√3J e 1 − (dκ/d˜σz)0J e 1 2 J 2 (8.10) 67 8.2.2 Elastic Moduli When fitting the thermodynamically admissible and inadmissible models to the same elastic characterization data, the coupling of hydrostatic and deviatoric re-sponses in isotropic materials requires that the apparent pressure sensitivity of elastic moduli be altered so that the two responses are the same in their axial component, the component most often measured. To see this, consider the rate of stress in the admissible and inadmissible models ˙σ = κ˙ vδ + 2μ ˙γ + 2μ ˙ vγe + γ:γ˙ δ (8.11) ˙ σ = ˜κ˙ vδ + 2˜μ˙γ (8.12) Specializing Eq. (8.11) for the special case of triaxial strain deformation, where the strain rate is given by ˙ = (˙ , ˙ /c, ˙ /c, 0, 0, 0)T, the axial and lateral stress rate components are ˙σ A = c + 2 c κ + 4 3 c − 1 c μ + 2c2 − c − 1 c2 μ ˙ σ ˙ A = c + 2 c ˜κ + 4 3 c − 1 c ˜μ ˙ (8.13) ˙σ L = c + 2 c ˜κ − 2 3 c − 1 c ˜μ + c2 − 5c + 4 c2 ˜μ ˙ σ ˙ L = c + 2 c ˜κ − 2 3 c − 1 c ˜μ ˙ (8.14) In practice, constitutive models are calibrated to the measured axial response to deformation. Adopting equivalent bulk modulus functions in each model κ = ˜κ = κ, and assuming that the shear modulus is proportional to volume change (or, by the chain rule, proportional to pressure) μ = μ0 + μ0 J 1 , (8.15) equating axial responses of both models so that each fits the same calibration data requires that μ = 2 + c 3(1 + c) ˜μ . (8.16) 68 Consider the special case of uniaxial strain deformation where c = ∞. Taking the limit of Eq. (8.16) as c→∞, μ = 1 3 ˜μ . (8.17) As seen in Eq. (8.17), neglecting RDIA leads to an apparent first order shear modulus term that is a factor of three times larger in the inadmissible model than in the admissible model when fit to uniaxial strain data. When fit to data for other values of c (other types of triaxial strain deformation), the factor of three is changed, but the following conclusion remains unchanged: the measured pressure-dependence of the shear modulus is larger in the inadmissible model than in the admissible model due to neglecting the implicit coupling of isotropic and deviatoric responses predicted by RDIA. This fact will be shown to have significant consequences when considering plastic deformation. Returning now to the admissible and inadmissible models in uniaxial strain de-formation and making the substitution μ = 1/3˜μ , the axial and lateral stress rate components reduce to ˙σ axial = κ + 4 3 μ0 + ˜μ 0 ˙, σ ˙ axial = κ + 4 3 μ0 + ˜μ 0 ˙ , ˙σ lateral = κ − 2 3 μ0 ˙, σ ˙ lateral = κ − 2 3 μ0 + ˜μ 0 ˙ . (8.18) Even though the axial stress response, by design, is equivalent for each model, the lateral stress responses differ in the first order. Thus, if elastic moduli are fit to the uniaxial strain experimental data (common in large deformation mechanics) such that the two models predict the same axial response (as was done here), they differ in their first order predictions of the lateral response and, in fact, would differ in any other distortional loading mode. When deformed uniaxially, the influence of the deformation-varying shear modulus on the elastic stiffness is only on the axial component. Physically, shear stiffening occurs only in the direction of loading while remaining unchanged in the lateral directions, which is not the case in the isotropic stiffness. The behavior of C is in agreement with the thought experiment of Chapter 1 in which the collection of foam spheres stiffened in its axial response to deformation and not its lateral response. 69 8.2.3 Elastic Limit Of course, the equations of the previous section are valid only up to the elastic limit which, for uniaxial strain deformation, is often referred to as the Hugoniot Elastic Limit (HEL). In this case, J 1 = HEL, J 2 = 1/3 2HEL, and the predicted difference in pressure at the elastic limit can be found by taking the first order Maclaurin series of Eq. (8.10) and approximating HEL by σHEL/H0, where σHEL is the axial stress at the Hugoniot elastic limit and H0 = κ0 + 4/3μ0 is the initial constrained modulus φHEL z ≈ 2 √3 dμ d˜σz σHEL H0 . (8.19) As expected, if μ = 0, φHEL z = 0, and the admissible model reduces to the inadmissible model. The corresponding difference in σr at the Hugoniot elastic limit σHEL r is φHEL r = σHEL r − ˜σHEL r ˜σHEL r ≈ − A4κ0σHEL A1H0 φHEL z . (8.20) In the admissible model, the increase in pressure due to the hydrostatic/deviatoric coupling, causes an apparent increase in the shear strength for materials whose yield strength varies with A4. 8.3 Inelastic Response Beyond the elastic limit, path dependence of the inelastic response requires that incremental relations relating stress and strain be employed. In the following section, incremental equations relating stress and elastic strain for the admissible and inadmis-sible model will be shown. The formulation is based on the additive decomposition of the rate of strain into elastic and plastic parts and satisfying the Kuhn-Tucker conditions Lubliner (2008). Operating on the rate of elastic strain by the reduced elastic stiffness C in Eq. (5.10), the rate of stress is ˙σ = C:˙ e = ˙σ zEz + ˙σ rE∗r , (8.21) where ˙σ z, ˙σ r, E∗r are ˙σ z = √3 κJ ˙ e 1 + 2μ 0J˙ e 2 , (8.22) 70 ˙σ r = 2μ 2J ˙ e 2 + 2 μ 0 μ ˙ e v J˙ e 2 + 2 μ 0 μ ˙ e v 2 J e 2 , (8.23) E∗r = 2μ ˙γ e + 2μ 0˙ e vγe 2μ 2J ˙ e 2 + 2 μ 0 μ ˙ e v J˙ e 2 + 2 μ 0 μ ˙ e v 2 J e 2 . (8.24) For inelastic loading, the stress rate is found by first computing the trial elastic stress, using Eq. (8.22) - Eq. (8.24), and then projecting it back to the yield surface, as shown in the following sections. 8.3.1 Trial Stress Increment In classical plasticity, explicitly determining the elastic strain rates required for the evaluation of the stress rate in Eq. (8.21) is not necessary. Instead, the stress rate is found by projecting the trial stress onto the zero isosurface of the yield function such that the yield criterion is satisfied. Assuming that the strain rate over the entire step is elastic, the trial stress rate is given by ˙σ trial = C:˙ = ˙σ trial z Ez + ˙σ trial r E∗r trial, (8.25) with ˙σ trial z , ˙σ trial r , and E∗r trial given by Eq. (8.22) - Eq. (8.24), ˙σ trial z = √3 κJ ˙ 1 + 2μ 0γe:˙γ , (8.26) ˙σ trial r = 2μ 2J ˙ 2 + 2 μ 0 μ ˙ v ˙γ :γe + 2 μ 0 μ ˙ e v 2 J e 2 , (8.27) E∗r trial = 2μ ˙γ + 2μ 0˙ vγe 2μ 2J ˙ 2 + 2 μ 0 μ ˙ v ˙γ :γe + 2 μ 0 μ ˙ e v 2 J e 2 . (8.28) If instead, the elastic stiffness were presumed isotropic, σ ˙ trial z , σ ˙ trial r , and E ∗r trial would be given by σ ˙ trial z = √3κJ ˙ 1 , (8.29) σ ˙ trial r = 2˜μ 2J ˙ 2 , (8.30) E ∗r trial = 2˜μ˙γ 2˜μ 2J ˙ 2 . (8.31) 71 The relative differences between the r and z components of ˙σ trial and ˙ σ trial are φtrial ˙z = ˙σ trial z − σ ˙ trial z σ ˙ trial z , (8.32) φtrial ˙r = ˙σ trial r − σ ˙ trial r σ ˙ trial r , (8.33) which, for the case of uniaxial strain deformation, reduces to φtrial ˙z = 4 3 μ 0 κ ( e axial − e lateral) , (8.34) φtrial ˙r = (μ0 + μ 0 (2 e axial − e lateral))2 (μ0 + ˜μ 0 ( e axial + 2 e lateral))2 − 1. (8.35) For compressive loading, it is easily shown that φ˙ztrial > 0 and φ˙rtrial < 0 meaning that the admissible model predicts greater pressure and less shear stress at a given strain than the inadmissible model. This fact will be highlighted in Section 8.4 where specific numerical examples will be considered. 8.3.2 Return Direction Neglecting elastic-plastic coupling, the direction of return from the trial stress increment to the actual stress increment on the yield surface is given by p = C:m, (8.36) which, for the admissible model, gives p = pzEz + prEp r , (8.37) where pz = √3 κ trm + 2μ 0γe:m , (8.38) pr = 2μ 2Jm 2 + 2 μ 0 μ trm m :γe + 2 μ 0 μ trm 2 J e 2 , (8.39) Ep r = 2μm + 2μ 0 trmγe 2μ 2Jm 2 + 2 μ 0 μ trm m :γe + 2 μ 0 μ trm 2 J e 2 . (8.40) 72 For the inadmissible model, the return direction is given by ˜p = ˜pzEz + ˜prE˜p r , (8.41) where ˜pz = √3κ trm, (8.42) ˜pr = 2˜μ 2Jm 2 , (8.43) E˜p r = 2˜μm 2˜μ 2Jm 2 . (8.44) The angle between the admissible and inadmissible return directions can be approx-imated by θ ≈ cos−1 pzp˜z + prp˜r p2z + p2r ˜p2z + ˜p2r . (8.45) A convenient reduced form of Eq. (8.45) does not exist, even for the simple case of uniaxial strain deformation. In the next section the angle between return directions will be investigated further. 8.4 Numerical Examples The effects of deformation-induced anisotropy will now be shown for the case of uniaxial strain and triaxial strain loading to pressures up to 250 MPa. In the following examples, the bulk modulus for the admissible and inadmissible models will each be computed by using the Mie Gr¨uneisen equation of state κ = 1 ρ 1 − γ (ρ0 − ρ) 2ρ0 − dp d(1/ρ) H + γ 2 pH, (8.46) where pH is the pressure on the Hugoniot. The shear modulus is given by μ = μ0 + μ 0J 1 , (8.47) and was calibrated as previously outlined for each model. Yield behavior is by linear Drucker-Prager plasticity with the following yield function f (σ) = Jσ 2 − a1 + a4Jσ 1 . (8.48) Material properties used in each simulation are given in Table 8.1 and are similar to those used by Salari et al. (2004) except that nonlinear elasticity is permitted by allowing μ to vary with deformation. 73 8.4.1 Uniaxial Strain Compression 8.4.1.1 Elastic Limit The difference in the elastic limit predicted by Eq. (8.19) and Eq. (8.20) is negligible for this material. For titanium diboride, the difference in the z and r components at the elastic limit are 1% and .6%, respectively, as shown in Figure 8.1. 8.4.1.2 Trial Stress Increments Inelastic behavior is shown in Figure 8.2. In the figure, the path through stress space is shown with the trial stress increment for the admissible and inadmissible models at each iteration. A detailed step is shown in Figure 8.3. As predicted by Eq. (8.34) and Eq. (8.35), the trial stress increment is greater in the z component and smaller in the r component in the admissible model than in the inadmissible model. The result is that the trial stress increment is more tangential to the yield surface in the admissible model. This behavior is readily explained in the context of isotropic/deviatoric coupling in isotropic materials required by thermodynamics: the calibration procedure of the two models results in a decrease of the first order shear modulus term for the admissible model which in turn leads to a decrease in the increment in the trial stress deviator, while at the same time, shear induced bulking causes an increase in the isotropic part of the trial stress. As shown in Figure 8.3, the final resultant stress state is nearly the same in each model, though achieved through different paths. Thus, differences in predicted stress values are likely to be similar in both models. However, since the plastic strain increment is proportional to distance of projection, which is clearly different, the two models have significant differences in their parti-tioning of strain into elastic and plastic parts. The difference in plastic strain is then expected to lead to differences in predicted plastic work and hardening evolution behavior. 8.4.1.3 Permanent Volume Change Turning attention to Figure 8.3, the direction of return of the admissible trial stress increment to the yield surface is more hydrostatic than in the inadmissible model. 74 Despite this fact, the prediction of permanent volume change in the admissible model is less than in the inadmissible model, as shown in Figure 8.4, where permanent volume change is shown plotted versus σz. Referring to Figure 8.3, the trial stress increment in the admissible model is more tangential to the yield surface meaning that a larger portion of the strain increment is elastic, thus, the predicted amount of permanent volume change is correspondingly less than in the inadmissible model. When compared with experimental data, model predictions assuming a classical elastic-plastic framework and associative plastic flow often over-predict the magnitude of plastic work and permanent volume change realized by the material under consid-eration (Seseny et al., 1983). This observation has led to the adoption of nonclassical features in constitutive models to bring congruence between numerical prediction and experimental observation. To account for the over-prediction of permanent volume strain in pressure-dependent plasticity theories, nonassociative plastic flow rules are often adopted. In nonassociative plasticity, instead of taking the direction of plastic flow to be normal to yield function zero isosurface, it is taken as normal to another surface, known as the flow surface, which is taken to be a level set of a flow function ϕ. In nonassociative plasticity, the flow potential typically has the same functional form as the yield function. For example, a flow potential often used with the Drucker-Prager yield function in Eq. (8.48) is ϕ (σ) = Jσ 2 − a1 + βJσ 1 , (8.49) where the parameter β is chosen such that model predictions of the permanent volume change match what is experimentally observed. In metals plasticity, where permanent volume strain is typically smaller than experimental uncertainty (Richmond and Spitzig, 1980), β = 0 and the direction of plastic flow is purely deviatoric. For geomaterials, where permanent volume change is observed, 0 ≤ β ≤ a4. Salari et al. (2004), using similar material properties as in Table 8.1, chose β = a4/2 = 0.115 for best model performance. Using this value of β in the inadmissible model, the difference between permanent volume strain for associative and nonassociative flow is shown in Figure 8.5. Also shown is nonassociative flow in the admissible model except β = 0.14 so that permanent volume strain predictions from each model are in 75 agreement. As seen in Figure 8.4, and due to the decrease in projection distance from the trial to actual stress rate in the admissible model, the admissible model appears to be nonassociative when compared to the inadmissible model due to the decrease in predicted permanent volume strain, even though each model was associative and fit to the same elastic characterization data. When considering nonassociative flow, achieving comparable degrees of permanent volume strain in each model required that β be 22% larger in the admissible model than in the inadmissible model. Thus, even for nonassociative flow, the admissible model requires less nonassociativity than the inadmissible model to obtain similar permanent volume strain prediction. As a result of the reduced permanent volume strain predicted by the admissible model, the predicted plastic work Πp will also be reduced. Shown in Figure 8.6 is the ratio of the plastic work of the admissible model and inadmissible model, showing the reduction in predicted plastic work in the admissible model. The partitioning of plastic work into elastic and plastic parts is an active research area with impor-tant ramifications in the prediction of temperature response of plastically deformed materials (Kim et al., 1990; Hodowany et al., 1999; Rosakis et al., 2000). Often, predictions based on an isotropic stiffness over-predict the amount of plastic work performed. Con