NASA/TM-2006-214282 

A Thermodynamically Consistent Damage 

Model for Advanced Composites

Pere Maimí 

Escola Politècnica Superior, Universitat de Girona, Girona, Spain 

Pedro P. Camanho 

Faculdade de Engenharia, Universidade do Porto, Porto, Portugal 

Joan-Andreu Mayugo 

Escola Politècnica Superior, Universitat de Girona, Girona, Spain 

Carlos G. Dávila 

Langley Research Center, Hampton, Virginia 

 

March 2006 
The NASA STI Program Office . . . in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role. 

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for NASA’s 

scientific and technical information. The NASA STI 

Program Office provides access to the NASA STI 

Database, the largest collection of aeronautical and 

space science STI in the world. The Program Office is 

also NASA’s institutional mechanism for 

disseminating the results of its research and 

development activities. These results are published by 

NASA in the NASA STI Report Series, which 

includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results ... even 

providing videos. 

For more information about the NASA STI Program 

Office, see the following: 

• Access the NASA STI Program Home Page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 

at (301) 621-0134 

• Phone the NASA STI Help Desk at 

(301) 621-0390 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7121 Standard Drive 

 Hanover, MD 21076-1320
NASA/TM-2006-214282 

A Thermodynamically Consistent Damage 

Model for Advanced Composites

Pere Maimí 

Escola Politècnica Superior, Universitat de Girona, Girona, Spain 

Pedro P. Camanho 

Faculdade de Engenharia, Universidade do Porto, Porto, Portugal 

Joan-Andreu Mayugo 

Escola Politècnica Superior, Universitat de Girona, Girona, Spain 

Carlos G. Dávila 

Langley Research Center, Hampton, Virginia 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

March 2006 
Available from: 

NASA Center for AeroSpace Information (CASI) National Technical Information Service (NTIS) 

7121 Standard Drive 5285 Port Royal Road 

Hanover, MD 21076-1320 Springfield, VA 22161-2171 

(301) 621-0390 (703) 605-6000 
A Thermodynamically Consistent Damage

Model for Advanced Composites

P. Maim ́ı a, P.P. Camanho b, J.A. Mayugo a, C.G. D ́avila c

