Prepared by

Oak Ridge Institute for Science and Education

ORAU Maryland

4692 Millennium Drive, Suite 101

Belcamp, MD 21017

Under contract: 1120-1120-99

Army Research Laboratory

Implementation of a Mechanochemical Model for Dynamic

Brittle Fracture in SIERRA

by Adam Sokolow

ARL-CR-0742 August 2014

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the Army position unless so

designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official endorsement or approval of the use

thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5066

ARL-CR-0742 August 2014

Implementation of a Mechanochemical Model for Dynamic

Brittle Fracture in SIERRA

Adam Sokolow

Oak Ridge Institute for Science and Education

Approved for public release; distribution is unlimited.

Prepared by

Oak Ridge Institute for Science and Education

ORAU Maryland

4692 Millennium Drive, Suite 101

Belcamp, MD 21017

Under contract: 1120-1120-99
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering 

and maintaining the data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, 

including suggestions for reducing the burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson 

Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to 

comply with a collection of information if it does not display a currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To) 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

4. TITLE AND SUBTITLE 

5c. PROGRAM ELEMENT NUMBER 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

6. AUTHOR(S) 

5f. WORK UNIT NUMBER 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 

 REPORT NUMBER 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S) 

11. SPONSOR/MONITOR'S REPORT 

 NUMBER(S) 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF: 

19a. NAME OF RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 

17. LIMITATION 

OF ABSTRACT 

18. NUMBER 

OF PAGES 

19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18

August 2014 Final

Implementation of a Mechanochemical Model for Dynamic Brittle Fracture in

SIERRA

ARL-CR-0742

Approved for public release; distribution is unlimited.

September 2013 – April 2014

1120-1120-99

Adam Sokolow

U.S. Army Research Laboratory

ATTN: RDRL-WMP-B

Aberdeen Proving Ground, MD 21005-5066

This report discusses the implementation of a mechanochemical model for brittle fracture in the government code SIERRA. It is

largely based on damage models for which damage is an internal state variable. The model is introduced first in a linear elastic

context and then a similar model is developed using finite deformation theory. The model is then verified and some example

applications discussed.

brittle fracture, mechanochemical, dynamic failure, damage evolution, SIERRA

52

Adam Sokolow

Unclassified Unclassified Unclassified 410-306-2985

UU

ii
Contents

List of Figures v

List of Tables vi

Acknowledgments vii

1. Introduction 1

2. Mechanochemical Model: Linear Elastic Framework 3

3. Mechanochemical Model: Finite Deformation Formulation 5

3.1 Finite Deformations and the Log-Strain Tensor .......................................... 5

3.2 Mechanochemical Constitutive Model.................................................... 8

4. Implementing the Damage Evolution Equations 10

4.1 Dynamic Evolution of Damage ........................................................... 10

4.2 Instantaneous Damage ..................................................................... 11

4.3 Additional Details of the Implementation ................................................ 11

5. Model Verification 12

5.1 Instantaneous Damage Evolution, Confined-Compressive Deformation................ 12

5.2 Instantaneous Damage Evolution, Isochoric Shear Deformation ........................ 15

5.3 Dynamic Damage Evolution, Confined-Compressive Deformation..................... 17

iii
6. Example Problems 19

6.1 One-Dimensional Dynamic Confined-Compression ..................................... 19

6.2 Two-Dimensional Radial Damage Bands................................................. 20

7. Future Considerations 23

7.1 Removing Healing From the Mechanochemical Model ................................. 23

7.2 Decoupling Bulk Damage and Deviatoric Damage Functions........................... 24

7.3 Damage Growth Asymmetry in Compression and Tension.............................. 25

8. Conclusions 26

9. References 27

Appendix A. Analytical Solution to Kinetic Equation 29

Appendix B. SIERRA Implementation 33

Distribution List 41

iv
List of Figures

Fig. 1 Single element compression test. Comparison of the theoretical prediction (solid

lines) against the simulation results (symbols). Panels a and b hold ξ/μ = 0.1 and vary

cmin from 0.2 to 1. Panels c and d hold cmin = 0.3 and vary ξ/μ from 0.2 to 1. Stress

made positive for viewing purposes although its value is negative in compression.. ........ 14

Fig. 2 Single element shear test. Comparison of the theoretical prediction (solid lines)

against the simulation results (symbols). Panels a and b hold ξ/μ = 0.1 and vary cmin

from 0.2 to 1. Panels c and d hold cmin = 0.3 and vary ξ/μ from 0.2 to 1.................... 16

Fig. 3 Single element dynamic compression test. Comparison of the theoretical prediction

(solid lines) against the simulation results (symbols). Panels a and b hold ξ/μ = 0.1

and vary cmin from 0.2 to 1. Panels c and d hold cmin = 0.3 and vary ξ/μ from 0.2 to 1.. .. 18

Fig. 4 A One-Dimensional alignment is slowly compressed from its Z = 0 end. Each

panel plots the damage vs. position within the alignment, where the plots are on the

initial configuration. The time for each panel is shown in the upper right corner and the

time difference between each panel is not constant.. .......................................... 20

Fig. 5 Two-Dimensional instability and the development of radial damage bands. A 2-D

annular disc is slowly compressed from its outside edge r = Ro. Each panel plots the

damage in false-color as it varies with position within the disc. Plots are on the initial

configuration. The elapsed time between each panel is 175 s. At first the damage

accumulates symmetrically before a material instability produces the localization of

damage along radial lines.. ...................................................................... 22

v
List of Tables

Table 1 Summary of material parameters used in the 1-D stability example.. ................... 19

Table 2 Summary of material parameters used in the 2-D radial damage band instability

example.. ......................................................................................... 21

vi
Acknowledgments

We would like to thank Mike Scheidler and Rich Becker for offering their expertise during

various aspects of this project and of course Michael Grinfeld for introducing me to this problem.

This research was supported in part by an appointment to the Postgraduate Research Participation

Program at the U.S. Army Research Laboratory (ARL) administered by the Oak Ridge Institute

for Science and Education through an interagency agreement between the U.S. Department of

Energy and ARL.

vii
INTENTIONALLY LEFT BLANK.

viii
1. Introduction

The fracture and failure of brittle materials has remained a core focus in Army related sciences.

Of particular interest is the destabilization of an axisymmetric solution in problems that are

completely axisymmetric. This destabilization problem largely motivated Grinfeld and Wright1to

study a model for dynamic failure that depends on the interplay between two physical effects:

internal elastic energy, and the energy associated with breaking chemical bonds. Because of these

two considerations, the model is referred to as a mechanochemical model for fracture. The

motivation of the mechanochemical model is based on Stress Driven Rearrangements Instabilities

of phase interfaces Grinfeld2and is summarized in Grinfeld3and Kassner et al..4 Motivated by

the Stress Driven Rearrangements Instabilities, the model was then used in a different context to

form the mechanochemical model, Grinfeld and Wright,1 which was numerically implemented in

MATLAB to handle quasi-static loading cases. Grinfeld et al.5 used this implementation in

MATLAB to established the appearance of radially damaged zones and the destabilization of the

axisymmetric solution.

A clear and condensed presentation of damage theory is given by Kachanov,6as well by

Chaboche,7, 8 and a short introduction to damage as an internal state variable can also be found in

the book by Holzapfel.9In the types of damage models considered by Kachanov, the energy

density equation is modified by a reduction factor or a damage function that depends on an

internal state variable, the damage. As one might expect, the damage reduces the internal energy

thereby reducing the stress response of the damaged material. In the mechanochemical model

discussed here, the damage also contributes to the internal energy. This addition to the energy

potential creates an interplay between mechanical and chemical constituents and produces the

source of the instability.

The original intent was to implement the mechanochemical model discussed in Grinfeld and

Wright1into the government code SIERRA. The model discussed in Grinfeld and Wright1is

based on a modified linear elasticity. In this report, however, we focus on a finite deformation,

nonlinear formulation of the problem. This change was required by constraints introduced within

the finite element solver SIERRA. In SIERRA, constitutive model designs are split between rate

integrated and hyperelastic formulations. However, the features of the mechanochemical model

make the decision between the two nontrivial. The mechanochemical model largely depends on

the internal elastic energy, which makes it a good candidate for the hyperelastic approach. The

desired linear elastic response on the other hand, would be simplest to implement in the rate

1
integrated formulation. At first, a rate integrated version of the mechanochemical model was

implemented in SIERRA. This initial implementation determined the linear elastic energy

through numerical integration but was prone to error. As a compromise, an alternate formulation

of the model was developed that mimicked the linear elastic response using the log-strain tensor.

This method benefited from an accurate calculation of elastic energy and has proved to be more

stable. The model and its implementation in SIERRA are discussed in this report. Unlike the

previous implementation in MATLAB,5 which could only handle quasi-static loading, this

implementation allows for the dynamic problem of failure to be studied on a massively parallel

computational architecture.

The report is organized as follows. We discuss the general physical features and introduce the

system of equations that describe the damageable material in section 2. This is done in the

framework of the linear elastic theory and largely recapitulates the work of Grinfeld and Wright.10

In section 3 we introduce the log-strain tensor and derive a new formulation of the

mechanochemical model for finite deformations. Section 4 presents the numerical details of how

the damage evolution is handled and discusses some of its subtleties. In section 5 we verify the

implementation of our model in SIERRA using some single element tests. A few example

problems that illustrate applications of the mechanochemical model are discussed in section 6,

and possible future alterations to the model are discussed in section 7. We make concluding

remarks in section 8 and derive the equations of damage evolution in appendix A. Appendix B

includes the code used to implement the model in SIERRA.

2
2. Mechanochemical Model: Linear Elastic Framework

In this section we largely recapitulate the work by Grinfeld and Wright10 but with additional

comments to motivate the following section where we introduce a finite deformation

mechanochemical model. We begin by introducing the terms considered in our energy density

potential from which we derive the elastic stress and the equations that govern the chemical

kinetics.

We take the total energy density to be a function of the mechanical deformation (or strain) and the

damage. Specifically, we assume the energy can be decomposed additively into a mechanical part

and a chemical part, i.e.,

