Bertram Stier

RWTH Aachen University, D-52074 Aachen, Germany

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Jaan W. Simon and Stefanie Reese

RWTH Aachen University, D-52074 Aachen, Germany

Investigation of Micro-Scale Architectural

Effects on Damage of Composites

NASA/TM—2015-218740

May 2015
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated 

to the advancement of aeronautics and space science. 

The NASA Scientifi c and Technical Information (STI) 

Program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI Program provides access 

to the NASA Technical Report Server—Registered 

(NTRS Reg) and NASA Technical Report Server—

Public (NTRS) thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counter-part of peer-reviewed formal 

professional papers, but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings 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. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 757-864-6500

• Telephone the NASA STI Information Desk at

 757-864-9658

• Write to:

NASA STI Program

 Mail Stop 148

 NASA Langley Research Center

 Hampton, VA 23681-2199
Bertram Stier

RWTH Aachen University, D-52074 Aachen, Germany

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Jaan W. Simon and Stefanie Reese

RWTH Aachen University, D-52074 Aachen, Germany

Investigation of Micro-Scale Architectural

Effects on Damage of Composites

NASA/TM—2015-218740

May 2015

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Available from

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by technical management. 

NASA STI Program

Mail Stop 148

NASA Langley Research Center

Hampton, VA 23681-2199

National Technical Information Service

5285 Port Royal Road

Springfi eld, VA 22161

703-605-6000

This report is available in electronic form at http://www.sti.nasa.gov/ and http://ntrs.nasa.gov/
Abstract

This paper presents a three-dimensional, energy based, anisotropic, stiffness reduction, progressive damage model

for composite materials and composite material constituents. The model has been implemented as a user-defined

constitutive model within the Abaqus finite element software package and applied to simulate the nonlinear behavior

of a damaging epoxy matrix within a unidirectional composite material. Three different composite microstructures

were considered as finite element repeating unit cells, with appropriate periodicity conditions applied at the boundaries. Results representing predicted transverse tensile, longitudinal shear, and transverse shear stress-strain curves

are presented, along with plots of the local fields indicating the damage progression within the microstructure. It is

demonstrated that the damage model functions appropriately at the matrix scale, enabling localization of the damage

to simulate failure of the composite material. The influence of the repeating unit cell geometry and the effect of the

directionality of the applied loading are investigated and discussed.

1. Introduction

Progressive damage modeling of polymer matrix composites has received a great deal of attention in recent years

as predictive capabilities for the complex nonlinear behavior of these materials are sought. A significant outstanding

question is in regards to the proper, or most beneficial, scale on which the material damage model should function.

Modeling the composite at the scale of the constituents (i.e., fiber and matrix materials) enables use of models for

the matrix material that are initially isotropic and avoids ad-hoc coupling rules for the various damage components.

However, a validated micromechanics method is required, and often, the constituent scale data best suited for characterizing the damage model are unavailable. In contrast, the macro-scale approach, wherein the composite is treated as

an effective anisotropic material, will always be more computationally efficient, yet it requires a more complex damage

model to handle the extreme anisotropy of the composite material with its directionally dependent multiple damage

mechanisms. Further, the way in which the damage components interact (i.e., how much shear stiffness reduction

results from a given amount of normal stiffness reduction) must be characterized, or, more commonly, assumed.

Herein, an anisotropic, three-dimensional, progressive damage model, capable of modeling composites at the

macro (composite) scale, as well as at the micro (constituent) scale, is presented. The development follows the

approach outlined by Barbero [1] in the case of plane stress. The damage model uses the continuum damage approach

first introduced by Kachanov [2, 3] in which continuous damage variables are introduced to quantify the stiffness

reduction of a material in various directions. The theoretical background of this approach, and specifically, the elastic

stiffness reduction damage approach, is discussed in detail by Lemaitre [5] and Chaboche [6], Kruch and Chaboche

[18], Kajcinovic [7], Voyiadjis and Kattan [8], and Talreja [9, 10]. Further, Ladeveze [14, 15] and Pineda et al. [12]  ́

have reported progressive stiffness reduction damage models aimed at polymer matrix composites at the composite