aAMADE, Escola Polit`ecnica Superior, Universitat de Girona, Girona, Spain

bDEMEGI, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto

Frias, 4200-465 Porto, Portugal

cNASA Langley Research Center, Hampton, VA 23681,U.S.A.

Abstract

A continuum damage model for the prediction of damage onset and structural collapse of structures manufactured in fiber-reinforced plastic laminates is proposed.

The principal damage mechanisms occurring in the longitudinal and transverse directions of a ply are represented by a damage tensor that is fixed in space. Crack

closure under load reversal effects are taken into account using damage variables

established as a function of the sign of the components of the stress tensor. Damage

activation functions based on the LaRC04 failure criteria are used to predict the

different damage mechanisms occurring at the ply level. The constitutive damage

model is implemented in a finite element code. The objectivity of the numerical

model is assured by regularizing the dissipated energy at a material point using

Baˇzant’s Crack Band Model. To verify the accuracy of the approach, analyses of

coupon specimens were performed, and the numerical predictions were compared

with experimental data.

Key words: Fracture Mechanics, Continuum Damage Mechanics, Composite

Materials.

1 INTRODUCTION

The methodology for designing high-performance structures of composite materials is still evolving. The complexity of the response of composite materials

and the difficulty in predicting structural modes of failure result in the need

for a well-planned test program. The recommended practice to mitigate the

technological risks associated with such materials is to substantiate the performance and durability of the design in a sequence of steps known as the
Building Block Approach (BBA)[1]. The BBA ensures that cost and performance objectives are met by testing greater numbers of smaller less expensive

specimens, assessing technology risks early in the program, and building on the

knowledge acquired at a given level of structural complexity before progressing

to a level of more complexity.

Achieving substantiation of structural performance by testing alone can be

prohibitively expensive because of the number of specimens and components

required to characterize all loading conditions. BBA programs can achieve

significant cost reductions by seeking a synergy between testing and analysis.

The more the development relies on analysis, the less expensive it becomes.

The use of advanced analytical or numerical models for the prediction of the

mechanical behavior of composite structures can replace some of the mechanical tests and can significantly reduce the cost of designing with composites

while providing to the engineers the information necessary to achieve an optimized design.

Strength-based failure criteria are commonly used to predict failure in composite materials. A large number of continuum-based criteria have been derived to

relate stresses and experimental measures of material strength to the onset of

failure [2]-[4]. Failure criteria predict the onset of the different damage mechanisms occurring in composites and, depending on the material, the geometry

and the loading conditions, may also predict final structural collapse.

For composite structures that can accumulate damage before structural collapse, the use of failure criteria is not sufficient to predict ultimate failure.

Simplified models, such as the ply discount method, can be used to predict ultimate failure, but they cannot represent with satisfactory accuracy

the quasi-brittle failure of laminates that results from the accumulation of

several damage mechanisms.

The study of the non-linear response of quasi-brittle materials due to the

accumulation of damage is important because the rate and direction of damage propagation defines the damage tolerance of a structure and its eventual

collapse. To model the phenomena of damage propagation, non-linear constitutive models defined in the context of the mechanics of continuum mediums

have been developed and implemented in finite elements codes in recent years.

The formalism of the thermodynamics of irreversible processes is a rigorous

framework from which the constitutive models can be developed.

The simplest way to describe damage is using a single scalar damage variable

as proposed by Kachanov [5]. Damage can be interpreted as the creation of

microcavities, and the damage variables as a measure of the effective surface

density of the microdefects. Such a mechanical interpretation of damage assumes that the loads are resisted only by the undamaged ligaments in the

2
material. The stresses (σ ̃) in the ligaments, referred to as effective stresses,

continue to increase until all ligaments are severed and the material has failed.

The tensorial representation of damage is a formal and general procedure to

represent the directionality of micro-cracks, which can take any direction in a

medium depending on the load history, geometry, boundary conditions, and

material properties. After Kachanov’s pioneering work, several damage models

have been developed that describe damage as a second order tensor [6]-[9] or

as a fourth order tensor [10]-[12]. Second order tensors describe an initially

isotropic material as an orthotropic one when damage evolves, whereas fourth

order tensor models can remove all material symmetries and provide a more

general procedure to simulate damage [13].

The application of continuum damage models in orthotropic or transversely

isotropic materials, such as fiber-reinforced plastics (FRP), results in additional difficulties. The nature and morphology of a material induces some preferred directions for crack growth, i.e., crack orientations are not only induced

by the loads, geometry and boundary conditions, but also by the morphology

of the material. The interface between fiber and matrix is weaker than the

surrounding material and interfacial debonding is normally the first damage

mechanism to occur. Furthermore, residual thermal stresses occur in the composite plies due to different coefficients of thermal expansion of the fiber and

matrix (micromechanical residual thermal stresses) and due to the different

coefficients of thermal expansion in the longitudinal (fiber) and transverse

(matrix) directions (macromechanical residual thermal stresses).

Multiscale models and mesomodeling are two approaches used to evaluate

the elastic and inelastic response of a material. Using homogenization laws,

multiscale models define relations between a mesoscale, normally the scale

of the finite elements, where material is considered homogeneous, and the

microscale the scale of constituents, fiber and matrix. The constitutive models

are defined at the microscale, and the strain and stress fields at the microscale

and the mesoscale are related via transformation field tensors [14]-[18], or

solved using finite elements [18],[19]. To reduce the amount of computations

that need to be performed, periodicity of the material is invoked.

3
Mesomodeling is an alternative way to define damage models for composite

materials that is more appropriate for large scale computations. Mesomodels

treat the composite lamina [20]-[24] or sub-laminate [25] as a homogeneous

material. When diffuse damage localizes in a narrow band and becomes a

macro-crack, the response is dominated by the crack tip and its ability to dissipate energy. On the other hand, the material morphology, which is the main

basis of homogenization techniques, loses importance due to the loss of periodicity. Therefore, in structures exhibiting stable crack propagation, i.e., when

the macrocrack length does not increase under constant load, mesomodeling

is more appropriate to predict the structural collapse than multiscale models.

The main objective of the present paper is to develop a continuum damage model able to represent the quasi-brittle fracture of laminated composite

structures, from damage onset up to final structural collapse.

The majority of the material properties required by the present model can

be measured using standard test methods. Most of the material properties

that are required can be obtained from ply-based test methods. The use of

ply properties, rather than laminate properties, is an advantage because it

avoids the need to test laminates every time the lay-up or stacking sequence

is modified.

The proposed constitutive model accounts for crack closure under load reversal

effects, an important phenomenon in cases where a composite structure is

subjected to multiaxial loading.

One important issue regarding the numerical modeling of damage is that the

convergence of the solution through successive mesh refinement must be ensured. The objectivity of the numerical model is ensured by adjusting the

energy dissipated by each damage mechanism using a characteristic element

length. The constitutive model proposed herein can be integrated explicitly,

making it computationally efficient and, therefore, suitable to be used in large

scale computations.

This paper is organized as follows: a brief description of the damage mechanisms occurring in laminated composites is presented. Based on the mechanisms of damage identified, a new constitutive damage model is proposed. The

constitutive model relates the damage mechanisms with a set of internal variables. The constitutive model is implemented in a non-linear finite element

code, and it is adjusted using a procedure based on Baˇzant’s Crack Band

Model [26] that ensures the correct computation of the energy dissipated by

each damage mechanism. The accuracy of the model is assessed by comparing

the predictions with experimental data for an open-hole tension carbon-fiber

specimen.

4
2 MECHANISMS OF DAMAGE AND FRACTURE IN LAMINATED COMPOSITES

2.1 Longitudinal failure

In fiber-reinforced composites, the largest portion of the loads is resisted by the

fibers. When these fail under either tension or compression, the internal loads

must redistribute to other areas of the structure, and may cause a structural

collapse.

In composites with high fiber volume fraction and those where the strain to

failure of the resin matrix is higher than the one of the reinforcing fiber, such as

carbon-epoxy composites, longitudinal failures start by isolated fibre fractures

in weak zones. The localized fractures increase the normal and interfacial shear

stresses in the adjoining fibers, and the local stress concentrations promote

matrix cracking, fibre matrix debonding and, for ductile matrices, conical shear

failures [27]. When increasing the load further, additional fibre fractures occur,

leading to final collapse.

Failure under longitudinal tensile loading occurs in both constituents, and

fracture occurs along a plane that is parallel to the fibers and the thickness

direction (Figure 1). A simple non-interacting failure criterion based on maximum stress or maximum strain along the longitudinal direction can usually

provide an accurate measure of longitudinal tensile failure [3].

0o plies

90o plies

Fig. 1. Longitudinal failure in 0◦ plies [28].

Compressive failure of aligned fiber composites occurs from the collapse of

the fibers as a result of shear kinking and damage of the supporting matrix

[29],[30]. A kink band in a carbon-epoxy laminate resulting from compressive

longitudinal stresses is shown in Figure 2.

5
0.1mm

0o ply

Fig. 2. Kink band in a 0◦ ply [31].

Argon [32] was the first to analyze the kinking phenomenon, based on the

assumption of a local initial fiber misalignment. Fiber misalignment causes

shear stresses between fibers that rotate the fibers, which in turn increase the

shear stress and leads to instability.

Recently, the calculation of the critical kinking stress has been significantly

improved with a more complete understanding of the geometry of the kink

band as well as the incorporation of friction and material nonlinearity in the

analysis models [28],[33].

2.2 Transverse failure

Failure in the transverse direction encompasses both matrix cracking and fibermatrix debonding. Under the presence of transverse tensile stresses and inplane shear stresses, the combined effect of small defects present in a ply

such as small fiber-resin debonds, resin-rich regions, and resin voids, trigger a

transverse crack that extends through the thickness of the ply (Figure 3). The

transverse cracks are formed without disturbing the fibers: they occur at the

fiber-resin interface and in the resin.

90 ply o

Fig. 3. Transverse matrix crack in a 90◦ ply.

A fundamental issue that needs to be considered is the effect of ply thickness

on the ply strength, usually called the ’in-situ effect’. As shown in Parvizi [34]

and Chang’s [35] experiments, the constraints imposed by the neighboring

plies of different fiber orientations cause an apparent increase in the tensile

6
and shear strengths of a ply compared to those of an unconstrained ply. The

scale effect for the in-plane shear strength is represented in Figure 4 , where it

is apparent that the strength of a ply is a function of the number of 90◦ plies

stacked together using the model proposed by Camanho et al. [36].

The in-situ effect is a deterministic size effect that can be represented using

fracture mechanics models of plies containing defects [3],[37].

Fig. 4. In-situ shear strength of T300/1034-C CFRP.

Experimental results have shown that moderate values of transverse compression have a beneficial effect on the strength of a ply. This effect can be observed

in the experimental results obtained by Swanson [38],[39] shown in Figure 5.

The failure envelope calculated using the LaRC04 failure criterion [4] is also

shown in Figure 5.

When the in-plane shear stress is large compared to the transverse compressive

stress, the fracture plane is perpendicular to the midplane of the ply. This

mode of failure is referred to as Mode A in Figure 5. However, increasing

the compressive transverse stress causes a change in the angle of the fracture

plane, as shown in Figure 6. This mode of failure is referred to as Mode B

in Figure 5. Normally, for carbon-epoxy and glass-epoxy composites loaded in

pure transverse compression, the fracture plane is at an angle (fracture angle,

α0) of 53◦ ± 3◦ with respect to the thickness direction [40]. Therefore, matrix

cracking does not occur in the plane of the maximum transverse shear stress

(45◦).

7
Fig. 5. Strength as a function of transverse compression and in-plane shear.

+45o

0.1mm



90o

- o 45

Fig. 6. Matrix crack in a 45◦ ply created by in-plane compressive transverse stress

[31].

8
2.3 Delamination

Delamination is a common damage mechanism in multidirectional laminated

composites due to their weakness in the thickness direction. The three dimensional stress states that occur near geometric discontinuities such as ply

drop-offs, stiffener terminations, flanges, bonded and bolted joints, and access

holes promote delamination initiation. Delamination causes a reduction of the

bending stiffness of a composite structure and, when compressive loads are

present, promotes local buckling.

Delamination models [41]-[44] always represent a discrete crack separating elements as do, for concrete, the pioneering work of Dudgale-Barenblatt [45],[46].

The purpose of present work is to propose a Continuum Damage Mechanics

(CDM) model for the calculation of the initiation and propagation of intralaminar damage. Consequently, delamination damage is not considered here.

3 CONSTITUTIVE DAMAGE MODEL

The thermodynamics of irreversible processes is a general framework that can

be used to formulate constitutive equations. It is a logical framework for incorporating observations and experimental results and a set of rules for avoiding

incompatibilities. In this section, we present a constitutive damage model for

laminated composites that has its foundation in irreversible thermodynamics,

and that uses the LaRC04 criteria as damage activation functions.

3.1 Complementary free energy and damage operator

To establish a constitutive law, it is necessary to define a scalar function

corresponding to the complementary free energy density in the material. This

function must be positive definite, and it must be zero at the origin with

respect to the free variables (the stresses) [47]. The proposed definition for the

complementary free energy density is:

G = σ2

11

2 (1 − d1) E1

+

σ2

22

2 (1 − d2) E2

− ν12

E1

σ11σ22 +

σ2

12

2 (1 − d6) G12

+

+ (α11σ11 + α22σ22) ∆T + (β11σ11 + β22σ22) ∆M (1)

where E1, E2, ν12 and G12 are the in-plane elastic orthotropic properties of

a unidirectional lamina. The subscript 1 denotes the longitudinal (fiber) di9
rection, and 2 denotes the transverse (matrix) direction. The damage variable

d1 is associated with longitudinal (fiber) failure, whereas d2 is the damage

variable associated with transverse matrix cracking and d6 is a damage variable influenced by longitudinal and transverse cracks. α11 and α22 are the

coefficients of thermal expansion in the longitudinal and transverse directions,

respectively. β11 and β22 are the coefficients of hygroscopic expansion in the

longitudinal and transverse directions, respectively. ∆T and ∆M are the differences of temperature and moisture content with respect to the corresponding

reference values.

The stress tensor σ corresponds to the average stress tensor over a representative volume that is assumed to be much larger than the diameter of a

fiber.

To ensure the thermodynamically irreversibility of the damage process, the

rate of change of the complementary free energy G ̇ minus the externally supplied work to the solid  ̇σ : ε at constant strains, must not be negative:

G ̇ − σ ̇ : ε ≥ 0 (2)

This inequality corresponds to the positiveness of the dissipated energy and

has to be fulfilled by any constitutive model [47]. Expanding the inequality in

terms of the stress tensor and damage variables gives:



∂G

∂σ − ε



:  ̇σ +

∂G

∂d ·  ̇d ≥ 0 (3)

Since the stresses are variables that can vary freely, the expression in the

parenthesis must be equal to zero to ensure positive dissipation of mechanical

energy. Therefore, the strain tensor is equal to the derivative of the complementary free energy density with respect to the stress tensor:

ε = ∂G

∂σ = H : σ + α∆T + β∆M (4)

The lamina compliance tensor can be represented in Voigt notation as:

H = ∂2G

∂σ2 =

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎣

1

(1 − d1) E1

−ν21

E2

0

−ν12

E1

1

(1 − d2) E2

0

0 0

1

(1 − d6) G12

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎦

(5)

10
The closure of transverse cracks under load reversal, also known as the unilateral effect, is taken into account by defining four damage variables associated

with longitudinal and transverse damage. To determine the active damage

variables, it is necessary to define the longitudinal and transverse damage

modes as follows:

d1 = d1+

σ11

|σ11| + d1− −σ11|σ11|

d2 = d2+

σ22

|σ22| + d2− −σ22|σ22|

(6)

where x is the McCauley operator defined as x := (x + |x|) /2.

The present model tracks damage caused by tension loads (d+) separately

from damage caused by compression loads (d−). Depending on the sign of the

corresponding normal stress, a damage mode can be either active or passive.

The model also assumes that the shear damage variable, d6, is not affected

by the closure effect. Shear damage is caused mainly by transverse cracks and

these do not close under shear stresses (σ12). Transverse cracks are influenced

by transverse stresses (σ22) producing the closure of cracks and a friction retention whereas longitudinal cracks produce the same effect under longitudinal

stresses (σ11) [48]. The effect of friction is neglected in the present model.

3.2 Damage activation functions

The determination of the domain of elastic response under complex stress

states is an essential component of an accurate damage model. Based on the

previously described mechanisms of crack generation in advanced composites,

a strain space is considered where the material is linear elastic. In the present

model, it is assumed that the elastic domain is enclosed by four surfaces, each

of them accounting for one damage mechanism: longitudinal and transverse

fracture under tension and compression. Those surfaces are formulated by the

damage activation functions based on the LaRC03 and LaRC04 failure criteria. The LaRC03-04 failure criteria have been shown to represent accurately

the physical process of damage onset in laminated composites. The LaRC04

criteria represent an evolution of the LaRC03 criteria: some criteria such as

the one for fiber kinking, are more accurate in LaRC04. However, the improvement in accuracy is associated with a significant increase in the computational

effort. The present damage model uses a combination of both sets of criteria

to achieve a compromise between accuracy and computational efficiency. The

full details of the derivation and validation of the LaRC04 failure criteria are

presented in references [3],[4].

11
The four damage activation functions, FN , associated with damage in the longitudinal (N = 1+, 1−) and transverse (N = 2+, 2−) directions, are defined

as:

F1+ = φ1+ − r1+ ≤ 0 ; F1− = φ1− − r1− ≤ 0

F2+ = φ2+ − r2+ ≤ 0 ; F2− = φ2− − r2− ≤ 0

(7)

where the loading functions φN (N = 1+, 1−, 2+, 2−) depend on the strain

tensor and material constants (elastic and strength properties). The elastic

domain thresholds rN (N = 1+, 1−, 2+, 2−) take an initial value of 1 when

the material is undamaged, and they increase with damage. The elastic domain

thresholds are internal variables of the constitutive model, and are related to

the damage variables dM (M = 1+, 1−, 2+, 2−, 6) by the damage evolution

laws. The elastic domain threshold defines the level of elastic strains that can

be attained before the accumulation of additional damage.

3.2.1 Longitudinal tensile fracture

The LaRC04 criterion for fiber tension failure is a non-interacting maximum

allowable strain criterion defined as:

φ1+ = E1

XT

ε11 = σ ̃11 − ν12σ ̃22

XT

(8)

where the effective stress tensor σ ̃ is computed as σ ̃ = H0

−1 : ε. H0 is the

undamaged compliance tensor obtained from equation (5) using d1 = d2 =

d6 = 0.

3.2.2 Longitudinal compressive fracture

The LaRC03 failure criterion for longitudinal compressive fracture postulates

that a kink band is triggered by the onset of damage in the supporting matrix.

Under this circumstance, the fibres lose lateral support and fail under the

effect of longitudinal compressive stresses. The initial fiber misalignment and

the rotation of the fibers as a function of the applied stress state are the

parameters used in the damage activation function.

The damage activation function used to predict damage under longitudinal

compression ( ̃σ11 < 0) and in-plane shear (fiber kinking) is established as a

12
function of the components of the stress tensor σ ̃(m) in a coordinate system

(m) representing the fiber misalignment:

φ1− =

	

|σ ̃m

12| + ηLσ ̃m22


SL

(9)

where the longitudinal friction coefficient can be approximated as [3]:

ηL ≈ −

SL cos (2α0)

YC cos2 α0

(10)

The components of the effective stress tensor in the coordinate system associated with the rotation of the fibers are calculated as:

σ ̃m

22 =  ̃σ11 sin2 φC +  ̃σ22 cos2 φC − 2 |σ ̃12|sin φC cos φC

σ ̃m

12 = ( ̃σ22 − σ ̃11) sin φC cos φC + |σ ̃12|



cos2 φC − sin2 φC

 (11)

where the absolute value of the shear stress is taken because the misalignment

angle can be positive or negative.

The misalignment angle (φC) is determined using standard shear and longitudinal compression strengths, SL and XC [3]:

φC = arctan

⎛

⎜⎜⎝

1 −



1 − 4

 SL

XC + ηL

 SL

XC

2

 SL

XC + ηL



⎞

⎟⎟⎠ (12)

It should be noted that the LaRC03 failure criterion, derived to predict the

onset of damage in laminated composites, calculates the misalignment angle

as a function of the applied stress. The LaRC03 failure criterion is modified

here assuming a constant misalignment angle, corresponding to the rotation of

the fibres at failure under pure longitudinal compression. This modification of

LaRC03 failure criterion assures that φ1− is a monotonic increasing function

under any state of proportional loading.

It should be pointed out that two criteria are used in LaRC03 for fiber kinking:

Equation (9) for  ̃σm

22 ≤ 0 and a second equation for  ̃σm

22 ≥ 0. However, the

omission of the second equation results in a minor loss of accuracy because the

equation is the same as equation (13.a) (below) with the stresses transformed

into the misaligned coordinate frame.

13
3.2.3 Transverse fracture perpendicular to the laminate mid-plane (α0 = 0◦)

Transverse matrix cracks perpendicular to the mid-plane of the ply, i.e. with

α0 = 0◦ (Figures 6 and 5), are created by a combination of in-plane shear

stresses and transverse tensile stresses, or in-plane shear stresses and small

transverse compressive stresses. These conditions are represented by the following LaRC04 failure criterion:

φ2+ =

⎧

⎪⎪⎪⎨

⎪⎪⎪⎩



(1 − g)

σ ̃22

YT

+ g

σ ̃22

YT

2

+

σ ̃12

SL

2

if  ̃σ22 ≥ 0

1

SL

	

|σ ̃12| + ηLσ ̃22
if  ̃σ22 < 0

(13)

where g is the fracture toughness ratio defined as: g = GIc

GIIc .

3.2.4 Transverse fracture with α0 = 53◦

The LaRC04 matrix failure criterion for transverse compressive stresses consists of a quadratic interaction between the effective shear stresses acting on

the fracture plane:

φ2− =





τ ̃T

eff

ST

2

+



τ ̃L

eff

SL

2

if  ̃σ22 < 0 (14)

where the effective stresses  ̃τ T

eff and  ̃τ Leff are computed as [4]:

τ ̃T

eff = 	

−σ ̃22 cos (α0)



sin (α0) − ηT cos (α0) cos (θ)




τ ̃L

eff = 	

cos (α0)



|σ ̃12| + ηLσ ̃22 cos (α0) sin (θ)


 (15)

The sliding angle θ is calculated as [4]:

θ = arctan  − |σ ̃12|

σ ̃22 sin (α0)



(16)

The transverse shear strength and transverse friction coefficient can be approximated as:

ST = YC cos (α0)



sin (α0) + cos(α0)

tan(2α0)



ηT = −1

tan(2α0)

(17)

14
The fracture angle α0 is approximately 53◦ in uniaxial compression. With

increasing amounts of in-plane shear, the fracture angle diminishes up to about

40◦ (Mode B in figure 5) and then abruptly switches to 00 (Mode A in figure

5). To find the correct angle of fracture, a maximization of the LaRC03-04

failure criteria as a function of α should be performed. However, in order to

improve the computational efficiency of the present model, it is assumed that

the fracture angle can only take one of two discrete values: 0◦ or 53◦.

The elastic domain in the  ̃σ11,  ̃σ22,  ̃σ12 space represented by the LaRC04 failure

criteria is shown in Figure 7.

Fig. 7. Elastic domain in the  ̃σ11,  ̃σ22,  ̃σ12 space.

3.3 Dissipation

The rate of energy dissipation per unit volume resulting from the evolution of

damage is given by:

Ξ = ∂G

∂d1

 ̇d1 +

∂G

∂d2

 ̇d2 +

∂G

∂d6

 ̇d6 = Y1  ̇d1 + Y2  ̇d2 + Y6  ̇d6 ≥ 0 (18)

The form of the complementary free energy defined in equation (1) assures

that the thermodynamic forces (YM) conjugated to their respective damage

variables (dM) are always positive:

15
Y1 = ∂G

∂d1 = σ2

11

2(1−d1)

2E1

≥ 0

Y2 = ∂G

∂d2 = σ2

22

2(1−d2)

2E2

≥ 0

Y6 = ∂G

∂d6 = σ2

12

2(1−d6)

2G12

≥ 0

(19)

Therefore, the condition of positive evolution of damage variables (  ̇dM ≥ 0) is

a sufficient condition for the fulfillment of the second law of thermodynamics.

It is important to note that the proposed model does not generate spurious

energy dissipation, i.e., the loss or gain of mechanical energy, under crack

closure or opening. At load reversal, the time derivative of the damage variable

is non-zero (  ̇dM = 0). Considering equation (18), the thermodynamical forces,

YM, associated with the damage variables, dM, must be zero to avoid spurious

energy dissipation at load reversal [49]. This condition is trivially satisfied in

the present model (equation (19)).

Damage evolution without energy dissipation is physically inadmissible. Therefore, it is necessary to avoid damage evolution when the corresponding conjugated thermodynamic force is zero. Consider the load history represented

in Figure 8: the material is loaded in transverse tension and shear to t1 and

then loaded to t2. At time t2, the damage variable d2 evolves because the corresponding damage activation function is activated. However, the corresponding

thermodynamic force is zero (σ22 = 0, Y2 = 0).

12

~

22

~ YT

YC

SL

t1

t2

r =1 2+

r >1 2+

r =1 2-



Fig. 8. Evolution of elastic domain in the  ̃σ22 − σ ̃12 space.

This non-physical response is avoided by modifying the longitudinal and transverse damage activation functions. The transverse damage activation function

is modified using the following equation:

16
φ2− = min

⎧

⎪⎨

⎪⎩





τ ̃T

eff

ST

2

+



τ ̃L

eff

SL

2

,

σ ̃22

Ω

⎫

⎪⎬

⎪⎭

(20)

where the constant Ω is equal to  ̃σ22 when φ2−=φ2+ (Figure 8).

The longitudinal damage activation function is modified by taking into account

that under shear dominated loads, matrix cracking is the first form of damage

to occur. After matrix cracking, the transverse and shear stresses are zero and

the fiber misalignment angle φ tends to π/4. Under these circumstances, the

kink band criteria, equation (9), reads:

φ1− = min

⎧

⎨

⎩

	

|σ ̃m

12| + ηLσ ̃m22


SL

,

ηL − 1

2SL

σ ̃11

⎫

⎬

⎭

(21)

3.4 Damage evolution

The evolution of the threshold values rN is mathematically expressed by the

Kuhn-Tucker conditions:

r ̇N ≥ 0 ; FN ≤ 0;  ̇rN FN = 0 (22)

Neglecting viscous effects, the damage activation functions, equations (7), always have to be non-positive. While the damage activation function FN is

negative, the material response is elastic. When the strain state activates a

criterion, FN = 0, it is necessary to evaluate the gradient φ ̇

N . If the gradient

is not positive, the state is one of unloading or neutral loading. If the gradient

φ ̇

N is positive, there is damage evolution, and the consistency condition has

to be satisfied:

F ̇

N = φ ̇N − r ̇N = 0 (23)

Two important characteristics of the model proposed here are that the threshold values are a function of the damage variables, and that the loading functions depend on the strain tensor. Under these conditions, it is possible to

explicitly integrate the constitutive model [11],[12].

In the definition of the constitutive model, it is necessary to represent the

relation between active and inactive elastic domains. The evolution of an active

elastic domain is defined by the consistency condition, i.e., it is defined with

17
respect to the corresponding damage activation function. However, it is also

necessary to specify how the inactive elastic domain evolves if other damage

modes are active. It is assumed that the longitudinal and transverse domains

are not coupled. On the other hand, compression damage is coupled with

tension damage, as explained in the next section.

3.4.1 Transverse loading

As previously described, transverse damage in the form of matrix cracks

can have different orientations as a result of tension, shear, or compressiondominated loads. Under load reversal, transverse cracks, which are perpendicular to the ply mid-plane, do not affect the compression response: elastic

domain and the compressive damage variable (d2−) are unaffected by r2+.

On the other hand, matrix cracks at a fracture angle α0 = 53◦ caused by high

compressive transverse stresses have the same effect as cracks perpendicular to

the mid-plane (α = 0◦) when the load is reversed from compression to tension.

Therefore, the evolution of the transverse tensile elastic domain threshold (r2+)

is governed by both damage mechanisms.

Based on the above considerations, the evolution of the elastic domain in the

transverse direction can be represented by the following equations:

Tension loading:  ̇r2+ = φ ̇

2+ and  ̇r2− = 0

Compression loading:  ̇r2− = φ ̇

2− and  ̇r2+ =

⎧

⎪⎨

⎪⎩

φ ̇

2− if r2+ ≤ r2−

0 if r2+ > r2−

The integration of the previous expressions results in:

r2+ = max #

1, max

s=0,t

$

φs

2−

%

, max

s=0,t

$

φs

2+%&

r2− = max #1, max

s=0,t

$

φs

2−

%& (24)

3.4.2 Longitudinal loading

Under longitudinal tensile stresses, the fracture plane is perpendicular to the

fiber direction. When reversing the load, the cracks close and can still transfer

load. However, the broken and misaligned fibers do not carry any additional

18
load. Therefore, the compressive stiffness is influenced by longitudinal damage.

However, the elastic domain is assumed to remain unchanged.

Under longitudinal compression, damaged material consisting of broken fibers

and matrix cracks forms a kink band, and there is not a unique orientation

for the damage planes. When the loads are reversed, the cracks generated in

compression open and the elastic domain threshold increases.

Therefore, the evolution of damage thresholds for longitudinal damage are

defined as:

Tension loading:  ̇r1+ = φ ̇

1+ and  ̇r1− = 0

Compression loading:  ̇r1− = φ ̇

1− and  ̇r1+ =

⎧

⎪⎨

⎪⎩

φ ̇

1− if r1+ ≤ r1−

0 if r1+ > r1−

(25)

The integration of the previous expressions results in:

r1+ = max #

1, max

s=0,t

$

φs

1+%

, max

s=0,t

$

φs

1−

%&

r1− = max #1, max

s=0,t

$

φs

1−

%& (26)

3.5 Damage evolution laws

The internal variables rN define the threshold of the elastic domains, and are

related to the damage state of each lamina, i.e., the damage variables depend

on the values of the internal variables. In order to fully define the constitutive

model, it is necessary to define the relation between the internal variables and

the damage variables.

When material is undamaged the internal variables rN take the initial value

of 1, and dN (rN = 1) = 0. Equations (24) and (26) define the evolution of the

internal variables assuring that  ̇rN ≥ 0. As shown in equations (18) and (19),

the condition for positive dissipation is satisfied if  ̇dN ≥ 0. The condition for

positive dissipation is automatically fulfilled if the damage evolution law satisfies the condition ∂dN /∂rN ≥ 0. When the material is completely damaged, a

fracture plane is created, the strains are localized in a plane in which rN → ∞

and the related components of the stiffness tensor are zero, dN (rN → ∞) = 1.

Matrix cracks are related to the internal variables r2+ and r2−. The internal

variable r2− accounts for compressive damage only, whereas r2+ accounts for

19
both, compressive and tensile, damage. Therefore, for positive transverse normal stresses damage is a function of d2+ (r2+) because both types of cracks

(α = 0◦ and α = 53◦) are open. Under transverse compressive loads, the

damage is only influenced by the inclined cracks, d2− (r2−).

Kink bands are related to the internal variable r1−. The internal variable r1+

accounts for kink bands and fiber tensile fracture. For positive longitudinal

normal stresses, the material loses stiffness as a result of both damage modes

because the cracks open. Therefore, the damage variable can be expressed as

d1+ (r1+).

When a lamina which is fully damaged in tension (d1+ = 1) is subjected to

load reversal and the crack closes, some of the original stiffness is recovered

because the tractions can be transmitted through the crack faces. However,

the broken fibers lose their alignment. Assuming that the fibers do not carry

any load, which can be considered as a limit case, the compressive stiffness can

be approximated by the rule of mixtures applied for components in parallel

as: (1 − d1−)E1 = VmEm. The generalization of the above arguments for an

intermediate damage state can be expressed as d1− = A±

1 d1+, with:

A±

1 ≈ b VfEf

VmEm + VfEf

≈ b

E1 − E2

E1

(27)

where Ef and Em are the fiber and matrix Young modulus, Vf and Vm the

corresponding volume fractions, and b is an adjustment parameter between 0

and 1. If b = 1 the stiffness recovery is due only by the matrix, and if b = 0,

the stiffness recovery is total and it is assumed that broken fibers do not lose

alignment under compressive loads, and that the initial stiffness is recovered.

Fiber damage (d1±) is not influenced by matrix cracking (r2±) as shown in experimental results carried by Carlsson and Pipes [50] and in micromechanical

models [51]-[54] of cracked composites. Therefore, the longitudinal stiffness is

not function of matrix transverse cracks.

The shear stiffness is reduced as a result of longitudinal and transverse cracks

regardless of their orientation. Under these circumstances, the damage variable

d6 is given by:

d6 = 1 − [1 − d6(r2+)](1 − d1+) (28)

The damage evolution laws used force strain-softening as soon as one damage

activation criterion is satisfied. Softening constitutive equations result in physically inadmissable responses: the damage is localized in a plane and fracture

occurs without energy dissipation. The numerical implementation of soften20
ing constitutive equations using the finite element method results in meshdependent results because the energy dissipated is a function of the element

size.

The solution normally used to ensure the correct computation of the energy

dissipated regardless of the refinement of the mesh is to adjust the damage

evolution laws using a characteristic dimension of the finite element. The definition of the damage evolution laws is therefore related to the computational

model, and will be described in the next section.

4 COMPUTATIONAL MODEL

Two different inelastic responses have to be taken into account in the numerical simulation of damage using Continuum Damage Mechanics. While the

stress-strain response of a material exhibits a positive-definite tangent stiffness tensor, the damage zone increases along all the directions. This initial

stage of damage is commonly referred to as diffuse damage, and the numerical

solution is independent of the numerical discretization. For example, matrix

cracking in a multidirectional composite is a form of diffuse damage if the

kinematics of laminate theory is used.

When the tangent stiffness tensor is not positive definite, damage localizes

in a narrow band, and the numerical solution depends upon the numerical

discretization: decreasing the element size in the localized zone decreases the

computed energy dissipated. Therefore, the structural response is not objective

because it does not converge to a unique solution with mesh refinement.

The proposed damage model uses a constitutive model that forces localization as soon as one of the damage activation functions associated with the

onset of transverse or longitudinal cracking is satisfied, i.e., when FN = 0.

In order to guarantee that the numerical solution is independent of the discretization a characteristic element length is used in the constitutive model

using a procedure based on the crack band model proposed by Baˇzant [26].

4.1 Damage laws in softening regime: crack band model

Baˇzant’s crack band model [26] assures the objective response of the global

finite element model by regularizing the computed dissipated energy using a

characteristic dimension of the finite element and the fracture toughness:

21
gM = GM

l∗ , M = 1+, 1−, 2+, 2−, 6 (29)

where GM is the fracture toughness, gM is the energy dissipated per unit

volume, and l

∗ is the characteristic length of the finite element. For square

elements, with an aspect ratio approximately equal to one, the characteristic

element length can be approximated by the following expression [26]:

l

∗ =

√AIP

cos (γ) (30)

where |γ| ≤ 45o is the angle of the mesh lines with the crack direction and AIP

is the area associated at each integration point. For an unknown direction of

crack propagation, the average of this expression can be used,  ̄l

∗ = π

4

' π

4

0 l

∗dγ =

1.12√AIP .

When the crack propagation path can be estimated in advance, it is recommended to align the mesh with the direction of crack propagation because

cracks tend to evolve along the mesh lines. If the crack propagation is aligned

with the mesh lines, the characteristic length must be the square root of the

area corresponding to an element integration point, i.e., γ = 0.

For triangular elements, the typical characteristic length is determined by the

expression:

l

∗ = 2

AIP

√3 (31)

A more accurate measure of the characteristic element length would be obtained using the element projections for both possible crack directions, transverse and longitudinal.

The crack band model assumes that the failure process zone can be represented

by a damaged finite element zone of one element width. This approximated

method for achieving the objectivity of the global response is appropriate for

the treatment of large structures under complex damage mechanisms, such as

the ones occurring in advanced composites.

The exponential damage evolution laws proposed here are expressed in the

following general form:

dM = 1 − 1

fN (rN )

exp {AM [1 − fN (rN )]} f (rK) (32)

22
where the function fN (rN ) is selected to force the softening of the constitutive relation. f (rK) is the coupling factor between damage laws and elastic

threshold domains.

The damage evolution laws for each damage variable are:

d1+ = 1 − 1

r1+ exp[A1+(1 − r1+)]

d1− = 1 − 1

r1− exp[A1−(1 − r1−)]f(A±

1 , r1+)

d2+ = 1 − 1

f2+(r2+) exp[A2+(1 − f2+(r2+))]

d2− = 1 − 1

r2− exp[A2−(1 − r2−)]

d6 (r2+)=1 − 1

r2+ exp[A6(1 − r2+)]

(33)

where f2+ (r2+) is a function of the same order as the damage onset function

in order to force softening:

f2+ (r2+) = 1

2g



g − 1 + ((1 − g)

2 + 4gr2

2+

(34)

The coupling factor A±

1 is used for the interaction of elastic domains in the

longitudinal (fiber) direction, and it is defined from equation d1− = A±

1 d1+

and (33) for r1− = 1 as:

f



A±

1 , r1+

= 1 − A±

1 + A±1

1

r1+

exp [A1+ (1 − r1+)] (35)

where A±

1 is a material parameter defined in equation (27).

The energy dissipated per unit volume for uniaxial stress conditions is obtained

by integrating the rate of dissipation, equation (18):

gM =

) ∞

0

YM  ̇dMdt =

) ∞

1

∂G

∂dM

∂dM

∂rM

drM, M = 1+, 1−, 2+, 2−, 6 (36)

Applying the crack band model, equation (29):

) ∞

1

∂G

∂dM

∂dM

∂rM

drM = GM

l∗ , M = 1+, 1−, 2+, 2−, 6 (37)

23
Using equations (29) and (36) and substituting in (37), it is possible to calculate the damage law parameters AM that assure that the dissipated energy

computed by the numerical model is independent of mesh refinement. It is

possible to obtain analytical closed form solutions for two of the equations

(37):

A1+ = 2l

∗X2

T

2E1G1+ − l∗X2

T

(38)

A6 = 2l

∗S2

L

2G12G6 − l∗S2

L

(39)

The remaining parameters, A2± and A1−, are calculated numerically using

the algorithm presented in Appendixes A and B. The adjusting parameters

complete the definition of the constitutive model.

The material response under load reversal cycles in the transverse and longitudinal direction resulting from the proposed constitutive model are illustrated

in Figures 9 and 10 respectively.

C

D

O

A

B

E

r =1 2+ max{r , } t1

r

t2

2+ 2+ 2- r

t1

2- r =1 2- r

t2

22

22

Fig. 9. Transverse load cycle: O-A-B-O-C-D-O-E.

24
C

D

O

A

B

F

r =1 1+ max{r , } t1

r

t2

1+ 1-

r =1 1- r1-

t2

11

11

r1+

t1

r1+

t1

E

Fig. 10. Longitudinal load cycle: O-A-B-O-C-D-E-O-F.

In the transverse load cycle shown in Figure 9, it can be seen that damage

created under tensile stresses (O-A-B-O) does not influence the compressive

behavior (O-C), the size of the elastic domain (r2− = 1), or the damage

variable (d2− = 0). On the other hand, damage created under compressive

loading (C-D-O) increases the elastic domain in tension (O-E) and affects the

damage variable d2+.

In the longitudinal load cycle shown in Figure 10, it can be seen that damage

created under tensile loads (O-A-B-O) influences the compressive behavior

under load reversal (O-C), corresponding to a decrease of stiffness due the

misalignment of the fibers (d1− = A±

1 d1+). However, the elastic domain size

remains unchanged (r1− = 1).

For damage created under compressive stresses, two regions can be distinguished. For r1− < r1+ (C-D), the tensile elastic domain threshold does not

change. For r1− ≥ r1+ (D-E), both elastic domains thresholds increase to

reflect the fracture of fibers in compression.

4.1.1 Critical finite element size

The constitutive model must not lead to a local snap-back in the stress-strain

relation. In other words, the elastic energy of an element at the onset of localization, which is X2

M (l

∗)

2 t/ (2EM) with M = 1±, 2±, 6, must be lower than

or equal to the fracture energy, GMl

∗t.

25
Therefore, the maximum size for the finite element for each damage law M is:

l

∗ ≤

2EMGM

X2

M

, M = 1±, 2±, 6 (40)

where EM, GM and XM are the Young modulus, fracture energies and strengths,

respectively.

If a finite element model consist of elements larger that the maximum size

prescribed by equation (40), there is an alternative to further refining the

mesh. The snap-back in the constitutive model can be avoided by reducing

the corresponding strength [26] while taking the parameter AM to infinity:

XM =



2EMGM

l∗ (41)

Under these circumstances, the damage variables take two possible values:

dM = 0 if rN = 1 or dM = 1 if rN > 1. It is clear that the modification of the

strengths, although assuring the correct calculation of the energy dissipated,

should not be performed in the elements representing the region where crack

initiation, which is controlled by the stress tensor, takes place. In the regions

of stress concentrations, where crack initiation is likely to take place, the

mesh should be sufficiently refined to avoid any adjustment of the material

properties.

4.1.2 Fracture toughness

Each damage evolution function includes one adjusting parameter, AM, M =

1±, 2±, 6, that needs to be calculated using the corresponding component of

the fracture toughness, GM, representing the energy dissipated by inelastic

processes in the fracture process zone.

G2+ and G6 correspond to the fracture toughness of a transverse crack in

mode I and II, respectively. The mode I component of the fracture toughness, G2+, can be measured using the Double Cantilever Beam (DCB) test

(ASTM-D5528). The mode II component of the fracture toughness, G6, can

be measured using the Four-Point End Notched Flexure (4-ENF) test specimen [28].

G1+ corresponds to the mode I component of the fracture toughness for a

longitudinal crack. There is no standard test method to measure this property.

The suggested test method to measure G1+ is the Compact Tension (CT) test

specimen proposed by Pinho et al. [28]

26
The measurement of the energy dissipated that is associated with longitudinal

compressive loading is far more complex because several complex dissipative

phenomena are involved, such as crack growth, crushing and friction. These

events occur sequentially throughout the fracture process and transverse inclined cracks and longitudinal kink bands are the result of the combination of

these dissipative mechanisms.

Baˇzant et al. [55] proposed the following expression to evaluate the energy

dissipated in a kink band, G1−:

G1− = w

s

G6 (42)

where w is the kink band thickness and s is the distance between two matrix cracks. This approximation requires a good knowledge of the kink band

geometry in advance, which is a function of the external loads, the geometry

of structure and the thickness confinement [56]. This approximation does not

take into account other dissipative mechanisms in the material such as sliding of the crack faces. An alternative procedure to measure G1− is to use the

compact compression (CC) specimen that triggers kink bands in laminated

composites, as proposed by Pinho [28].

The fracture toughness for transverse compression loading, G2− can be calculated approximately using the mode II component of the fracture toughness,

the fracture angle α0 and a term accounting for friction between the crack

faces:

G2− = G6

cos α0

+ atμYC cos α0 ≈

G6

cos α0

+ atηT YC cos α0 (43)

where α0 ≈ 53 ± 3◦, t is the lamina thickness, and a is an adjustment parameter between 0 (in an unidirectional laminate) and 1 (in a strongly confined

lamina).

4.1.3 Viscous regularization

Strain-softening constitutive models cause convergence difficulties when using

global solution methods, especially for damage in the longitudinal (fiber) direction. In order to improve the convergence of the numerical algorithm, an

artificial Duvaut-Lions viscosity model [57] is implemented. The time derivatives of the internal variables associated with longitudinal failure can be

defined as:

27
r ̇1+ =

	

max $φs

1+, φs1−

%

− r1+


η

and  ̇r1− =

	

φs

1− − r1−




η

(44)

where η is the viscous parameter. When η tends to zero and max $

φs

1+, φs1−

%

−

r1+ > 0, the mathematical definition of a derivative is obtained and the

functions (44) tend to the damage thresholds evolution functions defined by

equations (25).

A numerical algorithm needs to be implemented for the time integration of

the internal variables. Using a backward-Euler scheme, the internal variables

can be updated as:

rn+1

1− = max *

rn

1−, η

η + ∆t

rn

1− +

∆t

η + ∆t

φn+1

1−

+

rn+1

1+ = max *

rn

1+, rn+1

1− , η

η + ∆t

rn

1+ +

∆t

η + ∆t

φn+1

1+ + (45)

Although some materials exhibit time-dependent response, this regularization

is implemented with the objective of improving the numerical convergence of

the model.

Two undesirable consequences can occur when increasing the viscous parameter. The first one is that the tangent relation becomes positive definite at

first stage of damage therefore localization is not ensured at damage onset.

Secondly, the energy dissipated at a material integration point undergoing

damage evolution increases with the viscous parameter.

4.2 Material tangent constitutive tensor and algorithm

The fast convergence rate of the solution algorithm for the non-linear problem

requires the correct computation of the material tangent constitutive tensor,

CT:

σ ̇ = CT :  ̇ε (46)

where:

CT = H−1 : (I − M) (47)

28
H−1 is the secant constitutive tensor, I is the identity tensor and the tensor

M is defined as:

M =

⎡

⎢

⎢

⎢

⎢

⎢

⎣

σ11

(1−d1)

2E1

∂d1

∂ε11

σ11

(1−d1)

2E1

∂d1

∂ε22

σ11

(1−d1)

2E1

∂d1

∂ε12

σ22

(1−d2)

2E2

∂d2

∂ε11

σ22

(1−d2)

2E2

∂d2

∂ε22

σ22

(1−d2)

2E2

∂d2

∂ε12

σ12

(1−d6)

2G12

∂d6

∂ε11

σ12

(1−d6)

2G12

∂d6

∂ε22

σ12

(1−d6)

2G12

∂d6

∂ε12

⎤

⎥

⎥

⎥

⎥

⎥

⎦

(48)

The scalar components of the tensor M are presented in Appendix C.

The integration of the constitutive model is performed according to the following algorithm:

1 - Read the strain tensor at time t εt

2 - Compute the effective stress tensor σ ̃t = H−1

0 : εt

3 - Compute the loading functions φt

M(σ ̃t

)

4 - Compute the threshold values rt

M(rt-1M , φtM)

5 - Compute the damage variables dt

M(rtM)

6 - Compute the nominal stress tensor σt = ,

Ht

-−1

: εt

7 - Compute the tangent constitutive tensor Ct

T = ,

Ht

-−1 : (I − Mt

)

It should be noted that the numerical calculation of the adjusting parameters

AM is only performed once, and only if the damage variable at the integration

point is greater than zero. Also, the verification of the condition AM ≥ 0

is performed in the beginning of the analysis, and, if it is not satisfied, the

material strengths are reduced according to the equation (41).

5 VALIDATION

The model developed was implemented in ABAQUS non-linear finite element

code [58] using a user-subroutine UMAT.

The model is validated by comparing the predicted failure loads of quasiisotropic laminates containing a central hole and loaded in tension with the

corresponding experimental data.

29
The experimental data used was previously obtained by Tan [59] for the

open-hole tension test specimen shown in Figure 11. The material consists

of T300/1034-C carbon fiber reinforced epoxy with a nominal ply thickness of

0.1308 mm. The elastic properties and unidirectional strengths are shown in

Tables 1 and 2 respectively.

203.2mm

25.4mm

6.35mm 

x

y

Fig. 11. Configuration of the open-hole tension test specimen.

Table 1

T300/1034-C elastic properties [35],[59].

E1 (GPa) E2 (GPa) G12 (GPa) ν12

146.8 11.4 6.1 0.30

Table 2

T300/1034-C unidirectional strengths (MPa) [59],[35].

XT XC YT YC SL

1730.0 1379.0 66.5 268.2 58.7

The components of the fracture toughness for the different damage models

are required for both the determination of the in-situ strengths and of the

adjusting parameters A±

M.

30
The components of the fracture toughness associated with matrix cracking

used in the model were measured by Shahid and Chang [60] for T300/934

carbon-epoxy laminates. The components of the fracture toughness associated with longitudinal tensile and compressive fracture used were obtained by

Pinho [28] for a carbon-epoxy composite using the same fiber type (T300/913).

The values of fracture toughness used in the model are shown in Table 3.

Table 3

Fracture toughness (N/mm.).

G2+ G6 G2− G1+ G1−

0.23 0.46 0.76 89.83 78.27

The damage activation functions are established in terms of the in-situ strengths.

The in-situ strengths are calculated from the closed-form equations proposed

in [36], using a shear response factor β = 3.2 × 10−8 mm6/N3

. The resulting

in-situ strengths are shown in Table 4.

Table 4

In-situ strengths (MPa).

Embedded ply Outer ply Thick ply

YT 158.8 101.2 105.4

SL 109.5 89.8 73.4

The coefficients of thermal expansion of the material used are α11 = −1.0 ×

10−6/◦C and α22 = 26.0 × 10−6/◦C, and the temperature change is ∆T =

−152◦C [61]. The effects of moisture absorption are neglected in the model.

The lay-ups tested in specimens with the geometry shown in Figure 11 are

[0◦/[±45◦]3/90◦

3]s,[0◦/[±45◦]2/90◦5]s and [0◦/ ± 45◦/90◦7]s.

The finite element model created use 4-node shell elements. Based on the properties reported in Tables 1-4, the maximum element size that avoids lowering

the strength is 0.508 mm. The mesh in the vicinity of the hole, corresponding

to the region where damage occurs, consist of elements with sizes ranging from

0.127 mm. (hole edge) to 0.635 mm. (specimen edge).

It is worth noting that shell elements based on lamination theory assume

a linear strain field along laminate thickness expressed by two tensors: mid

plane strains (ε0) and curvatures (κ). This simplified kinematic description

may not be able to detect the localization of strains in a ply constrained by

sub-laminates. In other words, damage mechanisms, such as matrix transverse cracking in just one ply of a multidirectional laminate, correspond to

31
distributed or diffuse damage if lamination theory is employed. Localization

of damage occurs when the different damage mechanisms occur in all plies of

the laminate.

Before the localization of damage in the laminate, the ply constitutive models

that can be used are based on strain softening, such as in the model proposed

here, or on elastic analysis of cracked plies. The latter solution is normally only

obtainable for a periodic distribution of transverse matrix cracks in central 90◦

plies of rectangular laminates under constant stresses [53],[54].

Figure 12 shows the load-displacement relation of the three specimens simulated.

Fig. 12. Predicted load-displacement relations.

Table 5 presents the predicted and experimental failure stresses, σN , defined

using the failure load Pu, and the specimen width and thickness, w and t

respectively, as σN = P u

wt = 9.71Pu.

Table 5

Predicted and measured failure stress, σN (MPa).

Lay-up Experimental [59] Predicted Error (%)

[0◦/[±45◦]3/90◦

3]s 235.8 225.5 -4.4

[0◦/[±45◦]2/90◦

5]s 185.8 192.4 3.7

[0◦/ ± 45◦/90◦

7]s 160.0 139.3 -12.9

The predicted failure loads are in excellent agreement with the experimental

failure loads measured by Tan [59].

32
6 CONCLUSIONS

A new constitutive model for the prediction of damage onset, growth and

ultimate failure of composite structure under plane stress was proposed. The

onset of the different intralaminar damage mechanisms is predicted using a

simplification of the LaRC04 failure criteria.

The constitutive model proposed is based on four ply fracture planes and

accounts for the unilaterality of damage by its ability to represent complex

load histories, including tension-compression load reversals.

The constitutive law was implemented in a computational model that ensures that the computed dissipated energy is independent of the discretization.

Therefore, the numerical solution is objective with respect to mesh refinement.

The computational model developed was used in the simulation of open-hole

test specimens loaded in tension using different lay-ups. An excellent agreement between the predicted and measured failure stresses was obtained.

Acknowledgements

The research stays of the first author at the University of Porto, Portugal, were

funded by the University of Girona (BE-UdG-2004) and by the Spanish government under ”Acciones Integradas Hispano-Portuguesas” (HIP2004-0031).

The first and third authors also acknowledges the funds of the Spanish government through DGICYT under the contract MAT 2003-09768-C03-01.

The financial support of the Portuguese Foundation for Science and Technology (FCT) under the project PDCTE/50354/EME/2003 is acknowledged by

the second author.

33
References

[1] MIL-HDBK-17-3F, Military Handbook, Polymer Matrix Composites. U.S.

Department of Defense (2002).

[2] Soden P. D., Hinton M. J., Kaddour A. S., A comparison of the predictive

capabilities of current failure theories for composite laminates, Composites

Science and Technology, Vol. 58, No. 7, (1998) pp. 1225-1254.

[3] D ́avila C. G., Camanho P. P., Rose C. A., Failure criteria for FRP laminates,

Journal of Composite Materials, 39 (2005) pp. 323-345.

[4] Pinho S. T., D ́avila C. G., Camanho P. P., Ianucci L., Robinson P., Failure

models and criteria for FRP under in-plane or three-dimensional stress states

including shear non-linearity. NASA/TM-2003-213530 (2004).

[5] Kachanov L. M., Time of the rupture process under creep conditions. Isv. Nauk

SSSR, Otd. Tehk. Nauk, 8 (1958) pp. 26-31.

[6] Chaboche J-L., A continuum damage theory with anisotropic and unilateral

damage. La Recherche A ́erospatiale 2 (1995) pp. 139-147.

[7] Carol I., Rizzi E., Willam K., On the formulation of anisotropic elastic

degradation. I. Theory based on a pseudo-logarithmic damage tensor rate. Int.

J. Solid Structures 38 (2001) pp. 491-518.

[8] Carol I., Rizzi E., Willam K., On the formulation of anisotropic elastic

degradation. II. Generalized pseudo-Rankine model for tensile damage. Int. J.

Solid Structures 38 (2001) pp. 519-546.

[9] Luccioni B., Oller S., A directional damage model. Comput. Methods Appl.

Mech. Engrg. 192 (2003) pp. 1119-1145.

[10] Ortiz M., A constitutive theory for inelastic behaviour of concrete. Mechanics

of Materials 4 (1985) pp. 67-93.

[11] Simo J. C., Ju J. W., Strain and stress-based continuum damage models-I.

Formulation. Int. J. Solids Structures 23(23) (1987) pp. 821-840.

[12] Simo J. C., Ju J. W., Strain and stress-based continuum damage models-II.

Computational aspects. Int. J. Solids Structures 23(23) (1987) pp. 841-869.

[13] Cauvin A., Testa R. B., Damage mechanics: basic variables in continuum

theories. Int. J. Solid Structures 36 (1999) pp. 747-761.

[14] Chaboche J-L., Kruch S., Maire J. F., Pottier T., Towards a micromechanics

based inelastic and damage modeling of composites. Int. J. of plasticity 17 (2001)

pp. 411-439.

[15] Fish J., Yu Q., Two-scale damage modeling of brittle composites. Comp. Science

and Technology 61 (2001) pp. 2215-2222.

34
[16] Fish J., Yu Q., Shek K., Computational damage mechanics for composite

materials based on mathematical homogenization. Int. J. Numer. Meth. Engng.

45 (1999) pp. 1657-1679.

[17] Voyiadjis G. Z., Deliktas B., A coupled anisotropic damage model for the

inelastic response of composite materials. Comput. Methods Appl. Mech. Engrg.

183 (2000) pp. 159-199.

[18] Car E., Zalamea F., Oller S., Miquel J., O ̃nate E., Numerical simulation of

fiber reinforced composite materials — two procedures. Int. J. Solid Structures

39 (2002) pp. 1967-1986.

[19] Oller S., Miquel Canet J., Zalamea F., Composite material behavior using a

homogenization double scale method. Journal of Engineering Mechanics 31(1)

(2005) pp. 65-79.

[20] Ladev`eze P., Allix O., De ̈u J-F., L ́evˆeque D., A mesomodel for localisation and

