Use of Volume Average Stresses for Predicting Failure

FEA Final Project Report

Introduction

Modern composites are not only lightweight and strong but also chemical and corrosion resistant giving them an obvious edge over conventional materials like metals. They are widely used in many industries including aerospace, automobile, marine and sports. Thus accurate estimation of composite properties like yield strength, fracture strength and cycles to failure has become important. The prediction of failure of composites under arbitrary load states for arbitrary laminates is a challenging problem. One may extract volume average stresses of the composite and apply failure criteria to the constituents directly [1][2][3][4]. It is known that the strain energy computed using volume average quantities may not account for all the strain energy in the composite. This inequality arises due to the inhomogeneous material properties of the composite thus having a direct implication on the accuracy of the failure theories that depend on volume average quantities to predict failure. Failure predicting capabilities of such theories may be further improved by augmenting them with this missing energy which is termed as interaction energy. The focus of this work is to study the variation of interaction energy and its dependence on ber volume fraction. Another aspect of this work is to examine the contributions of composite constituents (ber and matrix) to the interaction energy.

Theoretical Motivation

Consider a general composite material consisting of a ber and matrix phase. Let denote the strain energy of the composite under an arbitrary load state. Assuming linear elasticity, the strain energy of the composite can be represented mathematically as [5] = 1 2
V

T c

c dV

1 T 2 c

Vc

(1)

where the terms in the brackets denote volume average quantities and Vc is the volume of the composite. Equation (1) may be separated into contributions from the ber and matrix as = f + m = 1 2
Vf

(2) T m
m dVm

T f

f dVf

1 2
Vm

(3)

When dealing with constituent average stresses and strains, all of the strain energy is not accounted for. Denoting this unaccounted energy by , eqs.(3) can be rewritten as f + m = 1 T 2 f
f

Vf +

1 T 2 m

Vm +

(4) (5)

= f Vf + m Vm

where f and m are the unaccounted constituent energies of the ber and the matrix respectively and is the interaction energy. The inhomogeneties of the material properties give rise to uctuation in strain . This may be represented as which in turn give rise to uctuation in stress f =
f

f = f f 1

(6)

m =

m = m m

(7)

Substituting eqs.(6) and (7) in eq.(3) yields the expression for f and m as f = 1 2
Vf

T f f dVf

or

f =

1 2
Vf

ij f ij f dVf

(8)

m =

1 2
Vm

T m m dVm

or

m =

1 2
Vm

ij m ij m dVm

(9)

Considering the indicial form of eqs.(8) and (9), using Hookes law to write stress in terms of strains and then expanding we get f = Af + Bf + Cf (10) where
2 2 Af = C11f 2 11f + C22f 22f + C33f 33f

Bf = 2C12f 11f 22f + 2C13f 11f 33f + 2C23f 22f 33f Cf = C44f Similarly for the matrix m = Am + Bm + Cm where
2 2 Am = C11m 2 11m + C22m 22m + C33m 33m 2 12 f

(11)

+ C55f

2 13 f

+ C66f

2 23 f

(12)

Bm = 2C12m 11m 22m + 2C13m 11m 33m + 2C23m 22m 33m Cm = C44m
2 12 m

(13)

+ C55m

2 13 m

+ C66m

2 23 m

3
3.1

Modeling
Overview

The goal of this work is to investigate the behavior of interaction energy with ber volume fraction. This was established by building and simulating a model in the commercially available FEA software ABAQUS. For the purpose of analysis, a representative volume element (RVE) having hexagonal ber packing was considered. The RVE has one entire ber section at its center and one quarter section of the ber at each vertex. The rest of the RVE consists of the matrix. Refer g.1 in appendix B.

3.2

Materials

The RVE contains a single type of ber and matrix. The ber material is carbon and the matrix material is epoxy. Refer Table 1 for detailed material properties. Table 1. Baseline material properties Material Carbon Epoxy Type Othrotropic Isotropic E1 (GPa) 235 4.8 E2 (GPa) 14 4.8 G12 (GPa) 28 1.8 12 0.2 0.34 23 0.25 0.34

3.3

Type of elements and meshing

For all the models, a single type of element C3D8R (an 8-node linear brick, reduced integration, hourglass control) was used. The RVE had to be partitioned appropriately to achieve a perfectly symmetric mesh. For the ber sections, sweep mesh technique was used and the medial axis algorithm was employed to reduce the mesh transition, while for the matrix sections, structured mesh technique was used. For maximum accuracy and to reduce computation time, the mesh size was varied with volume fraction. The mesh is comparatively coarse at lower ber volume fractions and becomes ne at higher ber volume fractions.

3.4

Boundary conditions

Periodic boundary conditions were applied to the RVE such that respective nodes on opposite vertices, edges and faces have equal displacements in any direction. This was established this by extracting all the nodes from the RVE and separating them into three ordered sets containing nodes on the vertices, nodes on the edges and nodes on the faces. It must be noted that, the edge node set does not contain the nodes on the vertices and similarly, the face node set neither contains the nodes on the vertices nor those on the edges. The nodes from these sets were then reordered appropriately and stored in unordered sets. Equations were written on the unordered sets to enforce periodicity. Six loads were applied to the RVE in separate steps. These loads were in the form of strains, which were achieved by providing displacements to control nodes which render movement to faces in three directions. Table 2 in appendix A lists the prescribed strains. Table 2. Prescribed strains Load Strain
xx xx xx xx xx xx

0.01

0.01

0.01