level. These models are motivated by the microstructure, but cannot consider micro scale effects like random fiber

distribution and void content/morphology. Pineda et al. [11] and Bednarcyk et al. [13] presented stiffness reduction

NASA/TM—2015-218740 1

Investigation of Micro-Scale Architectural 

Effects on Damage of Composites 

Bertram Stier 

Institute of Applied Mechanics 

RWTH Aachen University 

Mies-van-der-Rohe-Strae 1 

D-52074 Aachen, Germany 

Brett A. Bednarcyk 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Jaan W. Simon and Stefanie Reese 

Institute of Applied Mechanics 

RWTH Aachen University 

Mies-van-der-Rohe-Strae 1 

D-52074 Aachen, Germany 
damage models operating at the micro (fiber/marix) scale, but these studies did not examine micro scale architectural

effects. Micromechanical investigations considering random fiber distributions of unidirectional (UD) fiber reinforced

plastics (FRP) using an elasto-plastic isotropic damage model for the matrix material with the focus on fiber-matrix

debonding were carried out by Melro et. al. [16, 17]. Chen and Ghosh [19] studied the strain rate dependent

deformation and failure of composite materials at the micro-scale for different ideally packed unit cells under different

strain rates and loading conditions. However, none of these works studied the influence of the coupling of the damage

variables.

The progressive damage model presented below is characterized for a polymer resin based on neat matrix tensile data. The model has been implemented within Abaqus [4] and applied to model the matrix material within a

unidirectional polymer matrix composite. Three different composite repeating unit cells (with constant fiber volume

fraction), representing three different fiber architectures, have been considered: hexagonal dense packing (i.e., equilateral hexagonal packing), and 5 and 15 fiber repeating unit cells (RUCs) with random fiber packing. The boundary

conditions necessary to enforce periodicity at the RUC boundaries, while preserving the ability to apply arbitrary

global (far-field) strains, are also presented. Results are shown for the three RUC microstructures for applied transverse tension, longitudinal shear, and transverse shear as a function of a damage model parameter that controls the

extent of shear damage. The functionality of the damage model and its ability to predict damage localization, are

highlighted, as are the effects of the directionality of the applied transverse tension and longitudinal shear loading,

with respect to the RUC geometries.

2. Progressive Damage Model

A stiffness reduction progressive damage model has been developed and implemented based on the two-dimensional

approach given by Barbero [1]. For application to the matrix material within a composite, a three-dimensional extension of the model is needed due to the local three-dimensional stress state within a composite material. For a pure

monolithic matrix, it is reasonable to consider the material to be initially isotropic, however when subjected to deformation, the formulation enables the material to exhibit anisotropic behavior depending on the activation (influence) of

the damage mechanisms in the different directions. To allow for arbitrary damage, the in situ material model should

have the capability to allow for different coupling of the damage variables, to represent, for example, penny shaped

cracks or cavity type damage. The consequence is a very general formulation of the damage model, wherein the physical interpretation of the material parameters is maintained. The model is based on the principle of energy equivalence

[2] for small strains, and utilizes six scalar damage variables associated with the Young’s and shear moduli. Thus, the

proposed Gibbs free energy is given by,

χ = 1

2ρ

[ σ2

11

(1 − D1)2E ̃1

+ σ2

22

(1 − D2)2E ̃2

+ σ2

33

(1 − D3)2E ̃3

− (

ν ̃21

E ̃2

+ ν ̃12

E ̃1

) σ11σ22

(1 − D1)(1 − D2)

− (

ν ̃31

E ̃3

+ ν ̃13

E ̃1

) σ11σ33

(1 − D1)(1 − D3)

− (

ν ̃32

E ̃3

+ ν ̃23

E ̃2

) σ22σ33

(1 − D2)(1 − D3)

+ σ2

12

(1 − D12)2G ̃

12

+ σ2

13

(1 − D13)2G ̃

13

+ σ2

23

(1 − D23)2G ̃

23

]

(1)

where ρ is the density of the material, σi j are the stress components, D = [D11D22D33D12D13D23] are the damage variables, and E ̃i, G ̃i j and  ̃νi j are the pristine (undamaged) Young’s moduli, shear moduli, and Poisson’s ratios,