damage computation in laminates. Comput. Methods Appl. Mech. Engrg. 183

(2000) pp. 105-122.

[21] Barbero E. J., Lonetti P., A damage model for composites defined in terms of

available data. J. Mech. Compos. Mater. Struct. 8(4) (2001) pp. 299-316.

[22] Barbero E. J., Devivo L., A constitutive model for elastic damage in fiberreinforced PMC laminae. J. Damage. Mech. 10(1) (2001) pp. 73-93.

[23] Matzenmiller A., Lubliner J., Taylor R. L., A constitutive model for anisotropic

damage in fiber-composites. Mechanics of Materials 20 (1995) pp. 125-152.

[24] Oller S., Botello S., Miquel J., O ̃nate E., An anisotropic elastoplastic model

based on an isotropic formulation. Engineering Computations 12 (1995) pp.

245-262.

[25] William K. V., Vaziri R., Poursartip A., A physically based continuum damage

mechanics model for thin laminated composite structures. Int. J. Solids and

Structures 40 (2003) pp. 2267-2300.

[26] Baˇzant Z. P., Oh B. H., Crack band theory for fracture of concrete. Materials

and structures 16 (1983) pp. 155-177.

[27] Daniel I. M., Ishai O., Engineering mechanics of composite materials. Oxford