0.01

0.01

0.01

Results and discussions

For the RVE with hexagonal ber packing, the ber volume fraction was varied from 0.05 to 0.85. Volume average quantities (stresses, strains and stiness) of the composite were extracted from the model and then the total strain energy c was calculated using eq.(1). Similarly, volume average quantities (stresses, strains and stiness) were extracted both from the ber and the matrix sections and strain uctuations were computed using eqs.(6) and (7). Substituting these strain uctuations and the extracted stiness in eqs.(11) and (13) and using eqs. (10) and (12) yields f and m . The interaction energy can be calculated using eq.(5). The interaction energy for all the load cases was computed individually. This procedure was repeated for various ber volume fractions of 0.05 to 0.85. The interaction energy fraction / was calculated for each load case and for all the ber volume fractions. Plot 1 shows / as a function of ber volume fraction. It can be seen that interaction energy reaches it peak at a ber volume fraction of 0.7 and is apparent that at the peak, for the load case shear-xy around 27% of the strain energy is unaccounted for by using just the volume average quantities. Another important trend that can be noticed is that the interaction energy rst increases with ber volume fraction, reaches its peak and then decreases. This behavior can be better understood by analyzing eqs.(10) through (13). Since the stinesses of both the constituents remain constant, the constituent interaction energies are proportional only to the strain uctuations in the constituents. These strain uctuations give rise to stress uctuations in the constituents thereby giving rise to interaction energy. A homogeneous material does not have any strain uctuation and so does not have any interaction energy. At a ber volume fraction of 0, the composite material is essentially composed of all matrix and at a ber volume fraction of 1, the composite material is composed of all ber. Thus at these ber volume fractions the interaction energy of the composite must zero. This explains why the interaction energy rst increases and then decreases with ber volume fraction.

The interaction energy is always greatest for the load case of shear-xy and negligible for load case of tension-xx for all ber volume fractions. This is because strain uctuations, which give rise to stress uctuations, are highest in load case shear-xy and negligible for load case tension-xx. This can be seen from gures 2 and 3 in appendix B. In the remaining load cases, the uctuations in strain (and stress) are less than shear-xy but more than tension-xx. As a result of this, the interaction energy for these load cases is always less than that in shear-xy and more than than tension-xx for all ber volume fractions. This can be seen from gures 4 and 5 in appendix B.

4.1

Contribution of matrix to the interaction energy

Finally, the contributions of matrix to the interaction energy for all the ber volume fractions was computed. Plot 2 shows the matrix contribution to interaction energy plotted against ber volume fractions. It can be observed that up to reasonable ber volume fractions, the matrix accounts for the bulk of the interaction energy.

Summary and conslusions

In this study, an expression for interaction energy was obtained and a parametric study was performed in which the ber volume fraction was varied and the interaction energy was calculated. The nature and dependence of the interaction energy on ber volume fraction was examined. Finally, the matrix contribution to interaction energy was calculated. It was concluded that interaction energy peaks around a common ber volume fraction. It is minimum for composite loading in the ber direction because of minimal strain uctuations and maximum for composite shear loading in the transverse direction because of maximum strain uctuations. The matrix is the dominant contributor to interaction energy irrespective of ber volume fraction. These results are of vital importance when considering failure in composite materials, and that, the interaction energy in the composite must be accounted for. In order to get an accurate estimate of failure we need to augment constituent-level failure theories with this interaction energy. Since it can be seen from the results, that the matrix is the major contributor to this energy, we may ignore the ber almost entirely (for typical carbon-epoxy systems) and concentrate our eorts on augmenting the matrix failure theory to better predict failure.

References
[1] Ray S. Fertig III, An Accurate and Ecient Method for Constituent-Based Progressive Failure Modeling of a Woven Composite. Minerals, Metals and Materials Society, Feb 2010 [2] Ray S. Fertig, III and Douglas J. Kenik, Predicting Composite Fatigue Life Using Constituent-Level Physics. Firehole Composites, Laramie, WY, 82070 [3] Ray S. Fertig, III, Bridging the gap between physics and large-scale structural analysis: a novel method for fatigue life prediction of composites. Firehole Composites, Laramie, WY, 82070 [4] Ray S. Fertig, III, A Computationally Ecient Method For Multiscale Modeling Of Composite Materials: Extending Multicontinuum Theory To Complex 3D Composites. Firehole Composites, Laramie, WY, 82070 [5] C. T. Sun and R. S. Vaidya, Prediction of composite properties, from a representative volume element, Compos. Sci. Technol., vol. 56, no. 2, pp. 171179, 1996.

Appendix A
Table 1. Baseline material properties Material Carbon Epoxy Type Othrotropic Isotropic E1 (GPa) 235 4.8 E2 (GPa) 14 4.8 G12 (GPa) 28 1.8 12 0.2 0.34 23 0.25 0.34

Table 2. Prescribed strains Load Strain

0.01

0.01

0.01

0.01

0.01

0.01

Appendix B
Figure 1. Hexagonal ber packing with a ber volume fraction of [Link] circles in gray are the bers which are unidirectional along the x-axis. They are surrounded by the matrix which is depicted in dark green.

Figure 2. Meshing of the RVE at ber volume fractions of 0.4 and 0.75

Appendix B
Figure 3. Stress plot for Load case tension-xx

Figure 2. Stress plot for Load case shear-xy

Appendix B
Figure 5. Stress plot for Load case tension-yy

Figure 6. Stress plot for Load case shear-yz