respectively. The material compliance tensor, S, is then given by the second derivative of the Gibbs free energy with

respect to the stress tensor,

S = ρ

∂2χ

∂σ2 (2)

NASA/TM—2015-218740 2
The formulation relies on a damage surface, g(Y(D), γ(δ)), and a damage potential, f(Y(D), γ(δ)), which are

functions of thermodynamic forces, Y, and a hardening/softening function, γ(δ), with δ being the damage hardening/softening variable. Moreover, it is associative, which requires f = g. When damage is accumulating,

f = g = gˆ − γˆ = 0 (3)

with

gˆ = √

YHY (4)

γˆ = γ + γ0 = c1(e

δ

c2 − 1) + γ0 (5)

where c1 and c2 are material constants, γ0 is the damage threshold, and H is a matrix of material parameters

responsible for the interaction of the damage variables’ effect on the stiffness reduction of the material. If f is less

than zero, damage does not evolve.

The vector of thermodynamic driving forces is given by

Y = ρ

∂χ

∂D (6)

and the damage evolution is

D ̇ = λ ̇ ∂ f

∂Y (7)

where the Lagrangian (damage) multiplier λ ̇ can be interpreted as the damage rate magnitude and ∂ f

∂Y is the direction of damage evolution in Y-space. Further, the damage rate magnitude affects the hardening variable rate as

δ ̇ = λ ̇ ∂g

∂γ (8)

Making use of the consistency condition

g ̇ = ∂g

∂Y

Y ̇ +

∂g

∂γ γ ̇ = 0 (9)

the Lagrange (damage) multiplier can be obtained

λ ̇ = −

∂g

∂Y

∂Y

∂ε

∂g

∂Y

∂Y

∂D

∂ f

∂Y + ∂γ∂δ

ε ̇ (10)

where ε is the strain tensor. Then, the damage variable evolution can be calculated from (7). Finally, from (3), (5)

and (8), the hardening variable evolution is

δ ̇ = −λ ̇ (11)

A standard return mapping algorithm with an elastic predictor/damage corrector strategy was used for the numerical implementation into the UMAT user material routine of the commercial FE software Abaqus [4].

3. Damage Interaction

Based on (7) and (3-5), the evolution of one damage component is

D ̇ 1 = λ ̇ ∂gˆ

∂Y = λ ̇2 ˆg

[2H11Y1 + (H12H21)Y2 + (H13H31)Y3 + (H14H41)Y4 + (H15H51)Y5 + (H16H61)Y6] (12)

Thus, a uniaxial stress, which will produce only one component of Y (see (6) and (1)), denoted by Yj, will induce

damage in 1-direction, provided H1 j is non-zero. Thus, to simulate a fracture mechanics type interaction (penny

NASA/TM—2015-218740 3
shaped crack), H11, H14, and H15 can be chosen as non-zero. Then, σ11, σ12, and σ13, would induce damage that

would cause a stiffness reduction in the 11-, 12- and 13-components. Conversely, if the damage evolution of the

material is assumed to be isotropic, then the components Hi j can be chosen as,

H11 = H22 = H33 = Hnn

H14 = H41 = H24 = H42 = H15 = H51 = H35 = H53 = H26 = H62 = H36 = H63 = Hsn

H44 = H55 = H66 = Hss

(13)

In addition to specification of these H matrix terms, c1, c2 and γ0 must be chosen to fully characterize the damage

model.

4. Damage Model Characterization

The damage model was characterized for tension through correlation with a cyclic tension test on a neat epoxy

resin material (RIMR 135/RIMH 137, mass mixing ratio 10:3, 8 hr curing time at 40◦ C). The epoxy was vacuum

injected into a standard dogbone mold to minimize voids. Cyclic tension was achieved using a constant crosshead

speed of 10mm/s, with strains measured via a digital image correlation system.

Figure 1: Model correlation with epoxy tensile test data

The model characterization is shown in Figure 1. Note that the test specimen failed suddenly during the final

reloading cycle at the point indicated. Furthermore, the nonlinear resin response appears to be due to both inelastic