Univerity Press (1994).

[28] Pinho S. T., Modelling failure of laminated composites using physically-based

failure models. PhD Thesis, Department of Aeronautics, Imperial College

London, U.K. (2005).

[29] Fleck N. A., Liu D., Microbuckle initiation from a patch of large amplitude fibre

waviness in a composite under compression and bending. European Journal of

Mechanics - A/Solids, 20(1) (2001) pp. 23-37.

[30] Schultheisz C. R., Waas A. M., Compressive failure of composites, Part 1:

Testing and micromechanical theories. Progress in Aerospace Sciences, 32

(1996) pp. 1-42.

35
[31] Camanho, P. P., Application of numerical methods to the strength prediction

of mechanically fastened joints in composite laminates. PhD Thesis, Centre for

Composite Materials, Department of Aeronautics, Imperial College London,

U.K. (1999).

[32] Argon A. S., Fracture of Composites, Vol. 1. Academic Press, New York, (1972).

[33] Sutcliffe M. P. F., Fleck N. A., Microbuckle propagation in fibre composites.

Acta mater. 45(3) (1997) pp. 921-932.

[34] Parvizi A., Garrett K., Bailey J., Constrained cracking in glass fibre-reinforced

epoxy cross-ply laminates. Journal of Material Science, 13 (1978) pp. 195-201.