e = emechanical + echemical . (1)

Holding the damage fixed, the derivative of the energy with respect to the strain gives the stress,

∂e

∂ strain









damage constant

= stress. (2)

Similarly, holding the strain constant and taking the derivative of energy with respect to the

damage gives the driving force for damage evolution,

∂e

∂κ









strain constant

∝

∂

∂tdamage . (3)

Now we define the mechanical and chemical terms for a linear elastic material. We assume the

mechanical energy is reduced by a damage function φ, which depends on the damage κ so that the

mechanical energy is given by

emechanical = φ(κ) eelastic . (4)

Here emechanical is the mechanical energy of a damageable material, and eelastic is the elastic energy

associated with a hypothetical undamageable material undergoing the same deformation. We also

assume the damage function φ(κ) reduces the elastic energy with increasing damage. Specifically,

we take our damage function φ(κ) to be

φ(κ; κ

∗

, cmin) = (

1 − (1 − cmin)

κ

κ∗ : κ ≤ κ

∗

cmin : κ > κ∗

. (5)

The functional form of equation 5 is a simple linear ramp between undamaged and a maximally

3
damaged material. The parameters κ

∗

is the maximal damage value. If while updating the value of

κ the damage exceeds κ

∗

, it is simply set to κ

∗

(see section 4). The parameter cmin is the lower

limit to the reduction factor of a fully damaged material.

Using a linear elastic theory for the elastic energy term, we obtain the following energy density

per unit reference volume:

emechanical = μφ(κ)



ν

1 − 2ν

u

i

.|iu

j

.|j + u(i|j)u

ij

.|.



. (6)

Here uij is the linear elastic strain tensor,

uij ≡

1

2



∂ui

∂xj

+

∂uj

∂xi



, (7)

where uiis the displacement vector of the deformation.

The second term in the energy density is the chemical energy. This energy is assumed to be

quadratic and of the form

echemical =

ξ

2

(κ − κ

o

)

2

. (8)

Here ξ is a chemical constant with dimensions of stress and κ

o

the damage associated with

minimum energy.

Combining these two terms into an energy density, we have the so-called Kachanov-Lifshitz

function.

e(ui|j, κ) = μφ(κ)



ν

1 − 2ν

u

i

.|iu

j

.|j + u(i|j)u

ij

.|.



+

ξ

2

(κ − κ

o

)

2

. (9)

Thus, by considering changes in the energy density, and either holding the strain or the damage

constant, we obtain:

∂e

∂uik









κ constant

= 2μφ(κ)



ν

1 − 2ν

ullδik + uik≡ σik (10)

and

∂e

∂κ









uik constant

=

∂φ

∂κeelastic + ξ (κ − κ

o

) (11)

=

∂φ

∂κμ



ν

1 − 2ν

u

m

. |mu

n

.|n + u(m|n)u

mn

.|.



+ ξ (κ − κ

o

) . (12)

The stress tensor retains its usual form, i.e., it is still defined in terms of the derivative of the

4
energy with respect to the strain tensor; however, it is modified by the damage function φ(κ).

Similarly, a generalized-stress is derived from the chemical energy. It has units of stress and

drives changes in the damage. This generalized-stress is referred to in the literature as a chemical

potential. The following two equations result from combining the stress equation 10, with

momentum balance, and the chemical potential 11 with a kinetic reaction law:

ρ

∂

2ui

∂t2

=

∂σik

∂xk

=

∂

∂xk



2φ(κ)μ



ν

1 − 2ν

ullδik + uik , and (13a)

∂κ

∂t = −K



∂φ

∂κμ



ν

1 − 2ν

u

m

.|mu

n

.|n + u(m|n)u

mn

.|.



+ ξ (κ − κ

o

)



. (13b)

The dimensions of K are:

[K] = 1

[Time][Stress]

. (14)

The specific choice of equations 13 contain two physically questionable features, which we

discuss in later sections. Namely, a degrading bulk modulus and reversible damage. Reversible

damage used in this report is not consistent with the work by Grinfeld and Wright10 where the

derivative of the damage is forced to remain positive. See section 7.1 for details regarding the

minor alteration in the constitutive model that disables healing.

3. Mechanochemical Model: Finite Deformation Formulation

As discussed in the introduciton, the original goal of the research was to implement the

mechanochemical model for brittle fracture described in the previous section by equations 5

and 13 as a material model in SIERRA. However, the choice to implement the constitutive model

in SIERRA carried with it a number of constraints. In this section we derive a constitutive relation

for an isotropic hyperelastic material (hyperelastic for fixed damage) that captures the same

physical features described in the previous section. This does not represent the most general

nonlinear model, and in fact is a special case that was chosen to mimic the linear elastic

mechanochemical model when used in the small strain limit.

3.1 Finite Deformations and the Log-Strain Tensor

The following briefly covers some modern continuum mechanics and introduces the log-strain

tensor. Two good references on modern continuum mechanics are the text books by Holzapfel,9

and Truesdell and Noll.11 The Cauchy stress tensor T is a linear transformation of a direction

vector to a traction force in the current configuration. It can be written in terms of the spherical

5
part, i.e., the pressure p, and its deviatoric part T

∗

:

T = −pI + T

∗

, p ≡ −

1

3

tr T . (15)

We use the convention that normal stress components are positive in tension. The deviatoric part

of a tensor, denoted by a superscript ∗, is defined as follows:

A∗ ≡ A −

1

3

tr [A] I , (16)

where I is the identity tensor. The deformation gradient F is a linear transformation from the

reference configuration to the current configuration, which can be expressed in terms of a

properly orthogonal rotation matrix R and positive definite symmetric matrices U or V , called

the right- and left-stretch tensors, respectively, using the polar decomposition

F = RU = V R . (17)

Both U and V share the same eigenvalues λi called the principal stretches. The Jacobian J is the

determinant of the deformation gradient,

J ≡ detF = detV , (18)

and is equal to the ratio of the current specific volume to the reference specific volume. The right

and left Cauchy-Green deformation tensors, respectively, are defined as

C = F

TF = U2

and (19)

B = F FT = V

2

. (20)

Since C and B are symmetric and positive definite, the eigenvalues of C and B are λ

2

i

. The

principal invariants of a second-order tensor A are given by

IA = tr(A), IIA = tr(A2), and IIIA = tr(A3). (21)

If two tensors share the same eigenvalues one can equate the principal invariants of one to the

other. Thus, for isotropic elastic materials it can be shown that the energy density can only depend

on deformation gradient through the three invariants of B (or equivalently C). In this manner the

energy density of an elastic material can be expressed as some function,

e = eB(IB, IIB, IIIB). (22)

6
There are many choices of finite deformation strain tensors that are consistent with the linear elastic

strain tensor in the infinitesimal limit. These are typically expressed in terms of the stretch tensors,

U or V , or the Cauchy-Green deformation tensors, C or B. We choose to use the log-strain tensor,

also called Hencky strain.12 This particular choice of strain tensor has some desirable features when

compared with real materials at moderate strains13, 14 and is often considered to behave similarly

at moderate strains to an infinitesimal strain measure. The log-strain tensor L is defined as,

L ≡ lnV =

1

2

ln B . (23)

Scheidler15 considered in depth the log-strain tensor and some of its properties. Notably that

tr L = ln J and hence is a measure of the volumetric strain. Furthermore, the deviatoric part, L

∗

,

of L is independent of J and is a measure of shear strain. Since IL = trL = ln J and L

∗

is

traceless, the energy density is some function e

0

L

of another set of invariants:

e = eL(IL, IIL, IIIL) = e

0

L

(ln J, IIL∗ , IIIL∗ ). (24)

Scheidler15 using the log-strain tensor, also showed that the following general relationships for an

isotropic hyperelastic material can be derived for the pressure and the deviatoric stresses:

p = −

∂e(J, L

∗

)

∂J









L∗

, and (25a)

J T

∗ =

∂e(J, L

∗

)

∂L

∗









J

. (25b)

The log-strain tensor is work conjugate to the Kirchoff stress, JT . Work conjugacy follows from

a general kinematic relationship, from which it is shown the diagonal components of the time

derivative of the log-strain tensor are equal to the diagonal components of the rate of deformation

tensor in the principal basis for the left-stretch tensor16 and the assumption of an isotropic elastic

material. An isotropic elastic implies that the Cauchy stress is also diagonal in the principal basis

for the left-stretch tensor, ensuring that the inner product of the Cauchy stress and the rate of

deformation tensor is equal to the inner product of the Cauchy stress and the time derivative of the

log-strain tensor.

Only the second and third principal invariants of L

∗

are nonzero, so that the deviatoric stress can

7
be written:

T

∗ =

2

J



∂e

∂IIL∗



J,IIIL∗

L

∗ +

3

J



∂e

∂IIIL∗



J,IIL∗

[(L

∗

)

2

]

∗

. (26)

The deviatoric stresses depend linearly on the deviatoric strains through the second invariant, but

their dependence on the third invariant is quadratic in the deviatoric strains.15 Thus, in the small

shear-strain limit and when the shear modulus is independent of density, one can conclude:

T

∗ ≈

2μ

J

L

∗

, and (27a)

p ≈ p(J). (27b)

Assuming that the pressure depends only on the Jacobian J, Scheidler15 showed that the strain

energy can be decoupled additively into a function of J only and a function of shear only,

e = evol(J) + eiso(L

∗

). (28)

Taking both of these assumptions—specifically, that the pressure depends only on the Jacobian,

and the deviatoric stress depends linearly on the deviatoric strains—is an enabling step in the

formulation of the finite deformation analog to the linear elastic mechanochemical model. This

formulation is similar to the linear elastic theory since we have effectively eliminated the

pressure-shear coupling that would otherwise exist.

3.2 Mechanochemical Constitutive Model

We now define the decoupled energy density for the mechanochemical model. As before, we take

the energy density per reference volume and additively decouple it into a mechanical and

chemical part:

e = emechanical + echemical , (29)

and express the mechanical energy in terms of a damage function φ and the elastic energy:

emechanical = φ(κ) eelastic . (30)

The elastic energy is now given by equation 28 and the total energy is

e = φ(κ) [evol(J) + eiso(L

∗

)] + echem(κ). (31)

8
We assume the equation of state evol(J) is given by:

evol(J) = K (J ln J − J + 1) , (32)

where K is the bulk modulus; other equations of state could be used in the future.† The energy

associated with the deviatoric deformation is taken to be

eiso(L

∗

) = μ tr (L

∗

)

2



(33)

so that the total energy density function is given by

e = φ(κ)



K (J ln J − J + 1) + μ tr (L

∗

)

2

	 +

ξ

2

(κ − κ

o

)

2

. (34)

Equation 34 is the analog of the Kachanov-Lifshitz function from equation 9 consistent with finite

deformations. Using equation 25 and assuming that the shear strains are small gives

p = −φ(κ)K ln J , and (35)

T

∗ = φ(κ)

2μ

J

L

∗

. (36)

Thus, the total Cauchy stress tensor is

T = φ(κ)



K ln JI +

2μ

J

L

∗



. (37)

The dynamical system analogous to equations 13 for the finite deformation model is

ρ

∂

2u

∂t2

= ∇ · T , and (38a)

−

1

K

∂κ

∂t =

∂φ

∂κ (eelastic(J, L

∗

)) + ξ (κ − κ

o

) . (38b)

Note that eelastic again refers to the elastic energy associated with a hypothetical undamageable

material, recall equation 30. Again there has been no restriction on the sign of the derivative of κ ̇

so that the damage is reversible.

†Typically the (J ln J − J + 1) term in the expression for the energy that relates to the volumetric contribution

does not include the “+1”. In fact, in a hyperelastic model, the functional form of the energy is what is most important

within an arbitrary constant since the stresses depend on derivatives of the energy. However, since we now would like

to relate this energy to a damage growth mechanism, we must ensure that there is no damage growth in the undeformed

state from a non-zero energy (i.e., setting J = 1 must give a zero energy).

9
4. Implementing the Damage Evolution Equations

In this section we discuss how the coupled dynamic system given by equation 38 is handled in the

code and how it is different from the case where the left hand sides of equation 38 are identically

zero, where the damage develops instantaneously. Implementation of the model is in SIERRA

SolidMechanics (Presto/Adagio 4.28; Sandia National Laboratory) a Lagrangian finite element

solver. SIERRA contains both an explicit solver Presto and an implicit, nonlinear preconditioned

conjugant gradient solver, Adagio. The implicit code, Adagio, can run into some difficulties since

the solution can become nonunique as the damage increases. In other words, there can be multiple

deformations that correspond to the same stress, and the solver will no longer converge. Thus,

while the code works with Adagio, its use is not recommended. Since SIERRA is a solid

mechanics solver, the damage model must be introduced without disrupting the numerical

stability. Therefore, there are some additional assumptions discussed here that are needed to solve

the damage evolution alongside the mechanics.

4.1 Dynamic Evolution of Damage

As mentioned in section 2, there are dimensions of time in the kinetic reaction law parameter K.

Thus, there are two time scales that need to be considered when solving the equations of motion.

The first time scale is related to wave propagation speeds and the characteristic dimensions of the

body and the second is how quickly damage accumulates. The stability of the explicit code

depends critically on the elastic wave speeds and the minimum element size, which dictates the

time step used during the calculations. This places some additional requirements on the

calculation of the kinetic reaction law since it too needs to be updated at each time step. If the

kinetic reaction law is also implemented in an explicit manner, there would be a second critical

time step. Since the solutions to the kinetic reaction law are exponentials (see appendix A), it

could become an incredibly stiff system to solve, i.e., sensitive to error. We therefore solve the

kinetic reaction law for an incremental change in time so that we can directly solve for the

updated damage at each time step.

We find that during a small time interval t → t + ∆t, the updated damage κ(t + ∆t) can be

related to a convolution of the undamaged elastic energy with a derivative of the damage function

(see appendix A):

κ(t + ∆t) = κ(t)e

−Kξ∆t − κo



e

−Kξ∆t − 1



− K

Z t+∆t

t

∂φ

∂κ(τ ) eelastic(τ )e

−Kξ(t+∆t−τ)

dτ . (39)

10
We assume that ∆t is chosen by SIERRA such that the elastic energy is varying slowly. This

would appear to be a safe assumption since the time step is chosen for stability. Thus, we pull the

elastic energy out of the integral. If the damage function is given by equation 5, we can also pull

the derivative of the damage function out of the integral and evaluate the result to arrive at the

following:

κ(t + ∆t) ≈ κ(t)e

−Kξ∆t − κo



e

−Kξ∆t − 1



+

1

ξ



1 − e

−Kξ∆t

 1 − cmin

κ

∗

eelastic . (40)

The updated damage κ

n+1 is then:

κ

n+1 = κn

e

−Kξ∆t−κo



e

−Kξ∆t − 1



+

1

ξ



1 − e

−Kξ∆t

 1 − cmin

κ

∗



K (J ln J − J + 1) + μ tr (L

∗

)

2

n+1

,

(41)

where superscripts denote the time step. This updated damage can then be fed into the stress

calculation:

T = φ(κ)



K ln JI +

2μ

J

L

∗



. (42)

The minor alteration to the code to handle the case where the damage is irreversible is discussed

in section 7.1.

4.2 Instantaneous Damage

A special case of damage evolution is when the left hand side of equation 38b is zero. This can be

achieved by taking the limit as K goes to infinity in equation 40, or by directly solving

equation 38b,

κ = κ

o +

1

ξ

1 − cmin

κ

∗



K (J ln J − J + 1) + μ tr (L

∗

)

2

 . (43)

4.3 Additional Details of the Implementation

SIERRA provides the left-stretch tensor V and rotation tensor R from the polar decomposition of

the deformation gradient tensor F at the so-called next time step n + 1. From this updated

configuration, one must determine the Cauchy stress in the unrotated (reference configuration).

As discussed in section 3, the model depends on the log-strain, which requires taking the natural

logarithm of a tensor. We use built-in functions in SIERRA to perform the spectral decomposition

of V , from which the log can be taken. The overall algorithm is outlined as follows:

1. Calculate J and ln J from V

n+1

.

2. Determine eigenvalues/vectors for V

n+1

.

11
3. Calculate the log of the eigenvalues and rotate log V

n+1 from its eigenbasis back to current

basis.

4. Calculate (L

∗

)

n+1 from lnV

n+1

.

5. Update the current undamaged elastic energy, e

n+1

elastic.

6. Solve for the new damage value κ

n+1, restrict its range to 0 < κ < κ∗ using equation 41

or 43 as selected by the user.

7. Determine the value of the damage function φ(κ

n+1).

8. Determine T

n+1 from equation 42.

9. “Unrotate” T

n+1 to reference configuration by calculating RTT R.

10. Repeat steps 1–9 for all elements in the material.

5. Model Verification

This section compares theoretical predictions of the model against simulation results from single

element tests. The three tests considered here are a compression test where the volume is changed

and the damage evolves instantaneously to the equilibrium state so that equation 43 is satisfied.

The second test is a shear test where the volume does not change and again the damage is always

instantly updated. The last test is a compression test where the damage evolution is time

dependent, i.e., according to equation 41. In all tests, we chose K = 1 Pa, and μ = 0.75 Pa.

5.1 Instantaneous Damage Evolution, Confined-Compressive Deformation

We compare the results of the simulation on a single element for a confined compression test. The

physical components of the deformation gradient F, and the log-strain tensor lnV imposed are

[F ] =







α 0 0

0 1 0

0 0 1





 , [lnV ] =







ln α 0 0

0 0 0

0 0 0





 . (44)

12
The Jacobian J = detF = α, so that the deviatoric part of the log-strain tensor is

[L

∗

] =







2

3

ln α 0 0

0 −

1

3

ln α 0

0 0 −

1

3

ln α





 . (45)

From J and L

∗ we can solve for the damage, i.e., solving equation 43 using the deformation

given in equation 44. The solution to the damage is

κ = min 1, κ0 −

μ

ξ



K (α ln α − α + 1) + 4 μ ln α

3



(cmin − 1). (46)

This results in the magnitude of the xx-component of the Cauchy Stress to be given by

Txx =



K ln α +

4μ ln α

3α

 min 

1, κ0 −

μ

ξ



K (α ln α − α + 1) + 4μ ln α

3



(cmin − 1)(cmin − 1) + 1.

(47)

The results of the theoretical predictions are shown in figure 1 as solid lines and the symbols are

values taken from simulation results. Panels a and c plot the damage κ as a function of the

strain L

∗

xx, and panels b and d plot the xx-component of stress as a function of the strain L

∗

xx. In

these plots, κ

0 = 0.02 and κ∗ = 1. The top two plots explore the result of holding μ/ξ = 0.1 and

varying cmin from 0.2 to 1, while the bottom two plots explore the result of holding cmin = 0.3

and varying ξ so that μ/ξ = 0.2 to 0.1. In all cases the simulation results agree completely with

the theoretical predictions.

The stress-strain curves in panels b and d show the key feature required for the physical

instability. Curves that exhibit a zero slope in the stress followed by a decrease in stress at

increased strain are candidates for exploring the physical instability. We leave this largely without

discussion in this work since it will be developed further in subsequent work.

13
0 0.5 1

0

0.5

1 a)