deformation and damage mechanisms as some permanent deformation is observed after unloading, albeit small. In the

present formulation all nonlinear behavior is assumed to be associated with damage only. Since the damage model is

based purely on stiffness reduction, upon unloading, there is zero permanent strain. This is a reasonable approximation

of the actual epoxy tensile response as only a small amount of permanent strain was evident upon unloading. Based

on the correlation shown, the isotropic elastic properties and damage model parameters obtained for the epoxy are

shown in Table 1.

The Poisson’s ratio of the epoxy given in Table 1 was not measured, but set to a typical value. The influence of the

shear damage parameter, Hss, will be examined via parametric studies on a unidirectional composite repeating unit

cell.

NASA/TM—2015-218740 4
Table 1: Material parameters for RIMR 135/RIMH 137 epoxy matrix

Matrix

E 3000 MPa

ν 0.38

c1 27.5 MPa

c2 −2.8

Hnn 0.75

γ0 0.45 MPa

5. Finite Element Models

The three-dimensional damage model presented above has been implemented into Abaqus [4] as a UMAT userdefined material constitutive routine. To represent a unidirectional composite material, three RUC geometries have

been discretized using Abaqus, as shown in Figure 2.

Figure 2: Abaqus finite element model RUCs used to represent a 0.49 fiber volume fraction unidirectional composite. Left: HDP, center: 5 fiber

random, right: 15 fiber random

Figure 2 left shows a microstructure with hexagonal dense packing (HDP) meshed with 2,166 C3D8 elements,

Figure 2 center contains 5 fibers with random locations meshed with 4,420 C3D8 elements, and Figure 2 right contains

15 fibers with random locations meshed with 13,202 C3D8 elements. In all cases, the fiber volume fraction is 0.49.

The fiber is assumed to be carbon; the transversely isotropic elastic properties are shown in Table 2

Table 2: Material parameters for carbon fiber; 1-direction is longitudinal fiber direction

Fiber

E1 290GPa

E2 20GPa

ν12 0.2

G12 20GPa

G23 9GPa

The influence of external boundaries is removed by imposing periodic boundary conditions. Within Abaqus, this

is accomplished using constraint equations, as follows for a parallelepiped model. First, node sets for every node

on each face, edge, and vertex must be generated. The nodes at the vertices must not be included in the node sets

NASA/TM—2015-218740 5
associated with edges or faces. Similarly, the nodes at the edges are excluded from the face node sets. Next, the node

pairs (two nodes that are coupled via a constraint equation) for the three sets of faces, and the node groups for the

three sets of edges, must be identified. The node pair for a face node is always located on the opposite face. For

example, referring to Figure 3, a node on face 4 with the coordinates (x = x1, y = 0,z = z1) has its partner node on

face 2 with (x = x1, y = b,z = z1), where b is the length of the RUC in the y-direction. Similarly, the node groups

for edges must be determined. The edges are categorized based on direction; for instance, the x-direction edge set

includes edges 9, 10, 11 and 12. Then the node groups within each set of edges are found at the same coordinate along

the edge direction. For example, a node on edge 9, with the coordinates (x = x1, y = 0,z = 0) is grouped with nodes

on edges 10 (x = x1, y = b,z = 0), 11 (x = x1, y = b,z = c) and 12 (x = x1, y = 0,z = c), where c is the length of the

RUC in the z-direction.

Figure 3: Reference numbering for application of periodic boundary conditions on a parallelepiped finite element model using equation constraints

within Abaqus. Left: faces, center: vertices, right: edges

In Abaqus, when coupling nodes via equation constraints, it is only possible to couple the displacements (as

opposed to forces). Using three reference points (RPs), designated as RPX at the position of vertex 1, RPY at the

position of vertex 6, and RPZ at the position of vertex 8, the complete (arbitrary) far field strain state can be specified

by applying displacement boundary conditions to these three RPs.

The equation constraints must be defined such that, for every node pair/group, only one degree of freedom remains

in the system, without ever eliminating the same degree of freedom more than once. Denoting the displacement in

coordinate direction i and vertex, edge, face, or reference point j as uj

i , one possibility that accomplishes this is