[35] Chang F. K., Chen M. H., The in-situ ply shear strength distributions in

graphite/epoxy laminated composites. Journal of Composite Materials. 21

(1987) pp. 708-733.

[36] Camanho P. P., D ́avila C. G., Pinho S. T., Iannucci, L., Robinson, P., Prediction

of in-situ strengths and matrix cracking in composites under transverse tension

and in-plane shear. Composites- Part A, in press (2005).

[37] Dvorak G. J., Laws N., Analysis of progressive matrix cracking in composite

laminates II. First ply failure. Journal of Composite Materials. 21 (1987) pp.

309-329.

[38] Swanson S. R., Messick M. J., Tian Z., Failure of carbon epoxy lamina under

combined stresses. J. Composite Materials, 21 (1987) pp. 619-630.

[39] Swanson S. R., A micro-mechanical model for in-situ compression strength of

fiber composite laminates. Transactions of the American Society of Mechanical

Engineers. Series H, Journal of Engineering Materials and Technology 114

(1992) pp. 8-12.

[40] Puck A., Sch ̈urmann H., Failure analysis of FRP laminates by means of

physically based phenomenological models. Composites Science and Technology.

58 (1998) pp. 1045-1067.

[41] Camanho P. P., D ́avila C. G., Mixed-mode decohesion finite elements for the