κ

L∗

xx

0 0.5 1

0

1

2

3 b)

Txx (Pa)

L∗

xx

0 0.5 1

0

0.5

1 c)

κ

L∗

xx

0 0.5 1

0

1

2

3 d)

Txx (Pa)

L∗

xx

cmin = 0.2

cmin = 0.4

cmin = 0.6

cmin = 0.8

cmin = 1

ξ/μ = 0.2

ξ/μ = 0.4

ξ/μ = 0.6

ξ/μ = 0.8

ξ/μ = 1

Fig. 1 Single element compression test. Comparison of the theoretical prediction (solid lines) against the

simulation results (symbols). Panels a and b hold ξ/μ = 0.1 and vary cmin from 0.2 to 1. Panels c and d

hold cmin = 0.3 and vary ξ/μ from 0.2 to 1. Stress made positive for viewing purposes although its value is

negative in compression..

14
5.2 Instantaneous Damage Evolution, Isochoric Shear Deformation

Here we compare the results of the simulation on a single element for an isochoric shear

deformation in the Y-direction with instantaneous damage evolution. The physical components of

the deformation gradient F and left-stretch tensor V are

[F ] =







1 0 0

0 1 α

0 0 1





 , [V ] =







1 0 0

0

2+α

2

√

4+α2