For edges:

u1

i − u4i − uRPXi = 0

u2

i − u4i − uRPXi − uRPYi = 0

u3

i − u4

i − uRPYi = 0

u5

i − u8i − uRPX

i = 0

u6

i − u8i − uRPX

i − uRPZi = 0

u7

i − u8

i − uRPZ

i = 0

u10

i − u9

i − uRPX

i = 0

u11

i − u9

i − uRPX

i − uRPYi = 0

u12

i − u9

i − uRPY

i = 0

For faces:

NASA/TM—2015-218740 6
u1

i − u3

i − uRPX

i = 0

u2

i − u4i − uRPYi = 0

u6

i − u5i − uRPZ

i = 0

For vertices:

u1

i − uRPXi = 0

u2

i − uRPXi − uRPYi = 0

u3

i − uRPX

i − uRPYi − uRPZi = 0

u4

i − uRPXi − uRPZi = 0

u5

i = 0

u6

i − uRPY

i = 0

u2

i − uRPYi − uRPZi = 0

u8

i − uRPZ

i = 0

where these equation constraints are imposed on paired and grouped nodes, as described above. Thus, in order

to apply a far-field strain, only the displacements of the reference points (RPs), according to the strain state, must be

specified as boundary conditions. Figure 4 shows several examples of deformed RUC plots, where the above periodic

boundary condition application strategy has been employed. As can be seen, the approach maintains periodicity, even

in the presence of significant damage localization.

Figure 4: Example deformed plots (deformation scale factor = 5) showing that periodicity is maintained, even in the presence of significant damage

localization

It is well-known that, when applying damage models such as this in finite element analysis, pathological meshdependence will be present once localization occurs [21, 22]. Thus, for each RUC geometry, it was desired to keep

the element size constant during the study. Further, the global element size variation was kept small by restricting

the element volume during the discretization process in Abaqus. The element volume and was set to the same value

as used in the material calibration step. In the future, length scaling will be added to the model to ensure consistent

energy dissipation in the presence of varying mesh sizes [24].

NASA/TM—2015-218740 7
6. Results and Discussion

The impact of the RUC geometry (see Figure 2) on the predicted local and global response of the composite

has been examined while also varying the value of the shear damage parameter Hss. For the cases examined, the

shear-normal interaction damage parameters,Hsn, has been set to zero, eliminating local shear stiffness reduction due

to normal stresses as well as normal stress reduction due to shear stresses. Figure 5 shows the predicted transverse

tensile response of the composite as a function of the RUC geometry (see Figure 2) and the direction (x or y) of the

applied loading for the case of Hss = 5. Maximum principal strain fields, which highlight the location of damage

within the RUC, are plotted in Figure 6 for five values of the applied global strain. Note that these global strain values

are not exact since adaptive time steps were employed in the finite element simulations.

Figure 5: Comparison of the predicted effective composite transverse tensile response as a function of RUC geometry and direction of loading

(solid line = loading applied in x-direction, dashed line = loading in the y direction) for the case of Hss = 5

While all cases predict nearly identical initial elastic stiffness, once damage initiates, there are significant differences in the global deformation (i.e., damage) response due to microstructural variation in the three RUC representations. For the HDP RUC, there is a large ( 22%) difference between tensile loading in the x-direction and the

y-direction in terms of the stress at which departure from linearity occurs and at which a stress drop occurs. This stress

drop is associated with the onset of damage localization within the RUC. As evident in Figure 6, the local fields are

different (even initially), which leads to initiation of damage localization much earlier when the HDP RUC is loaded

in the x-direction. This is in spite of the fact that, elastically, the RUC behaves transversely isotropically, implying

that the average elastic fields are identical for tension in either direction. The final damage patterns are also quite

different, with localization progressing completely through the RUC much earlier in the case of y-direction loading as

compared to loading in the x-direction.

For the 5 fiber random RUC, Figure 5 shows that, while the difference in the x- and y-direction tensile curve is

reduced compared to the HDP RUC, it is still significant ( 18%). The y-direction loading results in earlier onset of

damage localization because of a concentration that arises in the matrix between two fibers that are close together in