simulation of delamination in composite materials. NASA/TM-2002-211737

(2002).

[42] Turon A., Camanho P. P., Costa J., D ́avila C. G., An interface damage model

for the simulation of delamination under variable mode ratio in composite

materials. NASA TM 213277 (2004).

[43] Camanho P. P., D ́avila C. G., Ambur D. R., Numerical simulation of

delamination growth in composite materials. NASA-TP-2001-211041 (2001).

[44] Allix O., Daudeville L., Ladev`eze, P., Delamination and damage mechanics.

Mechanical Engineering Publications, London. (1991) pp. 143-157.

[45] Dudgale D. S., Yielding of steel sheets containing slits. Journal of Mechanics

and Physics of Solids 8 (1960) pp. 100-104.

36
[46] Barenblatt G. I., Mathematical theory of equilibrium cracks in brittle fracture.

Advances in Applied Mechanics 7 (1962).

[47] Malvern L. E., Introduction to the mechanics of a continuous medium. PrenticeHall, Englewood Cliffs (1969).

[48] Chaboche J-L., Maire J. F., A new micromechanics based CDM model and its

application to CMC’s. Aerospace Science and Technology 6 (2002) pp. 131-145.

[49] Carol I., Willam K. Spurious dissipation/generation in stiffness recovery models

for elastic degradation and damage. Int. J. Solid Structures 33 (1996) pp. 2939-

2957.

[50] Carlsson L. A., Pipes R. B., Experimental characterization of composite

materials. Prentice-Hall (1987).

[51] Dvorak G. J., Laws N., Hejazi M., Analysis of progressive matrix cracking in

composite laminates I. thermoelastic properties of a ply with cracks. Journal of

Composite Materials. 19 (1985) pp. 216-234.

[52] Laws N., Dvorak G. J., Hejazi M., Stiffness changes in unidirectional composites

caused by crack systems. Mechanics of Materials. 2 (1983).

[53] Nuismer R. J., Tan S. C., Constitutive relations of a cracked composite lamina.

Journal of Composite Materials. 22 (1988) pp. 306-321.

[54] Tan S. C., Nuismer R. J., A theory for progressive matrix cracking in composite

laminates. Journal of Composite Materials. 23(1989) pp. 1029-1047.

[55] Baˇzant Z. P., Kim J-J. H., Daniel I. M., Becq-Giraudon E., Zi G., Size effect

on compression strength of fiber composites failing by kink band propagation.

International Journal of fracture. 95 (1999) pp. 103-141.

[56] McGee J. D., Nemat-Nasser S., Dynamic bi-axial testing of woven composites.

Material Science and Engineering A317 (2001) pp. 135-139.

[57] Duvaut G., Lions J. L., Les in ́equations en m ́ecanique et en physique. Dunod,

Paris (1972).

[58] ABAQUS 6.5 User’s Manual, ABAQUS Inc., Pawtucket, RI, U.S.A. (2005).

[59] Tan S. C., A progressive failure model for composite laminates containing

openings. Journal of Composite Materials. 25 (1991) pp. 556-577.

[60] Shahid I. S., Chang F. K., An accumulative damage model for tensile and

shear failures of laminated composite plates. Journal of Composite Materials.

29 (1995) pp. 926-981.

[61] Springer G. S., Environmental effects on composite materials. Technomic

Publishing Company, Inc., Westport, U.S.A. (1981).

37
Appendix A: Integration Algorithm

In order to calculate the adjusting parameters AM used in the damage evolution laws, it is necessary to integrate the following equation numerically:

gM =

) ∞

1

 1 − ν12ν21

1 − (1 − dM) ν12ν21 2

σ ̃2

M

2EM

∂dM

∂rN

drN (A-1)

The Simpson method of numerical integration approximates the solution using

quadratic polynomials. The general form of the polynomial can be expressed

as:

gM 


h

3



f 0

M + ... + 4f odd

M + 2f evenM + ... + fnM



(A-2)

where h is the step increment, and fi

M =

 1−ν12ν21

1−(1−di

M)ν12ν21 2 (σ ̃i

M)

2

2Ei

M

∂di

M

∂ri

N

is defined

between r = 1 and r → ∞.

Since the damage laws selected tend to zero, it is necessary to define a point to

stop the integration. When the stress becomes less than K times the stress at

the onset of localization, the remaining energy can be neglected. The increment

(h) can be selected defining the number of steps (n) as:

hi ≈ −

1

nA ln  1K



(A-3)

The algorithm is implemented as:

38
1- Select the parameters n and K

2- Initialize r=1 , G=0 and CONT=0

3- Compute step size h

4- WHILE CONT<n

DO I=1:3

f (I) =  1−ν12ν21

1−(1−dM)ν12ν21 2 σ ̃2

M

2EM

∂dM

∂rN

r=r+h

END DO

r=r-h

g = g + h

3 (f (1) + 4f (2) + f (3))

CONT=CONT+1

END WHILE

39
Appendix B: Secant method to determine the parameters AM

To find the adjustment parameters for the damage law, it is necessary to

integrate the stress-strain relation in terms of the unknown parameter AM. The

integration can be done numerically with the algorithm presented in Appendix

A. To iterate on the value of AM, the secant method is used. The problem to

be solved can be expressed as:

gM (AM) − GM

l∗ = 0 (B-1)

To select the two parameters to start the iteration, the following approximation is used:

A1

M = 2l

∗X2

M

2EMGM − l∗X2

M

and A0

M = 0.5A1M (B-2)

The function gi (Ai) is only defined for positives values of A. Defining the

minimization function as:

ln 

Aj+1

M



= ln Aj

M



−

.

ln gj

M



− ln GM

l∗

/ ln 

Aj

M



− ln Aj−1

M



ln gj

M



− ln gj−1

M

 (B-3)

the following algorithm is proposed:

1- Initialize A1

M = 2l

∗X2

M

2EMGM−l∗X2

M

, A0

M = 0.5A1M and j=1

2- Integrate numerically g0

M , see Appendix A

WHILE

0

0

0gM (AM) − GM

l∗

0

0

0 ≤ tol

Integrate numerically gj

M , see Appendix A

Aj+1

M = exp

ln(Aj

M)−



ln(gj

M)−ln GM

l∗

 ln(Aj

M)−ln(Aj−1M ) ln(gj

M)−ln(gj−1

M )

gj+1

M = gjM

Aj

M = Aj−1

M

j=j+1

END WHILE

40
Appendix C: Material tangent stiffness tensor

If artificial viscosity is considered, incremental damage laws have to be implemented:

∂d1

∂r1

= ∂dn+1

1

∂rn+1

1

∂rn+1

1

∂rn

1

= ∂dn+1

1

∂rn+1

1

 η

∆t + η

∂rn

1

∂rn−1

1

+

∆t

∆t + η



(C-1)

Applying the chain rule, the damage evolution can be written as:

.

∂d1

∂ε11

∂d1

∂ε22

∂d1

∂ε12 /

= ∂d1

∂r1

.

∂r1

∂ε11

∂r1

∂ε22

∂r1

∂ε12 /

.

∂d2

∂ε11

∂d2

∂ε22

∂d2

∂ε12 /

= ∂d2

∂r2

.

∂r2

∂ε11

∂r2

∂ε22

∂r2

∂ε12 /

.

∂d6

∂ε11

∂d6

∂ε22

∂d6

∂ε12 /

= (1 − d6 (r1)) ∂d6(r2)

∂r2

.

∂r2

∂ε11

∂r2

∂ε22

∂r2

∂ε12 /

+

+ (1 − d6 (r2)) ∂d6(r1)

∂r1

.

∂r1

∂ε11

∂r1

∂ε22

∂r1

∂ε12 /

Derivation of softening damage laws:

∂dN

∂rM

= 1+rMAN

r2

M

exp [AN (1 − rM)]

∂d1−

∂r1−

=

⎧

⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

1+r1−A1−

r2

1−

exp [A1− (1 − r1−)] f



A±

1 , r1+

if r1+ > r1−

1+r1−A1−

r2

1−

exp [A1− (1 − r1−)] K1− if r1+ = r1−

∂d2+

∂r2+

= 4 gr2+ exp(A2+(1−K2+))(g(2+A2+)−A2+(1−K2+))

K2+(g−1+K2+)

∂d6 (r1)

∂r1

= 1+r1+A1+

r2

1+

exp (A1+ (1 − r1+))

Where:

K1− = f



A±

1 , r1+

(1 + r1−A1−) + A±

1

r1− (1 + r1−A1+) exp [A1+ (1 − r1−)]

f



A±

1 , r1+

= 1 − A±

1 + A±1

1

r1+

exp [A1+ (1 − r1+)]

and K2+ =

(

(1 − g)

2 + 4gr2

2+

41
Evolution of damage thresholds respect the strains:

.

∂r1+

∂ε11

∂r1+

∂ε22

∂r1+

∂ε12 /

= E1

XT

.

100 /

.

∂r1−

∂ε11

∂r1−

∂ε22

∂r1−

∂ε12 /

=

⎧

⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

1

SL

.

0 ηL sign ( ̃σm

12)

/

[R] [C0] if |σ ̃m

12|+ηLσ ̃m22

SL < −E1YC ε11

E1

YC

.

100 /

if |σ ̃m

12|+ηLσ ̃m22

SL > −E1YC ε11

.

∂r2+

∂ε11

∂r2+

∂ε22

∂r2+

∂ε12 /

=

⎧

⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

1

r2+ .

0 1−g

2YT + gσ ̃22

Y 2

T

σ ̃12

S2

L

/

[C0] if  ̃σ22 > 0

1

SL

.

0 ηL sign ( ̃σ12)

/

[C0] if  ̃σ22 < 0

.

∂r2−

∂ε11

∂r2−

∂ε22

∂r2−

∂ε12 /

=

⎧

⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎩

1

r2−

σ ̃T

eff

S2

T

∂σ ̃T

eff

∂ε + σ ̃L

eff

S2

L

∂σ ̃L

eff

∂ε 

if σ ̃T

eff

ST

2

+

σ ̃L

eff

SL

2

< ε22

E2

YT

.

010 /

if σ ̃T

eff

ST

2

+

σ ̃L

eff

SL

2

> E22

Y T ε22

Where:

[R] =

⎡

⎢

⎢

⎢

⎢

⎢

⎣

cos2 φC sin2 φC 2 cos φC sin φCsign ( ̃σ12)

sin2 φC cos2 φC −2 cos φC sin φCsign ( ̃σ12)

− cos φC sin φC cos φC sin φC 

cos2 φC − sin2 φC



sign ( ̃σ12)

⎤

⎥

⎥

⎥

⎥

⎥

⎦

∂σ ̃T

eff

∂ε =

.

0 − cos α0



sin α0 − ηT cos α0

K



+ ηT

K3

 σ ̃12

tan α0 σ ̃22 2

− ηT σ ̃12

σ ̃22K3 tan2 α0

/

[C0]

∂σ ̃L

eff

∂ε =

.

0 ηL cos2 α0

K3 sin3 α0

|σ ̃12|

3

σ ̃3

22

ηL cos α0

K tan α0 + sign ( ̃σ12) cos α0



1 − ηL cos α0

K3 sin3 α0

σ ̃2

12

σ ̃2

22

 /

[C0]

and K =



1 + σ ̃2

12

σ ̃2

22 sin2 α0

42
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The 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 of information. Send comments regarding this burden estimate or 

any other aspect of this collection of information, including suggestions for reducing this 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) 

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER 

5b. GRANT NUMBER 

5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 5d. PROJECT NUMBER 

5e. TASK NUMBER 

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. SPONSORING/MONITOR'S ACRONYM(S) 

11. SPONSORING/MONITORING

REPORT NUMBER 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

15. SUBJECT TERMS 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 

ABSTRACT 

18. NUMBER 

OF 

PAGES 

19a. NAME OF RESPONSIBLE PERSON 

a. REPORT b. ABSTRACT c. THIS PAGE 

19b. TELEPHONE NUMBER (Include area code) 

Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std. Z39-18