√ α

4+α2

0 √ α

4+α2

√

2

4+α2





 . (48)

Using MuPad (MathWorks), one can calculate the deviatoric log-strain tensor L

∗

, the damage κ,

and the stress T from the previous expressions for F and V . Since the explicit forms of these

variables are quite complicated, they are not reproduced here. Instead, figure 2 compares the

simulation results (symbols) against the theoretical prediction (solid lines). This figure is very

similar to the previous one with the exception that the component of the stress and strain that is

plotted is the yz-component. The figure clearly shows agreement between theory and simulation

and that the physical instability can manifest itself in shear as well.†

†The level of shear in figure 2 likely exceeds the range of validity for equation 27a.

15
0 0.1 0.2 0.3 0.4

0

0.5

1 a)

κ

L∗

yz

0 0.1 0.2 0.3 0.4

0

0.2

0.4

b)

Tyz (Pa)

L∗

yz

0 0.1 0.2 0.3 0.4

0

0.5

1 c)

κ

L∗

yz

0 0.1 0.2 0.3 0.4

0

0.2

0.4

0.6 d)

Tyz (Pa)

L∗

yz

cmin = 0.2

cmin = 0.4

cmin = 0.6

cmin = 0.8

cmin = 1

ξ/μ = 0.2

ξ/μ = 0.4

ξ/μ = 0.6

ξ/μ = 0.8

ξ/μ = 1

Fig. 2 Single element shear test. Comparison of the theoretical prediction (solid lines) against the simulation

results (symbols). Panels a and b hold ξ/μ = 0.1 and vary cmin from 0.2 to 1. Panels c and d hold cmin =

0.3 and vary ξ/μ from 0.2 to 1..

16
5.3 Dynamic Damage Evolution, Confined-Compressive Deformation

The last verification step presented here considers the time-dependent damage evolution. We

impose a very rapid compression, identical to that in section 5.1 except that the damage is now

evolved dynamically according to equation 41. Here the initial value of the damage is equal to the

equilibrium damage, i.e., κ(0) = κ

0 = 0.02, and the kinetic coefficient K = 0.0001 (1/Pa·s).

For the following analytical calculation, we ignore the initial transient behavior. We assume the

body is initially deformed and that the damage starts off at the damage associated with minimum

energy κ(0) = κ

0

. Under these conditions we find the solution to the damage evolution to be

κ(t) = min 1, κ(0)e−K t ξ − κ

0



e

−K t ξ − 1



−

μ (cmin − 1)

ξ



1 − e

−K t ξ



4

3

μ log α + K (α log α − α + 1)

(49)

with a corresponding magnitude of the xx component of the stress evolution given by:

Txx =



K log α +

4μ log α

3α



φ(κ(t)). (50)

Figure 3 shows good agreement between the simulation and the theoretical prediction. The only

discrepancies arise from the initial transient behavior that we ignored (see the symbols near the

vertical axis panels b and d).

17
0 5000 10000

0

0.5

1 a)

κ

Time

0 5000 10000

0

2

4

6

8 b)

Txx (Pa)

Time

0 5000 10000

0

0.2

0.4

0.6

0.8 c)

κ

Time

0 5000 10000

0

2

4

6

8 d)

Txx (Pa)

Time

cmin = 0.2

cmin = 0.4

cmin = 0.6

cmin = 0.8

cmin = 1

ξ/μ = 0.2

ξ/μ = 0.4

ξ/μ = 0.6

ξ/μ = 0.8

ξ/μ = 1

Fig. 3 Single element dynamic compression test. Comparison of the theoretical prediction (solid lines)

against the simulation results (symbols). Panels a and b hold ξ/μ = 0.1 and vary cmin from 0.2 to 1. Panels c and d hold cmin = 0.3 and vary ξ/μ from 0.2 to 1..

18
6. Example Problems

The previous sections introduced the equations of motion and the verifications of the model’s

implementation. Here, we present some examples of how the model could be used.

6.1 One-Dimensional Dynamic Confined-Compression

In this section we consider a one-dimensional (1-D) problem of a linear alignment of material that

is slowly (compared to the elastic wave speeds) compressed until it becomes unstable. This

problem is investigated by fixing one end of the material at Z = L while the end at Z = 0 is

slowly displaced according to the function:

δ(t) = (

δmax

t

tc

: t ≤ tc

δmax : t > tc

. (51)

This displacement function has a discontinuity in its derivative at 0 and at tc. The dimensions and

material parameters used for this simulation are presented in table 1.

Table 1 Summary of material parameters used in the 1-D stability example..

L (m) ρ (kg/m3) K (Pa) μ (Pa) ξ (Pa) K (1/Pa·s) cmin κ

0 κ∗

tc (s) δmax (m)

1 1e6 10 7.5 1 0.1 0.3 0.02 1 15000 0.29

Figure 4 shows the results of this simulation. Each panel of the figure plots the damage

distribution versus location for a given time. The difference in time between each panel is not the

same. At first, in panels a–e, the damage evolution is similar to that of what we would expect in

the quasi-static case, i.e., it is proportional to the strain applied to the system and shows no

localization of damage. However, in panels f and g the damage has started to accumulate at the

end that is displaced. This localization of damage culminates by panels h and i where the

saturation of damage to κ

∗

in panel i sends off a relaxation wave in the damage. This wave

propagates toward the fixed end (Z = L) where a reflection of the wave at the boundary causes

two additional localizations of damage evolution, panels j and k. At this point the damage

evolution slows down as the system reaches a new stable equilibrium, panels l–o.

19
0 0.5 1

0

0.2

0.4

0.6

0.8

1

Damage

a) 1211 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

b) 6127 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

c) 12350 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

d) 13100 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

e) 13230 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Damage

f) 13350 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

g) 13480 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

h) 13600 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

i) 13720 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

j) 13850 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Damage

Location (m)

k) 13970 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Location (m)

l) 14090 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Location (m)

m) 14220 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Location (m)

n) 17330 s

0 0.5 1

0

0.2

0.4

0.6

0.8

1

Location (m)

o) 24820 s

Fig. 4 A One-Dimensional alignment is slowly compressed from its Z = 0 end. Each panel plots the

damage vs. position within the alignment, where the plots are on the initial configuration. The time for each

panel is shown in the upper right corner and the time difference between each panel is not constant..

One subtle feature of the constitutive model is that the chemical energy term introduces a

chemical potential that drives the damage toward κ

0

. This means that the material can effectively

heal and that the damage is reversible. This feature is present in this example and can be seen by

comparing panels e and o. In panel e, the material is fairly uniformly damaged with a value

slightly larger than 0.6. In panel o, the material is split into fully damaged material κ = 1 and

material that is only partially damaged κ ≈ 0.5. In some cases this self-healing feature may be

undesirable, and a modification to the model is discussed in section 7.1.

6.2 Two-Dimensional Radial Damage Bands

In this section we consider a two-dimensional (2-D) problem of an annular disc that is slowly

(compared to the elastic wave speeds) compressed until it becomes unstable and forms radial

damage bands. This is similar to the problem considered by Grinfeld et al.,5except that the

mechanical system is not evolved quasi-statically. The outer radius Ro = 1 of an annular disc is

20
slowly displaced radially toward the center of the disc according to the piecewise function given

in equation 51. The inner boundary at Ri = 0.1 is stress free and the disc is fixed in the axial

coordinate. This type of deformation concentrates stresses near the inner radius. The simulation

used the material parameters given in table 2. Because of the choices of material parameters and

the imposed deformation, the maximal displacement of the inner radius is only 0.004 m towards

the center.

Table 2 Summary of material parameters used in the 2-D radial damage band instability example..

Ri (m) Ro ρ (kg/m3) K (kPa) μ (kPa) ξ (Pa) K (1/Pa·s) cmin κ

0 κ∗

tc (s) δmax (m)

0.1 1 1e6 10 7.5 1 0.1 0.3 0.02 1 3000 0.05

Each panel in figure 5 shows the damage κ in false-color plotted on the initial configuration

throughout time. The time elapsed between each panel is 175 s. Initially, the damage uniformly

accumulates throughout the material. This can be seen by noting the uniform color distribution in

panels a–f. The solution to this problem without damage is rotationally symmetric and

concentrates energy at the small opening. This gives insight into why the damage begins to localize

near the center and where the material instability is first manifested. The accumulation of damage

causes a local weakening near the opening of the disc and the germination of radial damage bands

(panels g and h)

†

. Six radial damage bands clearly form (panel i) leaving the surrounding material

intact. The local weakening around the damage bands causes a strain concentration ahead of the

existing band. This subsequently causes the instability to propagate and causes the band to grow

predominantly in the radial direction as can be seen in panels h–m. The length of the band grows at

a finite velocity, where one can estimate the speed of the tip of the vertical (90◦) band in panels h

and l, 0.001 m/s. The damage is arrested as the stress is alleviated and the system finds a new stable

configuration. No additional damage accumulates in panels m–o. The damage bands alleviate hoop

stresses that are concentrated at the center opening and that develop from the imposed deformation.

In future developments it may be necessary to introduce a mechanism for hoop damage bands to

develop from radial instabilities. Two possible alterations are discussed in sections 7.1 and 7.2.

†

In these damage bands the material is still intact. In an abstract way these bands can be thought of as a proxy for

radial cracks.

21
0.25

0.5

0.75

kappa_inter

0

1

1

0

κ

0.25

0.5

0.75

kappa_inter

0

1

a)

0.25

0.5

0.75

kappa_inter

0

1

b)

0.25

0.5

0.75

kappa_inter

0

1

c)

0.25

0.5

0.75

kappa_inter

0

1

d)

0.25

0.5

0.75

kappa_inter

0