the y-direction, as can be seen in Figure 6. Both 5 fiber RUC loading simulations exhibit earlier damage localization

compared to the HDP RUC simulations due to such concentrations between nearby fibers. In contrast, both 15 fiber

RUC loading simulations predict localization at nearly the same time (strain). The damage patterns shown in Figure

6 are qualitatively similar for each level of applied global strain, with the x-direction loading exhibiting a somewhat

greater extent of damage localization at the applied global strain of 1%. This is manifested in the somewhat more

compliant response after damage initiation in Figure 5. The fact that the x- and y-direction responses are becoming

NASA/TM—2015-218740 8
closer indicates that, as one would expect, as a larger RUC with more fibers is considered, the model provides a better

representation of the actual (transversely isotropic) material response.

Figure 7 shows the predicted x-direction tensile response of the 15 fiber RUC as a function of the value of the

shear damage parameter, Hss, while the corresponding maximum principal strain fields are plotted (at an applied

global strain level of approximately 1%) in Figure 8. It is clear from both figures that, by increasing Hss, which

causes an increase in the shear damage induced by local shear stresses, the progression of damage is accelerated in

the composite. This highlights the importance of capturing normal-shear coupling in the micromechanics solution in

the presence of damage and the value of shear test data for the neat resin.

Figure 9 and Figure 10 show the response and local maximum principal strain plots for applied global longitudinal

shear strains, γxz and γyz. All cases show good agreement in the initial linear regime, with the 5 fiber RUC loading in

yz shear exhibiting slightly higher stiffness.

In general, the effect of the orientation of the applied shear strain (xz vs. yz) has much less of an effect on

the longitudinal shear response compared to the effect of the orientation of the tensile loading (x vs. y) on the tensile

response. The effect is, however, still present, and, once again, it is minimal in the case of the 15 fiber RUC, suggesting

the need to model a good number of fibers within an RUC (to approach a converged representative volume element

(RVE)). For the HDP RUC, while localization initiates at roughly the same time for both xz and yz shear, the damage

progresses much more quickly under yz shear loading. It is clear from the damage localization patterns in Figure 10

that this is due to the more tortuous path of the damage in the case of xz loading, which leads to the tougher response

in Figure 9. In the 5 fiber RUC, the opposite trend is present, with the damage progressing more rapidly, with a

straighter path, under xz shear loading. The effect of the RUC geometry itself on the longitudinal shear response is

also still present, with nearby fibers leading to earlier damage initiation for the 5 fiber and 15 fiber RUCs. However,

the RUC effect is muted compared to the transverse tensile response.

Figure 11 and Figure 12 show the effect of the matrix shear damage parameter, Hss, on the predicted xz shear

response of the 15 fiber RUC. For Hss = 0, no shear damage can occur, and the effective response remains linear to

very high strains when damage due to local normal strains can arise.

As one would expect, increasing Hss, which increases the rate of local matrix shear stiffness reduction in response

to local shear stress, has a primary effect on the predicted composite longitudinal shear response. The effect is much

greater than the effect of Hss on the predicted transverse tensile response (see Figure 7). Figure 13 and Figure 14 show

the predicted transverse (xy) shear response of the composite, along with plots of the local maximum principal strain

fields. The 5 fiber RUC predicts slightly higher elastic transverse shear stiffness compared to the HDP and 15 fiber

RUCs, which are in close agreement in the linear region. The magnitude of the effect of the RUC geometry on the

predicted transverse shear stress-strain response of the composite is similar to its effect on the predicted longitudinal

shear stress-strain response. As before, local stress concentrations arising in the matrix between nearby fibers tend

to cause early onset of localization due to damage. Finally, all three RVE geometries give rise to straight damage

localization paths, leading to similar post-peak xy stress-strain responses.

7. Conclusion

A three-dimensional, anisotropic, progressive damage model, relying on six scalar damage variables, has been

developed following the approach outlined by Barbero [1] in the case of plane stress problems. By extending the

model to three-dimensional stress states, it is now capable of modeling the behavior of the constituents within a

composite, wherein the stress state is inherently three-dimensional. Further, because it is anisotropic, the model could

also be applied to any effective anisotropic material; for example the tows within woven composites (though not done

herein), which also experience three-dimensional stress states. The damage model is also arbitrary in terms of the

interaction among the six damage variables and the stress components, which is controlled by the chosen form of

the damage potential. The damage model has been implemented within Abaqus [4] as a user-defined constitutive

model. The model was characterized for a neat resin material and then applied to model the transverse tensile,

longitudinal shear, and transverse shear response of a unidirectional composite composed of that matrix and carbon

fibers. Finite element models representing RUC geometries with HDP, 5 randomly placed fibers, and 15 randomly

placed fibers were considered. Periodicity conditions were imposed on the boundaries of the finite element RUC

model to ensure representation of the composite material. Results showed that the damage model is capable of

NASA/TM—2015-218740 9
predicting damage localization within a periodic composite in response to both normal and shear loadings for all

three RUCs. Further, while the linear portion of the predicted stress-strain curves of all three repeating cells were

nearly coincident, significant differences were present once damage initiated due to variations in the represented

microstructures. The HDP and 5 fiber RUCs also exhibited significant differences based on the directionality of the

applied transverse tension and longitudinal shear loading. In contrast, the 15 fiber RUC showed only minor influence

of the loading direction. The influence of the shear damage parameter on the predictions of the 15 fiber RUC was also

investigated and shown to significantly soften the post damage initiation response for both applied transverse normal

and longitudinal shear loading.

[1] Barbero, E.J.. Finite element analysis of composite materials using abaqus. Boca Raton, FL 33487-2742: CRC Press Taylor & Francis

Group; 2013.

[2] Kachanov, L.M.. Introduction to continuum damage mechanics. Distributor Center P.O. Box 322, 3300 AH Dordrecht, NL: Kluwer

Academic Publishers; 1986.

[3] Kachanov, LM. On the time to failure under creep conditions, Izv. AN SSSR, Ot. Tekhn. Nauk 1985;8:26–31.

[4] Karlsson and Sorensen. Abaqus: Theory Manual. Providence, R.I: Hibbitt, Karlsson & Sorensen; 1992.

[5] Lemaitre, Jean and Chaboche, Jean-Louis. Mechanics of solid materials. The Edinburgh building, Cambridge CB2 8RU, UK: Cambridge

university press; 1990.

[6] Chaboche, JL and Kruch, S and Maire, JF and Pottier, T. Towards a micromechanics based inelastic and damage modeling of composites.

International Journal of Plasticity 2001;17:411–439.

[7] Krajcinovic, Dusan. Damage mechanics. : Elsevier; 1996.

[8] Voyiadjis, George Z. and Kattan, Peter I.. Damage Mechanics. Boca Raton, FL 33487-2742: CRC Press Taylor & Francis Group; 2005.

[9] Talreja, R. and DCAMM. and DTH.. A Continuum Mechanics Characterization of Damage in Composite Materials. Lyngby. Technical

university of Denmark, Danish center for applied mathematics and mechanics. Report: Technical University of Denmark, DCAMM; 1983.

[10] Talreja, Ramesh. Transverse cracking and stiffness reduction in composite laminates. Journal of composite materials 1985;19:355–375.

[11] Pineda, Evan J and Waas, Anthony M and Bednarcyk, Brett A and Collier, Craig S and Yarrington, Phillip W. Progressive damage and failure

modeling in notched laminated fiber reinforced composites. International journal of fracture 2009;158:125–143.

[12] Pineda, Evan J and Waas, Anthony M. Numerical implementation of a multiple-ISV thermodynamically-based work potential theory for

modeling progressive damage and failure in fiber-reinforced laminates. International journal of fracture 2013;182:93–122.

[13] Bednarcyk, Brett A and Aboudi, Jacob and Arnold, Steven M. Micromechanics modeling of composites subjected to multiaxial progressive

damage in the constituents. AIAA journal 2010;48:1367–1378.

[14] Ladeveze, P., Lubineau, G., Violeau, D.. A computational damage micromodel of laminated composites. International Journal of Fracture  ́

2006;137:139 – 150.