1

e)

0.25

0.5

0.75

kappa_inter

0

1

f) 0.25

0.5

0.75

kappa_inter

0

1

g) 0.25

0.5

0.75

kappa_inter

0

1

h) 0.25

0.5

0.75

kappa_inter

0

1

i) 0.25

0.5

0.75

kappa_inter

0

1

j) 0.25

0.5

0.75

kappa_inter

0

1

k)

0.25

0.5

0.75

kappa_inter

0

1

l)

0.25

0.5

0.75

kappa_inter

0

1

m)

0.25

0.5

0.75

kappa_inter

0

1

n)

0.25

0.5

0.75

kappa_inter

0

1

o)

Fig. 5 Two-Dimensional instability and the development of radial damage bands. A 2-D annular disc is

slowly compressed from its outside edge r = Ro. Each panel plots the damage in false-color as it varies

with position within the disc. Plots are on the initial configuration. The elapsed time between each panel is

175 s. At first the damage accumulates symmetrically before a material instability produces the localization

of damage along radial lines..

22
7. Future Considerations

In the previous sections we introduced and verified the implementation of a mechanochemical

model for failure. In this section we consider some possible alternatives to the choices of the

damage function φ and the energy.

7.1 Removing Healing From the Mechanochemical Model

In section 6.1 we commented that the quadratic term in the chemical energy can result in a

healing process. This healing process is driven by the chemical potential and may not be

physically realistic. Analogous to Grinfeld and Wright,10 we now restrict the sign of κ ̇ so that

damage is irreversible. Here we consider some ways to modify the current model and

implementation to correct for this feature.

The kinetic evolution of the damage is given by the rate equation 38b. Repeated here for clarity:

−

1

K

∂κ

∂t =

∂φ

∂κ (eelastic(J, L

∗

)) + ξ (κ − κ

o

) . (52)

For our specific choice of φ (and likely in more general cases), the derivative of the damage

function with respect to the damage is negative. Since eelastic > 0, κ ̇ will be negative if

eelastic(J, L

∗

) <

ξκ∗(κ − κ

o

)

1 − cmin

. (53)

This would certainly become an issue in a test where the loading is applied very rapidly and there

are elastic waves propagating back and forth. A wave propagating past would initially cause

damage, but as the wave continued on, the damaged material might be in a lower elastic energy

state and begin to self-heal (recall the results from section 6.1).

We can quite trivially remove the self-healing feature. The algorithm outlined in section 4

calculates the mechanical energy associated with the deformation. This energy could then used to

evaluate the inequality in equation 53. If the inequality holds true, then applying the traditional

update equation 41 would result in a smaller damage value in the next time step. However, instead

of using the update equation, one could set ∂κ/∂t = 0, and

κ

n+1 = κn

. (54)

23
This type of modification might be sufficient to observe radial and hoop damage bands a dynamic

test where healing would have been an issue before.

7.2 Decoupling Bulk Damage and Deviatoric Damage Functions

Here, we consider a simple extension to the current model where the only modification is that

there are two damage functions that can affect the volumetric and isochoric responses separately:

e = φ1(κ)evol(J) + φ2(κ)eiso(L

∗

) + echem(κ). (55)

The stress is no longer given by equation 42, but instead

T = φ1(κ)K ln JI + φ2(κ)

2μ

J

L

∗

. (56)

The damage still evolves in time; however, derivatives of φ1 and φ2 with respect to κ will enter

into the kinetic equation separately:

−

1

K

∂κ

∂t =

∂φ1

∂κ evol(J) + ∂φ2∂κ eiso(L

∗

) + ξ (κ − κ

o

) . (57)

This would be particularly useful if the material was more sensitive to damage in shear than in its

bulk response. One possible choice of damage functions is

φ1(κ) = 1 − (1 − cmin)

κ

κ

∗

and (58)

φ2(κ) = max h0, 1 − α

κ

κ

∗

i

. (59)

A dimensionless parameter α > 0 sets how rapidly the damage function φ2 drops to zero.∗ The

maximum operator is needed here since the original notion of κ

∗

is lost by introducing α which

would otherwise make φ2 negative when ακ > κ∗. In this example cmin is dropped from φ2 so

that the material completely loses its deviatoric response and behaves more like a soft, nearly

incompressible material when it is maximally damaged.

An alternate approach to this type of modification would be to introduce a second damage

parameter κ2, with its own damage evolution and maximal damage κ

∗

2

. In this case there would be

two internal state variables and two damage evolution equations κ ̇ ito solve. This also would

require the notion that damage carries a direction instead of simply being a scalar.

∗There is still no restriction on the sign of the time derivative of the damage. Thus, bulk and deviatoric damage can

still heal and that α merely adjusts the rate of damage accumulation to be preferentially bulk or preferentially shear.

24
7.3 Damage Growth Asymmetry in Compression and Tension

We introduce a pathway to add an asymmetric growth of the damage in tension and compression.

This type of model might be useful to ensure that a material does not lose strength in compression

and to promote damage growth in tension. The modified energy density is

e = φ(κ) [αtH(J − 1) + αcH(1 − J)] evol(J) + φ(κ)eiso(L

∗

) + echem(κ), (60)

where H is a Heaviside function:

H(x) = (

1 : x ≥ 0

0 : x < 0

. (61)

The construction in the square brackets in equation 60 acts as a switch between the two energy

scale factors αt and αc. H(J − 1) takes on the value 1 when the material is in tension, and if

αt > αc, the energy is increased by αt − αc. Assuming this form of the energy density gives the

desired behavior in the kinetic law

−

1

K

∂κ

∂t =

∂φ

∂κ [αtH(J − 1) + αcH(1 − J)] evol(J) + ∂φ∂κeiso(L

∗

) + ξ (κ − κ

o

) . (62)

The deviatoric stress is

J T

∗ =

∂e

∂L

∗









J

= 2μ L

∗

, (63)

and the pressure can be obtained from applying equation 25a to equation 60,

p = −

∂e

∂J









L∗

= −

∂

∂J ([αtH(J − 1) + αcH(1 − J)] K (J ln J − J + 1)) . (64)

Evaluating the derivatives and simplifying gives:

p = [αtH(J − 1) + αcH(1 − J)] K ln J . (65)

Now the slope of the pressure has a discontinuity as we cross J = 1 from above or below:

∂p

∂J = [αtH(J − 1) + αcH(1 − J)] K

1

J

+ [αt δ(J − 1) − αc δ(1 − J)] K ln J , (66)

where

lim

J→1−

∂p

∂J = αc K (67)

25
and

lim

J→1+

∂p

∂J = αt K . (68)

This is only a minor concern since the derivative of the stress response already has a discontinuity

from the damage function (see figures 1 and 2).

8. Conclusions

In conclusion we have developed a mechanochemical model for dynamic failure using a finite

deformation formulation analogous to the model developed by Grinfeld and Wright.1 We have

conducted single element tests to verify the model is working. We have also outlined some

example problems that the model was implemented to study. We also discussed some of the

limitations and possible future alterations that could be made to the model.

26
9. References

1. Grinfeld, M. A.; Wright, T. W. Thermodynamics of Solids: Recent Progress With

Applications to Brittle Fracture and Nanotechnology. In 23rd U.S. Army Science Conference,

Orlando, FL, 2002.

2. Grinfeld, M. A. Instability of the Interface Between a Nonhydrostatically Stressed Elastic

Body and a Melt. Akademiia Nauk SSSR Doklady 1986, 290, 1358-1363.

3. Grinfeld, M. A. Thermodynamic Methods in the Theory of Heterogeneus

Systems; Longman Scientific and Technical: New York, 1991.

4. Kassner, K.; Misbah, C.; Müller, J.; Kappey, J.; Kohlert, P. Phase-Field Modeling of

Stress-Induced Instabilities. Physical Review E. 2001, 63 (3), 036117.

5. Grinfeld, M. A.; McCauley, J. W.; Schoenfeld, S. E.; Wright, T. W. Failure Pattern

Formation in Brittle Ceramics and Glasses. In 23rd International Symposium on Ballistics,

April 2007, Tarragona, Spain, 2007.

6. Kachanov, L. M. Introduction to Continuum Damage Mechanics; Springer: New York, 1986.

7. Chaboche, J. L. Continuum Damage Mechanics: Part I General Concepts. Journal of

Applied Mechanics 1988, 55 (1), 59–64.

8. Chaboche, J. L. Continuum Damage Mechanics: Part II Damage Growth, Crack Initiation,

and Crack Growth. Journal of Applied Mechanics 1988, 55 (1), 65–72.

9. Holzapfel, G. A. Nonlinear Solid Mechanics: A Continuum Approach for

Engineering; John Wiley & Sons Ltd.: New York, 2000.

10. Grinfeld, M. A.; Wright, T. W. Thermodynamics of Brittle Fracture; ARL-TR-3659; U.S.

Army Research Laboratory: Aberdeen Proving Ground, MD, 2005.

11. Truesdell, C.; Noll, W. The Non-Linear Field Theories of Mechanics; Springer: New York,

2004.

12. Hencky, H. The Law of Elasticity for Isotropic and Quasi-Isotropic Substances by Finite

Deformations. Journal of Rheology 1931, 2, 169.

27
13. Anand, L. On H. Hencky’s Approximate Strain-Energy Function for Moderate

Deformations. Journal of Applied Mechanics 1979, 46, 78.

14. Anand, L. Moderate Deformations in Extension-Torsion of Incompressible Isotropic Elastic

Materials. Journal of the Mechanics and Physics of Solids 1986, 34 (3), 293–304.

15. Scheidler, M. On the Coupling of Pressure and Deviatoric Stress in Hyperelastic Materials.

In Proceedings of the 13th Army Symposium on Solid Mechanics, Aug 1993, Plymouth, MA,

1993.

16. Scheidler, M. Time Rates of Generalized Strain Tensors Part I: Component Formulas.

Mechanics of Materials 1991, 11 (3), 199–210.

28
Appendix A. Analytical Solution to Kinetic Equation

29
In the main text, we derived a kinetic law relating the evolution of the damage to the elastic

energy of the material,

∂κ

∂t = −K



∂φ

∂κeelastic + ξ (κ − κ

o

)



. (A-1)

For the special case discussed in the main text, equation A-1 readily integrates. This appendix

considers a more general approach that could be used for other damage functions φ.

For notational simplicity we write,

∂φ

∂κ = g(κ(t)). (A-2)

and rewrite equation A-1 as

dκ

dt

= −K (g(κ(t)) eelastic(t) + ξ (κ − κ

o

)) . (A-3)

Recall that the elastic energy is the hypothetical elastic energy of an undamageable material, i.e.,

the energy that only depends on the displacements and the elastic moduli. In general, g might not

be a simple function of the damage so a method of Laplace transforms is useful when attempting

to solve this differential equation. In the following, L[f(t)] = F(s) denotes the Laplace

transform of f(t). Similarly, L

−1

[F(s)] = f(t) denotes the inverse Laplace transform of F(s).

Applying the Laplace transform to both sides of equation A-3 gives

sL[κ(t)] − κ(0) = −KL[g(κ(t)) eelastic(t)] − KξL[κ(t)] + Kξκo

s

, (A-4)

which is an algebraic equation that can be solved for L[κ(t)],

L[κ(t)] = κ(0)

s + Kξ +

1

Kξ 

Kξκo

s

−

Kξκo

s + Kξ 

−

KL[g(κ(t)) eelastic(t)]

s + Kξ . (A-5)

Applying the inverse Laplace transform gives:

κ(t) = L

−1



κ(0)

s + Kξ 

+ L

−1



κ

o

s



− L−1



κ

o

s + Kξ 

− L−1



KL[g(κ(t)) eelastic(t)]

s + Kξ 

, (A-6)

κ(t) = κ(0)e

−Kξt − κo



e

−Kξt − 1



− KL

−1



L[g(κ(t)) eelastic(t)]L



e

−Kξt . (A-7)

The last term is the inverse Laplace transform of the product of two Laplace transforms, which is

the convolution of the two original functions:

κ(t) = κ(0)e

−Kξt − κo



e

−Kξt − 1



− K

Z t

0

∂φ

∂κ(τ ) eelastic(τ )e

−Kξ(t−τ)

dτ . (A-8)

30
Since, in general, ∂φ/∂κ will depend on κ, this is an integral equation for κ(t). It can be used to

update the damage parameter for any damage function as a convolution of ∂φ/∂κ and eelastic.

Considering a small time step ∆t, however, one can derive:

κ(t + ∆t) = κ(t)e

−Kξ∆t − κo



e

−Kξ∆t − 1



− K

Z t+∆t

t

∂φ

∂κ(τ ) eelastic(τ )e

−Kξ(t+∆t−τ)

dτ . (A-9)

If ∆t is chosen for stability in the explicit codes, the change in elastic energy will likely be small

over the time interval t → t + ∆t. Thus, we take the elastic energy outside of the convolution

integral.∗ We also assume that the damage is not changing rapidly during this time interval so that

∂φ/∂κ is slowly varying and can also be taken outside the integral. In the special case of φ given

by equation 5 this holds trivially, since ∂φ/∂κ is constant. Thus, we can write the incremental

change in the damage as:

κ(t + ∆t) ≈ κ(t)e

−Kξ∆t − κo



e

−Kξ∆t − 1



−

1

ξ



1 − e

−Kξ∆t

 ∂φ

∂κ eelastic (A-10)

Thus, the updated damage κ(t + ∆t) depends on the current damage κ(t) the equilibrium damage

κ

o

and the elastic energy.

Using a superscript n denotes that the variable is evaluated at the current time step, and n + 1 the

future time step gives

κ

n+1 = κn

e

−Kξ∆t − κo



e

−Kξ∆t − 1



−

1

ξ



1 − e

−Kξ∆t





∂φ

∂κn

e

n+1

elastic . (A-11)

When φ is given by equation 5,

∂φ

∂κ =

1 − cmin

κ

∗

, (A-12)

and

κ

n+1 = κn

e

−Kξ∆t − κo



e

−Kξ∆t − 1



+

1

ξ



1 − e

−Kξ∆t

 1 − cmin

κ

∗

e

n+1

elastic . (A-13)

In the text, the elastic energy term is used in two different ways. First, that of a linear elastic

material and later, a hyper-elastic material.

∗

In future models where ∂φ/∂κ is not constant, additional state variables may be needed in the constitutive model,

which keeps track of past elastic energies and past damage values. These might be used in some Simpson quadrature

rule to simplify the convolution so that it could be evaluated numerically.

31
INTENTIONALLY LEFT BLANK.

32
Appendix B. SIERRA Implementation

33
1 # i n c l u d e <m o del s / B r i t t l e E l a s t i c . h>

2 # i n c l u d e <cmath >

3 # i n c l u d e < k i n e m a t i c s / l am e _ g e n e r al _ m et h o d s . h>

4 # i n c l u d e < k i n e m a t i c s / LameUtil . h>

5 # i n c l u d e < k i n e m a t i c s / Ki n e m ati c s . h>

6 u s i n g name s pace std ;

7 name s pace lame {

8

9 Material * BrittleElastic : : createMaterial( c o n s t MatProps & props ) {

10 r e t u r n new BrittleElastic( props ) ;

11 }

12 / / ******************************************************************

13 / / T h e r e a r e t h r e e s t a t e v a r i a b l e s ( a l t h o u g h o nl y one c h a n g e s r i g h t now )

14 / / The f i r s t two a r e b a si c , BULK_MODULUS and SHEAR_MODULUS

15 / / The t h i r d which e v o l v e s i n tim e i s t h e damage p a r a m e t e r KAPPA

16 / / To make t h i s model work , howeve r , you nee d t o s p e c i f y

17 / / XI_CHEMICAL damage c h e mi c al c o n s t a n t , K_KINETIC , KAPPA_ZERO t h e s t e a d y s t a t e damage

18 / / Al s o w i l l nee d C_MIN and KAPPA_STAR which a p p e a r i n t h e damage f u n c t i o n p hi .

19 / / T he se a r e a l l c o n s t a n t t h r o u g h o u t t h e r u n and a r e p r o p e r t i e s .

20 / / T hi s i s t h e c o n s t r u c t o r f o r t h e B r i t t l e E l a s t i c model . The p r o p e r t i e s

21 / / i t r e a d s i n a r e t h e BULK_MODULUS and SHEAR_MODULUS. T he se

22 / / a r e s t o r e d i n t h e p r o p e r t i e s a r r a y .

23 / /

24 / / ******************************************************************

25 BrittleElastic : : BrittleElastic( c o n s t MatProps & props) :

26 Material(props) {

27 setFlag(USE_LEFT_STRETCH) ;

28 setHyper ( ) ;

29 / /

30 / / M a t e r i a l P r o p e r t y D e f i n i t i o n s

31 / / What i s r e a d i n f rom t h e i n p u t dec

32 / / YOUNGS_MODULUS, POISSONS_RATIO , XI_CHEMICAL, K_KINETIC , KAPPA_ZERO, C_MIN, KAPPA_STAR, ←-

INSTANT − 8 t o t a l p r o p e r t i e s

33

34 mat_name = "BRITTLEELASTIC" ;

35 num_material_properties = 8 ;

36 initializeProperties ( ) ;

37

38 setMaterialProperty( 0 , "YOUNGS_MODULUS" ,props) ;

39 setMaterialProperty( 1 , "POISSONS_RATIO" ,props) ;

40 / / The d e f a u l t b e h a v i o r i s l i n e a r e l a s t i c w i t h o u t damage

41 setMaterialPropertyDefault ( 2 , 0 . 0 ) ;

42 setMaterialProperty( 2 , "XI_CHEMICAL" ,props) ;

43 setMaterialPropertyDefault ( 3 , 0 . 0 ) ;

44 setMaterialProperty( 3 , "K_KINETIC" ,props) ;

45 setMaterialPropertyDefault ( 4 , 0 . 0 ) ;

46 setMaterialProperty( 4 , "KAPPA_ZERO" ,props) ;

47 setMaterialPropertyDefault ( 5 , 0 . 0 ) ;

48 setMaterialProperty( 5 , "C_MIN" ,props) ;

49 setMaterialPropertyDefault ( 6 , 1 . 0 ) ;

50 setMaterialProperty( 6 , "KAPPA_STAR" ,props) ;

51 setMaterialPropertyDefault ( 7 , 0 ) ;

52 setMaterialProperty( 7 , "INSTANTANEOUS" ,props) ;

53

34
54 LAME_FORTRAN(elastic_property_check) ( num_material_properties ,

55 &properties [ 0 ] ) ;

56

57 / /

58 / / S t a t e V a r i a b l e D e f i n i t i o n s

59 / / We a r e u s i n g 5 s t a t e v a r i a b l e s . The b ul k and s h e a r a r e two . t h e n damage ka p pa . t h e n two ←-

e n e r g i e s

60 / /

61

62 num_state_vars = 5 ;

63

64 set_state_variable_alias( "YOUNGS_MODULUS" , 0 ) ;

65 set_state_variable_alias( "POISSONS_RATIO" , 1 ) ;

66 set_state_variable_alias( "KAPPA_INTERNAL" , 2 ) ;

67 set_state_variable_alias( "MECHANICAL_ENERGY" , 3 ) ;

68 set_state_variable_alias( "CHEMICAL_ENERGY" , 4 ) ;

69

70 }

71

72 i n t BrittleElastic : : initialize( matParams * p ) {

73 c o n s t d o u bl e ym = properties [ 0 ] ;

74 c o n s t d o u bl e nu = properties [ 1 ] ;

75

76 originalYoungs = properties [ 0 ] ;

77 originalNu = properties [ 1 ] ;

78

79 xi = properties [ 2 ] ;

80 bigk = properties [ 3 ] ;

81 kzero = properties [ 4 ] ;

82 cmin = properties [ 5 ] ;

83 kstar = properties [ 6 ] ;

84 instant = properties [ 7 ] ;

85

86 c o n s t d o u bl e kappa = properties [ 4 ] ; / / i n i t i a l l y we a r e n e u t r a l / s t e a d y s t a t e damaged

87 d o u bl e * stateOld = p−>state_old ;

88 d o u bl e * stateNew = p−>state_new ;

89 f o r ( i n t i ( 0 ) ; i < p−>nelements ; ++i) {

90 / /

91 / / S et t h e d e f a u l t s t a t e t o t h e bl o c k e l a s t i c c o n s t a n t s

92 / /

93 stateOld [ 0 ] = stateNew [ 0 ] = ym ;

94 stateOld [ 1 ] = stateNew [ 1 ] = nu ;

95 stateOld [ 2 ] = stateNew [ 2 ] = kappa ;

96 stateOld [ 3 ] = stateNew [ 3 ] = 0 . 0 ;

97 stateOld [ 4 ] = stateNew [ 4 ] = 0 . 0 ;

98 stateOld += 5 ;

99 stateNew += 5 ; / / 5 i n t e r n a l s t a t e v a r i a b l e s

100 }

101

102 r e t u r n 0 ;

103 }

104

105 / / ******************************************************************

106 / /

107 / / The g e t S t r e s s method r e t u r n s t h e s t r e s s gi v e n t h e s t r a i n

35
108 / / and t h e tim e s t e p .

109 / /

110 / / ******************************************************************

111

112 i n t BrittleElastic : : getStress( matParams * p ) {

113 c o n s t d o u bl e dt = p−>dt ;

114 / / c o n s t d o u bl e * s t r a i n R a t e = p−> s t r a i n _ r a t e ;

115 d o u bl e * rotation = p−>rotation ;

116 c o n s t d o u bl e * leftStretch = p−>left_stretch ;

117 c o n s t d o u bl e * stressOld = p−> stress_old ;

118 c o n s t d o u bl e * stateOld = p−> state_old ;

119 d o u bl e * stressNew = p−> stress_new ;

120 d o u bl e * stateNew = p−> state_new ;

121

122 c o n s t d o u bl e onehalf = 1 . 0 / 2 . 0 ;

123 c o n s t d o u bl e onethird = 1 . 0 / 3 . 0 ;

124

125 d o u bl e devlnV [ 6 ] ;

126 d o u bl e evals [ 3 ] ;

127 d o u bl e evecs [ 9 ] ;

128 d o u bl e LQT [ 9 ] ;

129 d o u bl e resultmat [ 9 ] ;

130 d o u bl e tracelnV ;

131

132 f o r ( i n t i ( 0 ) ; i < p−>nelements ; ++i) {

133

134 c o n s t d o u bl e ym = originalYoungs ; / / s t a t e O l d [ 0 ] < 0 ? 0 : s t a t e O l d [ 0 ] ; / / ym must be ←-

p o s i t i v e

135 c o n s t d o u bl e pr = originalNu ; / / s t a t e O l d [ 1 ] < −0.9999999 ? −0.9999999 :

136 c o n s t d o u bl e bulk = ym / 3 . 0 / ( 1 . 0 − 2 . 0 * pr) ;

137 c o n s t d o u bl e shear = ym / 2 . 0 / ( 1 . 0 + pr) ;

138 c o n s t d o u bl e kappa = stateOld [ 2 ] > kstar ? kstar :

139 stateOld [ 2 ] < 0 ? 0 : stateOld [ 2 ] ; / / kee p ka p pa l e s s t h a n k s t a r and g r e a t e r t h a n 0

140

141 c o n s t d o u bl e J = leftStretch[xx ] * ( leftStretch[yy]*leftStretch[zz]−leftStretch[yz]*←-

leftStretch[yz ] )

142 −leftStretch[xy ] * ( leftStretch[xy]*leftStretch[zz]−leftStretch[zx]*leftStretch[yz ] )

143 +leftStretch[zx ] * ( leftStretch[xy]*leftStretch[yz]−leftStretch[yy]*leftStretch[zx ] ) ;

144 c o n s t d o u bl e lnJ = FastLog(J) ;

145 c o n s t d o u bl e sheartwooverJ = shear * 2 . 0 / J ;

146 c o n s t d o u bl e lnJdivJ = lnJ / J ;

147

148 devlnV[xx]=leftStretch[xx ] ;

149 devlnV[yy]=leftStretch[yy ] ;

150 devlnV[zz]=leftStretch[zz ] ;

151 devlnV[xy]=leftStretch[xy ] ;

152 devlnV[yz]=leftStretch[yz ] ;

153 devlnV[zx]=leftStretch[zx ] ;

154

155 LAME_FORTRAN(lame_eigen) ( 1 , 1 ,devlnV , evals , evecs) ; / / c a l c u l a t e e i g e n v a l u e s and v e c t o r s

156 / / Then u s e t h e s e v a l u e s t o c a l c u l a t e logV i n i t s e i g e n b a s i s .

157 / / t h e n c o n v e r t b a s i s bac k t o whe re we a r e u s i n g : QλQT = V

158 evals[ 0 ] =FastLog(evals [ 0 ] ) ;

159 evals[ 1 ] =FastLog(evals [ 1 ] ) ;

160 evals[ 2 ] =FastLog(evals [ 2 ] ) ;

36
161

162 / / make Lambda*Q^T , name i s LQT .

163 / / t h i s i s t a k i n g t h e rows o f Q^T and m ult by l am b d a _i f o r e a c h row

164 LQT[fxx] = evals[ 0 ] *evecs[fxx ] ; / /

165 LQT[fxy] = evals[ 0 ] *evecs[fyx ] ; / / f y x b e c a u s e i t i s Q^T , n ot Q

166 LQT[fxz] = evals[ 0 ] *evecs[fzx ] ;

167

168 LQT[fyy] = evals[ 1 ] *evecs[fyy ] ;

169 LQT[fyz] = evals[ 1 ] *evecs[fzy ] ;

170 LQT[fyx] = evals[ 1 ] *evecs[fxy ] ;

171

172 LQT[fzz] = evals[ 2 ] *evecs[fzz ] ;

173 LQT[fzx] = evals[ 2 ] *evecs[fxz ] ;

174 LQT[fzy] = evals[ 2 ] *evecs[fyz ] ;

175

176 / / now Q m a t r i x i s j u s t t h e e v e c s m at ri x , s o m u l t i p l y

177 / / lnV = e v e c s * LQT

178 / / t h i s lame c a l l t o f o r t r a n d o e s n ot work , s o I r e p l a c e d i t wit h a d i r e c t c a l c u l a t i o n i n s t e a d

179 / / LAME_FORTRAN( l a m e _ t e n s o r _ p r o d u c t 3 6 3 6 ) ( 1 , 1 , e v e c s ,LQT, r e s u l t m a t ) ; / / m u l t i p l y e v e c s *LQT

180 / / r e s u l t mat s h o ul d be s ymm et ri c s i n c e i t i s now l n V

181

182 resultmat[fxx] = evecs[fxx]*LQT[fxx]+evecs[fxy]*LQT[fyx]+evecs[fxz]*LQT[fzx ] ;

183 resultmat[fxy] = evecs[fxx]*LQT[fxy]+evecs[fxy]*LQT[fyy]+evecs[fxz]*LQT[fzy ] ;

184 resultmat[fxz] = evecs[fxx]*LQT[fxz]+evecs[fxy]*LQT[fyz]+evecs[fxz]*LQT[fzz ] ;

185

186 resultmat[fyx] = evecs[fyx]*LQT[fxx]+evecs[fyy]*LQT[fyx]+evecs[fyz]*LQT[fzx ] ;

187 resultmat[fyy] = evecs[fyx]*LQT[fxy]+evecs[fyy]*LQT[fyy]+evecs[fyz]*LQT[fzy ] ;

188 resultmat[fyz] = evecs[fyx]*LQT[fxz]+evecs[fyy]*LQT[fyz]+evecs[fyz]*LQT[fzz ] ;

189

190 resultmat[fzx] = evecs[fzx]*LQT[fxx]+evecs[fzy]*LQT[fyx]+evecs[fzz]*LQT[fzx ] ;

191 resultmat[fzy] = evecs[fzx]*LQT[fxy]+evecs[fzy]*LQT[fyy]+evecs[fzz]*LQT[fzy ] ;

192 resultmat[fzz] = evecs[fzx]*LQT[fxz]+evecs[fzy]*LQT[fyz]+evecs[fzz]*LQT[fzz ] ;

193

194 tracelnV = resultmat[fxx]+resultmat[fyy]+resultmat[fzz ] ;

195

196 devlnV[xx]=resultmat[fxx]−onethird*tracelnV ;

197 devlnV[yy]=resultmat[fyy]−onethird*tracelnV ;

198 devlnV[zz]=resultmat[fzz]−onethird*tracelnV ;

199 devlnV[xy]=resultmat[fxy ] ;

200 devlnV[yz]=resultmat[fyz ] ;

201 devlnV[zx]=resultmat[fzx ] ;

202

203 / / u p d at e t h e m e c h a ni c al e n e r g y

204 stateNew[ 3 ] =bulk*(J*lnJ−J+ 1 ) +shear*(devlnV[xx]*devlnV[xx]+devlnV[yy]*devlnV[yy]+devlnV[zz]*←-

devlnV[zz]+

205 2 . 0 *devlnV[xy]*devlnV[xy] + 2 . 0 *devlnV[yz]*devlnV[yz] + 2 . 0 *devlnV[zx]*devlnV[zx ] ) ;

206

207 / / u p d at e t h e new damage

208 c o n s t d o u bl e expbxt = (instant == 1 ) ? 0 : exp(−bigk*xi*dt) ;

209 d o u bl e kappaNew = kappa*expbxt−kzero*(expbxt−1.0 ) −(1.0−expbxt) *(cmin−1.0 ) / kstar*stateNew [ 3 ] / ←-

xi ; / / alphaNew h a s a mu i n i t

210 kappaNew = kappaNew > kstar ? kstar :

211 kappaNew < 0 . 0 ? 0 . 0 : kappaNew ;

212

213 c o n s t d o u bl e phi = 1 . 0 − ( 1 . 0 − cmin) *kappaNew / kstar ;

37
214

215 / / f i n a l l y c a l c u l a t e t h e ca uc h y s t r e s s

216 stressNew [ 0 ] = phi*bulk*lnJ + phi*sheartwooverJ*devlnV[xx ] ; / / xx

217 stressNew [ 1 ] = phi*bulk*lnJ + phi*sheartwooverJ*devlnV[yy ] ; / / yy

218 stressNew [ 2 ] = phi*bulk*lnJ + phi*sheartwooverJ*devlnV[zz ] ; / / zz

219 stressNew [ 3 ] = phi*sheartwooverJ*devlnV[xy ] ; / / xy

220 stressNew [ 4 ] = phi*sheartwooverJ*devlnV[yz ] ; / / yz

221 stressNew [ 5 ] = phi*sheartwooverJ*devlnV[zx ] ; / / zx

222

223 / / r o t a t e ( o r u n r o t a t e ) bac k t o r e f e r e n c e RT T R

224 / / i n S i e r r a T i s c a l l e d si gma

225 / / u s e r e s u l t m a t t o t e m p o r a r i l y s t o r e t h e i n f o

226 unRotate(stressNew , rotation) ;

227

228 stateNew [ 0 ] = ym*phi ; / / o r i g i n a l y o u n g s modulus i s r e d u c e d by p hi

229 stateNew [ 1 ] = pr ; / / p o i s s o n s r a t i o i s n ot a f f e c t e d by damage

230 stateNew [ 2 ] = kappaNew ;

231 stateNew[ 3 ] =stateNew[ 3 ] *phi ;

232 stateNew [ 4 ] = onehalf * (kappaNew − kzero) *(kappaNew − kzero) *xi ;

233 rotation +=9;

234 stressOld += 6 ;

235 stressNew += 6 ;

236 leftStretch +=6;

237 stateOld += 5 ;

238 stateNew +=5;

239 }

240 r e t u r n 0 ;

241 }

242 / /−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

243 i n t BrittleElastic : : getInitialElasticModuli( matParams * p ) {

244 c o n s t d o u bl e * stateOld = p−> state_old ;

245 f o r ( i n t i ( 0 ) ; i < p−>nelements ; ++i) {

246 p−>pcVars_old[PC_YOUNGS_MODULUS ] [ i] = stateOld [ 0 ] ; / / 9 . 0 * s t a t e O l d [ 0 ] * s t a t e O l d [ 1 ] / ( 3 . 0 * ←-

s t a t e O l d [ 0 ] + s t a t e O l d [ 1 ] ) ;

247 p−>pcVars_old[PC_POISSONS_RATIO ] [ i] = stateOld [ 1 ] ; / / ( 3 . 0 * s t a t e O l d [ 0 ] −2. 0* s t a t e O l d [ 1 ] ) ←-

/ ( 2 . 0 * ( 3 . 0 * s t a t e O l d [ 0 ] + s t a t e O l d [ 1 ] ) ) ;

248 }

249 r e t u r n 0 ;

250 }

251

252 i n t BrittleElastic : : pcElasticModuli( matParams * p )

253 {

254 c o n s t d o u bl e * stateOld = p−> state_old ;

255 d o u bl e * stateNew = p−> state_new ;

256 d o u bl e * oldE = p−>pcVars_old[PC_YOUNGS_MODULUS ] ;

257 d o u bl e * newE = p−>pcVars_new[PC_YOUNGS_MODULUS ] ;

258 d o u bl e * oldNu = p−>pcVars_old[PC_POISSONS_RATIO ] ;

259 d o u bl e * newNu = p−>pcVars_new[PC_POISSONS_RATIO ] ;

260

261 f o r ( i n t i ( 0 ) ; i < p−>nelements ; ++i)

262 {

263 oldE[i] = stateOld [ 0 ] ;

264 newE[i] = stateNew [ 0 ] ;

265 oldNu[i] = stateOld [ 1 ] ;

266 newNu[i] = stateNew [ 1 ] ;

38
267 stateOld += 5 ;

268 stateNew +=5;

269 }

270 r e t u r n 0 ;

271 }

272 v oi d BrittleElastic : : unRotate( d o u bl e * stressNew , d o u bl e * rotation)

273 {

274 d o u bl e resultmat [ 6 ] ;

275 / / RT SR

276 resultmat[xx] = rotation[fxx ] * ( rotation[fxx]*stressNew[xx] + rotation[fyx]*stressNew[xy←-

] + rotation[fzx]*stressNew[zx ] ) + rotation[fyx ] * ( rotation[fxx]*stressNew[xy] + ←-

rotation[fyx]*stressNew[yy] +

277 rotation[fzx]*stressNew[yz ] ) + rotation[fzx ] * ( rotation[fyx]*stressNew[yz] + rotation[fxx]*←-

stressNew[zx] + rotation[fzx]*stressNew[zz ] ) ;

278 resultmat[yy] = rotation[fxy ] * ( rotation[fxy]*stressNew[xx] + rotation[fyy]*stressNew[xy←-

] + rotation[fzy]*stressNew[zx ] ) + rotation[fyy ] * ( rotation[fxy]*stressNew[xy] + ←-

rotation[fyy]*stressNew[yy] +

279 rotation[fzy]*stressNew[yz ] ) + rotation[fzy ] * ( rotation[fyy]*stressNew[yz] + rotation[fxy]*←-

stressNew[zx] + rotation[fzy]*stressNew[zz ] ) ;

280 resultmat[zz] = rotation[fxz ] * ( rotation[fxz]*stressNew[xx] + rotation[fyz]*stressNew[xy←-

] + rotation[fzz]*stressNew[zx ] ) + rotation[fyz ] * ( rotation[fxz]*stressNew[xy] + ←-

rotation[fyz]*stressNew[yy] +

281 rotation[fzz]*stressNew[yz ] ) + rotation[fzz ] * ( rotation[fyz]*stressNew[yz] + rotation[fxz]*←-

stressNew[zx] + rotation[fzz]*stressNew[zz ] ) ;

282

283 resultmat[xy] = rotation[fxx ] * ( rotation[fxy]*stressNew[xx] + rotation[fyy]*stressNew[xy←-

] + rotation[fzy]*stressNew[zx ] ) + rotation[fyx ] * ( rotation[fxy]*stressNew[xy] + ←-

rotation[fyy]*stressNew[yy] +

284 rotation[fzy]*stressNew[yz ] ) + rotation[fzx ] * ( rotation[fyy]*stressNew[yz] + rotation[fxy]*←-

stressNew[zx] + rotation[fzy]*stressNew[zz ] ) ;

285 resultmat[yz] = rotation[fxy ] * ( rotation[fxz]*stressNew[xx] + rotation[fyz]*stressNew[xy←-

] + rotation[fzz]*stressNew[zx ] ) + rotation[fyy ] * ( rotation[fxz]*stressNew[xy] + ←-

rotation[fyz]*stressNew[yy] +

286 rotation[fzz]*stressNew[yz ] ) + rotation[fzy ] * ( rotation[fyz]*stressNew[yz] + rotation[fxz]*←-

stressNew[zx] + rotation[fzz]*stressNew[zz ] ) ;

287 resultmat[zx] = rotation[fxz ] * ( rotation[fxx]*stressNew[xx] + rotation[fyx]*stressNew[xy←-

] + rotation[fzx]*stressNew[zx ] ) + rotation[fyz ] * ( rotation[fxx]*stressNew[xy] + ←-

rotation[fyx]*stressNew[yy] +

288 rotation[fzx]*stressNew[yz ] ) + rotation[fzz ] * ( rotation[fyx]*stressNew[yz] + rotation[fxx]*←-

stressNew[zx] + rotation[fzx]*stressNew[zz ] ) ;

289

290 f o r ( u n si g n e d i n t j= 0; j<6;++j)

291 stressNew[j]=resultmat[j ] ;

292 }

293 } / / lame

39
INTENTIONALLY LEFT BLANK.

40
1

(PDF)

DEFENSE TECHNICAL

INFORMATION CTR

DTIC OCA

2

(PDF)

DIRECTOR

US ARMY RESEARCH LAB

RDRL CIO LL

IMAL HRA MAIL & RECORDS MGMT

1

(PDF)

GOVT PRINTG OFC

A MALHOTRA

4

(PDF)

DIR USARL

RDRL DPW

R COATES

P FROUNFELKER

R SPINK

M TEGTMEYER

1

(PDF)

DIR USARL

RDRL ROP L

F GREGORY

2

(PDF)

DIR USARL

RDRL CIH C

B HENZ

M VINDIOLA

10

(PDF)

DIR USARL

RDRL HR

P FRANASZCZUK

RDRL HRS C

S GORDON

W HAIRSTON

B LANCE

K MCDOWELL

J MCARDLE

K OIE

A PASSARO

M PETERSON

J VETTEL

7

(PDF)

DIR USARL

RDRL SLB W

D BOOTHE

A BREUER

N EBERIUS

P GILLICH

C KENNEDY

A KULAGA

W MERMAGEN

44

(PDF)

DIR USARL

RDRL WM

S KARNA

RDRL WML C

T PIEHLER

RDRL WML H

B SCHUSTER

RDRL WMM B

B LOVE

RDRL WMM G

L PIEHLER

N ZANDER

RDRL WMP

S SCHOENFELD

RDRL WMP B

A DAGRO

A DILEONARDI

A DWIVEDI

W EVANS

C GUNNARSSON

C HOPPEL

Y HUANG

M LYNCH

J MCDONALD

P MCKEE

B SANBORN

S SATAPATHY

A SOKOLOW

C WEAVER

T WEERASOORIYA

S WOZNIAK

T ZHANG

K ZIEGLER

RDRL WMP C

S BILYK

T BJERKE

D CASEM

J CLAYTON

D DANDEKAR

M GREENFIELD

B LEAVY

M RAFTENBERG

RDRL WMP D

R DONEY

J RUNYEON

RDRL WMP E

P SWOBODA

RDRL WMP F

E FIORAVANTE

A FRYDMAN

N GNIAZDOWSKI

R GUPTA

R KARGUS

RDRL WMP G

41
R BANTON

N ELDREDGE

S KUKUCK

42