[15] Ladeveze, P and Lubineau, G. On a damage mesomodel for laminates: micro–meso relationships, possibilities and limits. Composites `

Science and Technology 2001;61:2149–2158.

[16] Melro, AR and Camanho, PP and Pires, FM Andrade and Pinho, ST. Micromechanical analysis of polymer composites reinforced by

unidirectional fibres: Part II–Micromechanical analyses. International Journal of Solids and Structures 2013;50:1906–1915.

[17] Melro, AR and Camanho, PP and Pires, FM Andrade and Pinho, ST. Micromechanical analysis of polymer composites reinforced by

unidirectional fibres: Part I–Constitutive modelling. International Journal of Solids and Structures 2013;50:1897–1905.

[18] Kruch, Serge and Chaboche, Jean-Louis. Multi-scale analysis in elasto-viscoplasticity coupled with damage. International Journal of Plasticity

2011;27:2026–2039.

[19] Chen, Yuli and Ghosh, Somnath. Micromechanical analysis of strain rate-dependent deformation and failure in composite microstructures

under dynamic loading conditions. International Journal of Plasticity 2012;32:218–247.

[20] Hansen, NR and Schreyer, HL. A thermodynamically consistent framework for theories of elastoplasticity coupled with damage. International

Journal of Solids and Structures 1994;31:359–389.

[21] Bazant, Zdenek P and Cedolin, Luigi. Blunt crack band propagation in finite element analysis. Journal of the Engineering Mechanics Division

1979;105:297–315.

[22] Pietruszczak, ST and Mroz, Z. Finite element analysis of deformation of strain-softening materials. International Journal for Numerical

Methods in Engineering 1981;17:327–334.

[23] Bazant, Zden ˇ ek P and Oh, Byung H. Crack band theory for fracture of concrete. Mat ˇ eriaux et construction 1983;16:155–177.  ́

[24] Pineda, Evan J and Bednarcyk, Brett A and Waas, Anthony M and Arnold, Steven M. Progressive failure of a unidirectional fiberreinforced composite using the method of cells: discretization objective computational results. International Journal of Solids and Structures

2013;50:1203–1216.

NASA/TM—2015-218740 10

References 
Figure 6: Predicted local maximum principal strain plots during simulated transverse tensile tests at different applied global strains  ̄ε as a function

of RUC geometry and direction of loading for the case of Hss = 5. The deformation scale factor is 5

NASA/TM—2015-218740 11
Figure 7: Comparison of the predicted effective composite X-direction transverse tensile response for the 15 fiber RUC as a function of the shear

damage parameter, Hss

Figure 8: Predicted local maximum principal strain plots during simulated x-direction transverse tensile tests on the 15 fiber RUC as a function of

the shear damage parameter, Hss, at an approximate applied global strain of  ̄ε = 1%. The deformation scale factor is 5

NASA/TM—2015-218740 12
Figure 9: Comparison of the predicted effective composite longitudinal shear response as a function of RUC geometry and direction of loading for

the case of Hss = 5

NASA/TM—2015-218740 13
Figure 10: Predicted local maximum principal strain plots during simulated longitudinal shear tests at different applied global shear strains  ̄γ as a

function of RUC geometry and direction of loading for the case of Hss = 5. The deformation scale factor is 5

NASA/TM—2015-218740 14
Figure 11: Comparison of the predicted effective composite xz longitudinal response for the 15 fiber RUC as a function of the shear damage

parameter, Hss

Figure 12: Predicted local maximum principal strain plots during simulated xz longitudinal shear tests on the 15 fiber RUC as a function of the

shear damage parameter, Hss, at an approximate applied global shear strain of  ̄ gamma = 1.35%. The deformation scale factor is 5

NASA/TM—2015-218740 15
Figure 13: Comparison of the predicted effective composite transverse (xy) shear response as a function of RUC geometry for the case of Hss = 5

NASA/TM—2015-218740 16
Figure 14: Predicted local maximum principal strain plots during simulated transverse (xy) shear tests at different applied global shear strains  ̄γ as

a function of RUC geometry and direction of loading for the case of Hss = 5. The deformation scale factor is 5

NASA/TM—2015-218740 17