Microstructure Lectures

by

Ronald D. Kriz, Associate Professor
Engineering Science and Mechanics
Virginia Polytechnic Institute and State University
Blacksburg, Virginia 24061


Table of Contents:

  1. Introduction
  2. Cracks: Atomistic - at Interfaces - in a Continuum
  3. Cracks Near Interfaces Between Dissimilar Isotropic Materials
  4. Introduction to Mechanical Behavior of Anisotropic Laminates
  5. Laminate Singularities Caused by Anisotropy: "Free-Edge Problem"
  6. Laminate Singularities Caused by Ply Cracks
  7. Cracks Near Interfaces Between Dissimilar Anisotropic Materials

    ------ Extension to a Homogeneous Continuum ------

  8. Cracks in Homogeneous Isotropic Materials
  9. Cracks in Homogeneous Anisotropic Materials
  10. Wave Propagation in Homogeneous Isotropic/Anisotropic Materials
  11. References

1.0 Introduction:

The micro-scale, 10-3m, exists between nano, 10-9m, and macro, 100m. At the micro-scale, individual grains in a polycrystalline material or "lamina" in a fiber-reinforced composite laminate, would be consider typical micro-structures.


Figure 1. Microstructures of grains and lamina

Macroscopic mechanical behavior can be influenced by structures at the microscopic scales of 10-3m. For example elastic properties of a polycrystalline material are controlled by the orientation of individual grains sizes of 10 -3m, i.e. random orientation leads to isotropic macroscopic elastic properties. For composite laminates macroscopic properties of laminated composites are influenced by the elastic properties and fiber orientation in each laminate layer ("lamina"). Typical lamina thicknesses are 10-3m and hence are considered microscopic structures. Cumulative damage events, i.e. ply cracks, at the lamina-microscopic scale can also influence mechanical behavior and failure at the macroscopic laminate scale. Relationships between the microscopic and macroscopic scales embraces a wide variety of other topics and examples. It is not the intention of this section to provide a comprehensive overview of these micro/macro models but rather to introduce the students to examples of micro/macro models that is representative of the research conducted by Dr. Ron Kriz while working with computers at NIST (formerly NBS) from 1980 to 1990 and with computers at NCSA (National Center for Supercomputing Applications) from 1990 to 2000. During this time period high performance computing was rapidly evolving. CDC, Thinking Machine CM2, Cray, and SGI computers were all used to study the mechanical behavior of a variety of problems related to damage in laminated composites and wave propagation in highly anisotropic crystals and fiber-reinforced composites.

Throughout this section we provide links to examples where the reader can work with interactive computer programs, written by Dr. Kriz, as they relate to the theory developed in the sections below. Hence our objective is to create an interactive section that compliments the theory and allows the student to experiment with the various micro/macro models as they are introduced.


2.0 Cracks: Atomistic 10-9m - at Interfaces 10-3m - in a Continuum 100m

Cracks span the scales. The mechanical behavior of materials with cracks can be modeled to include: 1) the influence of individual atoms ("nano"), 2) the presence of interfacial layers separating different microstructures ("micro"), or 3) the idealized continuum model where atomic and interfacial effects are neglected ("macro"). Regardless of the scale, models make simplifying approximations. Understanding these approximations is key to understanding how some models bridge the scales. Therefore we will highlight these approximations below as we develope the theory and also when we develope the mathematical model for each of the scales. The images on the home page exemplify these three scales.


Figure 2. Schematic of Mode-I cracks spanning the length scales

In this section we focus our attention on modeling the mechanical behavior of materials that are influenced by cracks at or near interfacial layers or grain boundaries. The image in the center of Fig. 2 shows a schematic of a crack that has intersected a boundary between two dissimilar materials: region 1 and region 2. This boundary can be a grain boundary or interfacial layer separating two different lamina in a composite laminate. In either case the model for predicting the behavior of crack propagation is fundamentally the same but in the following sections we will focus mostly on interfaces in laminated composites.


3.0 Cracks Near Interfaces Between Dissimilar Isotropic Materials

Numerous models for interfacial cracks exist. In the most general case we assume a crack arrested at an interface separating two dissimilar anisotropic materials and then simplify the model to dissimilar isotropic materials. Such cracks are common in fiber-reinforced composite laminates. In Figure 3 we show a schematic and an acetate replica micro-graph of a 90 degree ply crack arrested at a 0/90 degree ply interface.

Figure 3. Schematic and acetate replica micro-graph showing a ply crack arrested at a 0/90 ply interface

If the crack continues to propagate, it can grow into the adjacent 0o ply, causing the fibers to fail or it can grow into an interfacial crack and debond along the 0/90 interface. Figure 4 shows experimental evidence of both types of crack growth, that can be controlled by the residual stress state, Ref.[1]. Statistical evidence confirms that laminate fracture strength is lower, if the ply crack extends into the 0o load bearing ply, and higher if the ply crack extends along the interface as an interfacial debond. Hence we have experimental evidence that macroscopic laminate fracture strength is controlled by damage at the microscopic lamina scale.

Figure 4. Ply crack debonds along the 0/90 interface or extends into the 0o ply.

With this experimental evidence the problem at hand is to create a model that will scale both the micro and macro and predict crack growth near an interface between two dissimilar materials. With such a model we can then study which parameters control crack growth. Does the residual stress state, anisotropy, or both control how the crack grows near an interface? The most general mathematical formulation is shown in Figure 5.

Figure 5. Schematic showing debond along the interface or self similar crack growth into the adjacent region.

(1)

From elementary fracture mechanics we can idealize the singular stresses near the ply crack tip as a linear superposition of Mode-I, Mode-II, Mode-III stress states. The superscripts and subscripts I, II, and III represent the three different fracture Modes. The exponent -α on the variable "r" (radius) has a value of -1/2 for isotropic and homogeneous anisotropic materials, but differs from -1/2 for dissimilar isotropic materials, which will be discussed below. The constants CI, CII, and CIII are functions of Stress Intensity and radius "r", and the variables FijI, FijII, and FijIII are functions of geometry only, i.e. location of boundaries and the angle, θ, measured at the crack tip. The remainder of this section will focus on deriving simple relationships for the special case of two dissimilar isotropic regions. After we learn more about anisotropy we will return to this problem where the dissimilar regions are anisotropic.

An introduction of the different fracture Modes is provided in classic fracture mechanics publications, Ref.[2,3] and at the Department of Materials Science and Engineering (MSE) Fracture Web Site (Student Team Project MSE2094, Spring 1997), see section on "Energy Methods " and "Stress Intensity"for an introduction to the different fracture Modes.

The first model of a crack near interfaces separating dissimilar isotropic materials was developed by Williams Ref.[4,5].


Interface Crack "Debond" Between Two Dissimilar Insotropic Materials

Williams Ref.[4] was the first to show that the singularity on "r" for an interfacial crack, "Debond" between two dissimilar materials, exhibited an oscillating singularity.

(2)

Derivation of this equation will be outlined and discussed in class. Although the solution is correct mathematically, the strong oscillation predicted as the observer approaches the crack tip is physically unrealistic. Research by Dan Post, Professor at Virginia Tech, experimentally confirm a stress reversal near a debond crack by high resolution Moire Interferometry but although this is an encouraging piece of information it is far from what we need to confirm a full oscillation. Whether or not such a oscillating singularity can exist in a real material and at what scale such an oscillation might exist is exactly the type of question that arises when we try to interpret and justify model predictions. Although we are at the micro level were the crack behavior is controlled by two "thin" dissimilar isotropic materials we still assume a continuum but this approximation breaks down as we approach atomic spacing ("nano"). Are such oscillating singularities a realistic mechanical behavior that can be modeled at the nano scale? Probably not. At this point this type of behavior is at best only an interesting mathematical solution to a interfacial debond solution at the microscale that suggests that this behavior extends into the nanoscale -- maybe.


Crack Normal to Interface, "Ply Crack", Between Two Dissimilar Isotropic Materials

Williams Ref.[5] also showed that the singularity for a crack acting normal to an interface between two dissimilar isotropic regions is a function of the ratio of the elastic shear moduli of these two regions. Here we derive an expression for the exponent, λ on the radius "r". The exponent, λ is calculated as an eigenvalue or root from the characteristic equation below.

(3)

Figure 6. Exponent singularity as a function of the ratio of shear moduli

DERIVATION OF EQUATIONS FOR PREDICTING THE STRESS SINGULARITY, λ


4.1 Introduction to Mechanical Behavior of Anisotropic Materials

Our previous discussion on laminated structures provides many microscale examples of mechanical behavior. Since most composite laminates are highly anisotropic it would be appropriate to first learn more about how anisotropy influences mechanical behavior and then we can proceed with an introduction to the mechanical behavior of anisotropic laminates.

In the early part of the 20th century most engineering materials were isotropic. In the later part of the 20th century, 1970s -> 1990s, the aerospace industry created the demand for new, high-strength, high stiffness and light-weight fiber-reinforced composite materials. Because these materials are also highly anisotropic, designers could optimize their designs, i.e. orient the direction of high stiffness and strength in the load direction. Other materials used in "Smart" or "Intelligent" structures, such as piezioelectric polymers are also highly anisotropic. Because the properties of these new materials could be tailored, designers, could for the first time, design the material for the application instead of designing an application around fixed material properties. Anisotropic materials are becoming the material of choice in a variety of engineering applications as we enter the 21st century. Studying the mechanical behavior of these highly anisotropic materials under static and dynamic loads at all scales is relevant when considering current engineering practice.

In the early part of the 20th century anisotropy was more a topic of scientific research then a property used in engineering design. Nye's text on "Physical Properties of Crystal", Ref.[6], gives an excellent introduction to anisotropy in crystals from a material scientist perspective. Lekhnitski's text on "Theory of Elasticity of an Anisotropic Body", Ref.[7], began to address issues of anisotropy in engineering design. In this section we will give a very brief introduction to anisotropy for a continuum (macroscale) and then apply what we learned to anisotropic laminates at the microscale. In section 10.0 we will also show how anisotropy can significantly control dynamic mechanical behavior in a continuum (macroscale).

Mechanical behavior of anisotropic materials is characterized in the most general notation as a fourth order stiffness tensor, Cijkl, which is contracted into a six by six stiffness matrix, Cα,β, relating six components, σα, of a symmetric stress tensor to the six components, eβ, of the symmetric strain tensor which are written here as column vectors.

(4)

Although this contracted notation shown above is convenient and will be used below to model mechanical behavior of anisotropic laminates at the microscale, a more complete description of elastic anisotropy will require that we return to the more general form where the stiffness matrix that relates strains to stress is written as a fourth order stiffness tensor Cijkl, where the stresses, σij, and Lagrangian strains, lij, must then be written as second order tensors. We will show that this more general tensoral form will allow us to create unique three dimensional geometric representation surfaces for each unique anisotropic elastic symmetry.

(5)

First we will show how stress can be written in the form of a second order tensor. The proof that stress is a second order tensor because it transforms as a second order tensor is given in Ref.[8].

Figure 7 shows how all possible components of stress acting as forces on a small differential element can be organized into a matrix format which can also be reduced to index notation where i and j are called "free" indices that are assigned values of 1, 2, or 3. All possible combinations of the i and j indices yield the desired three by three matrix.


Figure 7. Stress tensor components shown on a differential element

Writing all six components of strain into a similar matrix is not as straight forward as it was for stress where each stress components could be simply seen as a force acting on the surface of a cube. Strain implies displacement which requires much more development, see Ref.[8] Chapter 3, "Analysis of Deformation in a Continuum". Suffice it to say that strain like stress has 9 components which can be identified with the same indices which identifies a displacement direction or shearing in a plane associated with each strain component. From a Lagrangian viewpoint, for large displacement gradients, the nonlinear strain displacement relationships are shown below.

(6)

If we assume only small displacement gradients then the nonlinear terms can be neglected and the upper case "L" is reduced to the lower case "l" which is used to denote the linear lagrangian strain tensor, lij.

(7)

With strain and stress defined as second order tensors, it is not difficult to derive the desired linear anisotropic relationship, equation (5), from the first law of thermodynamics. One form of the first law of thermodynamics is shown below.

(8)

where u is the specific internal energy per unit mass, ρ is the density, eij is the Eulerian strain associated with Eulerian displacement gradients, c is the radiative heat transfer per unit time per unit mass, and bi is the heat transferred through a unit surface per unit time. Note d/dt is the total or comoving derivative. For a complete derivation see Ref.[8], Chapter 4, section 4.5.

If we assume no heat transfer by radiation or conduction and we assume that only small displacement gradients result in lagrangian strains, lij, equation (8) reduces to a form where u is now called the strain energy per unit mass since the only mechanism available for storing internal energy in the continuum is through deformation.

(9)

Next let's assume density is constant, or at least a function of strains and that stresses are related to strains by some unknown function (to be determined). Hence strain energy can be a function of stress or strain. For the purpose of discussion here we assume strain energy can be written as some function of strain.

(10)
hence,

(11)

Substituting equation (11) into equation (9) yields.

(12)

If this relationship is true for any arbitrary nonzero dlij, then the relationship in the parenthesis must vanish which yields an equation that shows the relationship between stress, strain and strain energy per unit mass.

(13)

Recall that for equation (10) we assumed some undefined functional form for strain energy as a function of strain. One possible function would be to expand the strain energy into a power series whose coefficients like strain are tensors but constants with indices that sum ("contraction") with the indices of the strain tensors such that each term in the series is a scalar component of the strain energy. The first term α is an arbitrary energy reference state, so we can set it to zero. The second term is a zero strain reference state. The third term is the linear elastic term and the fourth term is the nonlinear, but elastic, term. Higher order terms can be added here as needed to model the desired constitutive stress-strain relationship.

(14)

Substituting this function for strain energy into equation (13) yields the most general form of the anisotropic elastic constitutive law in tensor format.

(15)

NOTE: 1) it is perhaps more evident at this point that the coefficient βij represent the residual stresses when the strains, lij, go to zero, and 2) the coefficient δijklmn is the nonlinear elastic term that gives rise to nonlinear elastic stress-strain diagrams not associated with dissipative constitutive models, which we also set to zero. With these simplifications equation (15) reduces to the desired linear-elastic equation.

(16)

Where Cijkl are the terms that relate components of strain to the components of stress per unit mass. A similar expression without density can be derived where strain energy per unit volume is introduced into the derivation. Interestingly to remove density is not as straight forward as one might think, see Ref.[8].

(17)

Lets briefly look at an arbitrary geometric transformation of equation (17). For a full development of this type of transformation see Ref.[8]. It can be shown that both stress and strain are second order tensors because they transform as second order tensors and Cijkl transforms as a fourth order tensor.

(18)
(19)
(20)
where, aij is the direction cosine transformation matrix (not a tensor) for any arbitrary rotation.


Transformation Law of a vector


Figure 8. Vector drawn with respect to two different Cartesian axes

Introduce the direction cosines between the two sets of axes:

(a11, a12, a13) = direction cosines of the y'1-axis with respect to the (y1,y2,y3)-axes.
(a21, a22, a23) = direction cosines of the y'2-axis with respect to the (y1,y2,y3)-axes.
(a31, a32, a33) = direction cosines of the y'3-axis with respect to the (y1,y2,y3)-axes.


Figure 9. Direction Cosine Matrix With a Simple Rotation About the y3 axis.

Hence, the relation between the vector A' and A can be summarized as follows.

(21)

Note that the indices on A' can be replace as i=1,2,3 hence the three equations in equation (21) reduces to one equation with a "free" index "i".

(22)

Also note that when the number in each term in equation (22) are the same (1 1 in the first term, + 2 2 in the second term, and + 3 3 in the third term) we introduce a new idea of a "summation" notation where any such summation can be notationally simplified and replaced with repeated indices in one term. Hence we can now use this summation index to reduce equation (22).

(23)


The same direction cosine matrix, aij, can also be used in the transformation of stress, strain, and stiffness, in equations (18), (19), and (20) which when substituted into equation (17) reveals that the same equation exists in the arbitrary prime coordinate system, equation (24). The fact that equations (17) and (24) are identical demonstrates that equation (17) is invariant to any ARBITRARY transformation, hence equation (17) is a "tensor equation". The keyword here is arbitrary.

(24)

NOTE ON ARIBITRARINESS: Before we move on to understanding the different types of anisotropies, there is one very important idea that needs to be introduced about equations (17) and (24) which is that tensor equations by definition exhibits mathematical "invariance" to any ARBITRARY transformation. Although the two ideas of invariance and arbitrariness are mathematical ideas, these same ideas of invariance to an arbitrary transformation is also the same fundamental idea, indeed a requirement, of a physical law. Hence invariant arbitrary transformations are essential properties for both physical laws and tensor equations. Are all tensor equations by definition physical laws, is a question that remains unanswered but an interesting topic for discussion.

Equations (17) and (24) are the form of Hooke's Law we need to continue the discussion of anisotropic elasticity by Nye and also to show some interesting "glyphs" (geometric representations) of the fourth order stiffness tensor, Cijkl .

Now that we have introduced elastic anisotropy both as a fourth order tensor and a six-by-six matrix with reduced notation, the reader can now understand the comprehensive treatment of elastic anisotropy given by Nye, Ref.[6], in chapter 8, "Elasticity. Fourth-Rank Tensors", pp. 131-149, where both stiffnesses and compliances are described in terms of matrix and tensor notation. Our objective here is to summarize and highlight the different elastic anisotropies commonly encountered by todays materials engineers. Figure 10 summarizes the various elastic anisotropies in matrix notation.

KEY TO NOTATION
TRICLINIC (21)
MONOCLINIC (13)
ORTHORHOMBIC (9)                       CUBIC (3)    
(7)          TETRAGONAL          (6)
(7)             TRIGONAL             (6)
HEXAGONAL (5)                          ISOTROPIC (2)

Figure 10. Various anisotropies for Sij and Cij.

As shown in Figure 10, the different elastic anisotropies are related to the number of independent terms in the elastic property matrix where symbols represent either stiffnesses or compliances. Each of these unique matrices can be related to crystal symmetry classes but here we will discuss just a few of the major classes and symmetry conditions. Triclinic is the most general with 21 independent components. Starting with Triclinic the number of independent terms can be reduced using symmetry planes and invariant rotations. For example if we could imagine ourselves standing inside a crystal and observed that the crystal structure looked identical if we looked either in the +z or -z direction. With this observation we would conclude that symmetry exists about the x-y (1-2) plane and the transformation matrix, aij, would be the identity matrix but with the term a33 set equal to -1, see Figure 11. Using this transformation matrix in equation (17) components of the stress tensor are compared in the primed and unprimed coordinates.


Figure 11. Sample calculation of transformation for symmetry about the x-y (1-2) plane.
NOTE: transformation of stress and strain requires the use of tensor transformations
eqns. (18) and (19) to get the correct sign.

This comparison determines which matrix terms must vanish. For example in Figure 11 comparing σ1 with σ1' requires C14 and C15 vanish. Similar comparisons with the other stress components requires,

(25)

Hence this reduces the Triclinic with 21 independent terms to Monoclinic with 13 independent terms. Figure 10 shows two different types of Monoclinic matrices. If we continue with symmetry in the x-z and y-z plane additional terms must vanish and Monoclinic reduces to Orthotropic with 9 independent terms. Arbitrary rotation about the z axis (3-direction) assuming elastic properties are invariant reduces Orthorhombic to Hexagonal with five independent terms. Of course for isotropic there are only two independent terms or components of the stiffness or compliance fourth order tensor.

Elastic symmetries commonly encountered when working with fiber-reinforced materials commonly used in lightweight aircraft and aerospace structures, are called "orthotropic" which is the same as orthorhombic and "transversely isotropic" which is equivalent to hexagonal.

We continue the development of anisotropic mechanical behavior into two sections. First we focus on engineering elastic properties and represent anisotropy in terms of engineering properties: E, G, and ν, and second we focus on all components of the fourth order stiffness tensor. In both cases we ephasize the use of geometry to study the anisotropy.

Geometric Representations of Elastic Anisotropy: Single Stiffness Components: E, G, and ν

Another linear constitutive relationship, similar to equation (24), exists for compliances, Sijkl, that relate stresses with strains.

(26)

Compliances like stiffnesses can be considered fourth order tensors because they transform as fourth order tensors.

(27)

Hence the relationship between compliances and stiffness can be shown as an inverse relationship. This inverse relationship is shown here in reduced contracted notation. NOTE: although these matrices are inversely related, Sij and Cij are not tensors quantities because they do not transform as second order tensors.

(28)

Compliances, unlike stiffnesses, have a simple relationship with engineering elastic properties such as Young's modulus, shear modulus, and Poisson's ratio. For example the compliance along the z axis (3-direction), S33 is defined as the reciprocal of the Young's modulus in the same z-direction.

(29)

Similarly expressions for shear modulus and Poisson's ratio in the y-z (2-3) plane can be simply derived from terms in the compliance matrix.

,       (30), (31)

It is now possible to show how these three engineering elastic properties vary as a function of rotation from the z axes (3-direction) into the x-y (1-2) plane by using the general fourth order compliance tensor transformation of equation (27) where the direction cosine terms, aij, are expanded into sines and cosines of the rotation angle, θ, starting at zero along the z axis and rotating 90o into the x-y plane. This transformation although simply done as a fourth order transformation, becomes a tedious task in reduced contracted notation. After much algebra the final equations are simplified, see Ref.[9].

(32)

(33)

(34)

For all isotropic engineering materials equations (32), (33), and (34) will, by definition, reveal no variation in properties as we rotate from the z axis into the x-y plane, hence circles in a two-dimensional (2-D) plane or spheres in a three-dimensional (3-D) space are visual representations (manifestations) of isotropic elastic properties. Unidirectional fiber-reinforced composite materials, on the other hand, exhibit highly anisotropic behavior which is exhibited as significant deviations from the isotropic spherical geometry.

Table 1. Typical carbon-fiber/epoxy-matrix composite stiffness, GPa, Ref.[10]

Material C33 C11 C12 C13 C44 ν32 ν12
Epoxy-matrix (Kriz-Stinchcomb, 1979) 8.63 8.63 4.73 4.73 1.95 0.345 0.345
Carbon-fiber (Kriz-Stinchcomb, 1979) 235 20.0 9.98 6.45 24.0 0.279 0.490
Carbon-fiber (Dean-Truner, 1973) 240 20.4 9.40 10.5 24.0 0.35 0.497
Carbon-fiber (Smith, 1972) 221 19.4 6.60 5.80 20.3 0.230 0.320
Composite (60% Fiber) (Kriz-Stinchcomb), 1979) 144 13.6 7.0 5.47 6.01 0.284 0.497



Table 2. Stiffnesses in Table 1 converted to engineering properties, GPa, Ref.[10]

Material Poisson's Ratio(ν)    Modulus: Bulk(K), Shear(G), Young's(E)
Subscripts: Matrix(m), Fiber-Axis(L), Transverse-Plane(T)
Isotropic: two independent properties-> Em νm
Epoxy-matrix (Kriz-Stinchcomb, 1979) 5.35 0.345
Hexagonal: five independent properties-> EL ET GLT GTT KTT
Carbon-fiber (Kriz-Stinchcomb, 1979) 232 15.0 24.0 5.02 15.0

For example if we take elastic properties predicted in Ref.[10], see Table 1., for unidirectional carbon-fiber/epoxy composites and substitute these properties into equations (32), (33), and (34), both spherical and nonspherical geometries are predicted, see Figures 12, 13, and 14. In these figures the elastic properties are also shown as a function of fiber volume fraction where 0.0 fiber volume fraction yields pure isotropic epoxy properties, shown as circles, and 1.0 fiber volume yields pure anisotropic fiber properties and hence the geometry significantly deviates from the circular geometry. Even low fiber volume fractions yield significant anisotropy.

Click on image to start an interactive JWave session

Figure 12. Polar diagram of graphite/epoxy Young's modulus E3 for various fiber volume fractions, in GPa.

Click on image to start an interactive JWave session

Figure 13. Polar diagram of graphite/epoxy shear modulus G23 for various fiber volume fractions, in GPa.

Click on image to start an interactive JWave session

Figure 14. Polar diagram of graphite/epoxy Poisson's modulus ν32 for various fiber volume fractions, in GPa.
NOTE: yellow, green, and purple line indicates negative Poisson's ratio.

NOTE: since all fibers are aligned along the z axis (3-direction), the x-y (1-2) plane
             is isotropic and hence we have polar symmetry about the z axis and 3-D
             images of these elastic properties would also be symmetric about the z axis.
             We call this condition "transversely isotropic" which is similar to hexagonal
             symmetry. Because of polar symmetry, simple 2-D figures are sufficient to
             show this anisotropy.

A closer look at these figures reveals some obvious facts, but also some unexpected new information. If all fibers are aligned along the z axis (3-direction) then we would expect that the elastic stiffness would be largest along the z axis (3-direction) and less stiff in the x-y (1-2) plane. Indeed this trend is confirmed geometrically. However unexpected properties are observed such as a large negative (-0.354) Poisson's ratio, see purple line in Figure 14. A negative Poisson's ratio would require that the material would expand both in the load direction and transverse direction to the load direction. If the material is transversely isotropic the transverse expansion would be the same in all directions within the transverse plane. Putting this all together would require that such a material when loaded in the fiber direction would expand in all directions, hence it's volume would increase under a uniaxial tension load.

Perhaps the reader might find other interesting mechanical behavior by experimenting with either the links to "Interactive JWave session" shown above each of the figures (12,13,14), or below an alternate link is provided that provides the student access to a more comprehensive interactive session with access to: (1) computer programs, (2) sample data files, (3) numerical results, and (4) images.

The elastic properties: Young's Modulus, E, shear modulus, G, and Poisson's Ratio, ν, in each orthogonal plane can be used to classify the nine independent orthorhombic elastic constants in terms of engineering properties.


Figure 15. Transversely-isotropic carbon fibers and their composites.

With the mechanical behavior ("deformation") shown in each plane of Fig. 15 it is easier to see how each of these nine properties are independent. If all the fibers are aligned along the 3-axis, it is easy to see the 12-plane is transversely isotropic and the number of independent elastic properties reduce to 5 (hexagonal symmetry) and the isotropic relationships can be used in the 12-plane. The large negative Poisson's ratios predicted in Fig. 14 is most likely due to the nano-structure of the graphite fibers where basal planes become aligned parallel to the circumference when pyrolyzed. This pyrolyzation also creates a lamellar structure of basal planes aligned parallel to the fiber longitudinal axis ("L") and explains the unusually high elastic stiffness along the x3 axis. Together these two nano-structural observations explain the fiber has transversely- isotropic ("hexagonal") symmetry. Here we choose the fibers longitudinal axis to be parallel to the x3 axis because for hexagonal symmetry the x3 axis is the conventional "unique" symmetry axis. However in engineering design the x1 axis is choosen as the fiber-axis convention.

A brief literature review has revealed some interested hypothetical structures for pyrolyzed carbon fiber structures, see Figure 16. Although preliminary observations can explain how this nano-structure results in a large modulus along longitudinal fiber axis, much more research needs to be done before we can conclude any nano-structure property relationship that can explain the large negative Poisson's ratios.


Figure 16. Carbon fiber pyrolyzed structures Ref.[11] and Ref.[12].

Because so little is known about these transversely isotropic ("hexagonal") fiber properties, engineering designers make some gross approximations. For example the five independent elastic properties shown in Fig. 15 are further reduced to four independent elastic properties by equating all three shear moduli and similarly all three Poisson's ratios are also equated.

E3,       E2 = E1,
G23 = G13 = G12,
ν32 = ν31 = ν21

These same fiber-reinforced materials and their peculiar anisotropies will be used in the next section on "Laminated Plate Theory". As a result we can expect some unusual mechanical behavior of laminated plates fabricated with these materials at the micro level.

Geometric representations of engineering elastic properties, E, G, ν, although interesting, can only reveal variations of just that one property. The modulus E3 = S33-1 represents only one of nine components of a fourth order tensor and consequently does not represent the full anisotropy of all components of the stiffness tensor. Next we will show how to geometrically represent anisotropy that includes all possible terms (components) of the fourth order stiffness tensor.

Geometric Representations of Elastic Anisotropy: All Tensor Components, Cijkl

To include all components of a fourth order stiffness tensor in a geometric representation of anisotropy requires we look closer at the equations of motion for an anisotropic solid (dynamic mechanical behavior). Here our objective is to look at a spherical disturbance such as a dilitational pulse that expands equally in all directions: invoking Huygen's principal the reader can envision a very small plane wave exists on the surface of a very small sphere in the center of an anisotropic crystal. One can simply imagine that from Huygen's principal each of these plane waves travels in a specific direction, direction cosines, νi, at a speed that corresponds to elastic properties unique to the same direction. Hence, plane waves traveling in different directions will travel at different speeds and the continuous collection of all of these plane waves, although initially a sphere, soon deviates into a nonspherical shape simply because plane waves will travel faster in stiffer directions and slower in less stiff directions. This idea was modeled and simulated for a highly anisotropic Calcium-Formate crystal, "Spherical Wave Simulation in 2-D Anisotropic Media". The three-dimensional wave surface for a Calcium-Formate crystal was first plotted in Ref.[13]. We will show that this wave surface topology does indeed uniquely describe the complete set of components of the fourth order stiffness tensor. These resulting geometries are now being used as new sub-classification schemes within orthorhombic symmetry, Ref.[14].

First we start with the equations of motion for a continuum.

(35)

Recall equation (17)

(17)

and substitute the strain-displacement relationship, equation (7),

(7)

into equation (17) yields

(36)

Recall that the strain tensor is symmetric which further reduces equation (36).

(37)

Substitituting equation (37) into equation (35) yields the equation of motion in terms of displacements.

(38)

This equation can be further reduced if we assume the material is homogeneous.

(39)

With this simplification equation (38) reduces to a form where we can now assume a solution for displacement and substitute it into this equation.

(40)

Let us assume a small plane wave periodic disturbance, see Figure 17.


Figure 17. Mathematical model for a plane wave disturbance: oscillating plane written in exponential format in terms of wave velocity, v, propagation direction, νl, and particle vibration direction, αk

The periodic displacement shown in Figure 17 uses an exponential functional form instead of sines and cosines because it will be easier to differentiate. The reader can envision a plane wave in Figure 17 where the wave number coefficients, ki, are terms used to define a plane (this will perhaps require a review of basic descriptive geometry in freshman calculus) and the circular frequency, ω, when multiplied by time, t, replaces the constant for a plane with a term which models the oscillatory motion of that plane.

(41)

The terms in equation (41) can be replaced with wave velocity, v, and direction cosines of the propagation direction, νi which will allow the reader to interpret the results of this particular eigen-problem in a more meaningful way.

(42)

Substituting equation (42) into equation (40) and after much manipulation, equation (40) reduces to an eigen-value problem which is written below in both tensor and matrix notation.

tensor equation: (43)

matrix equation:

(43)

If we use the matrix form of equation (43) it is easier to see that the velocity terms along the diagonal, ρv2, are the eigen-values and the displacement direction cosines, αi are the eigen-vectors. Closer examination of the solution to this particular eigen-value problem reveals that both the eigen-values and eigen-vectors can only be functions of ALL stiffness tensor components for an assigned propagation direction, νi. Hence we conclude that if we solve for eigen-values (wave speeds) in all possible propagation directions, νi, and connect points to form a continuous surface, this would generate a 3-D shape of wave speed for each eigenvalue. Since there are three eigen-values, equation (43) predicts three wave surfaces: 1) longitudinal which is usually the fastest and observed as the largest shapes, and 2) two different shear ("transverse") waves which are usually slower and consequently observed as smaller shapes. The eigen-vectors which are particle vibration direction cosines (0 to 90 degrees) can be mapped onto the eigen-value surfaces as color. Hence the color would reveal the type of wave: longitudinal (0-degrees i.e. purple) and transverse (90-degrees i.e. red) where green would indicate a mode transition from longitudinal to transverse. With colors the observer can quickly determine the wave type and discover possible mode transitions. Graphical definitions of these and other properties associated with this solution are given by the problem statement / updated statement.

Some obvious geometric features that correspond to anisotropy are: 1) stiffer directions map as a larger distance from the geometric center which would be seen as "bulges" and 2) less stiff directions map as smaller distances from the geometric center which would be seen as "depressions". Highly anisotropic materials would appear as significantly deviating from the isotropic spherical geometry. Graphite-Epoxy and Calcium-Formate exemplifies all of these features, see Figure 18 and Figure 19.


Figure 18. Wave surfaces predicted for Graphite-Epoxy, Hexagonal symmetry
NOTE: light-blue/green indicates a mode transition -- a very rare anisotropy which
is only observed in Calcium-Formate, spruce-wood, and graphite/epoxy.

INTERACTIVE COMPUTER PROGRAM: Experiment with Hexagonal Anisotropy / update results



Figure 19. Wave surfaces predicted for Calcium-Formate, orthorhombic symmetry
NOTE: light-blue/green indicates a mode transition -- a very rare anisotropy which
is only observed in Calcium-Formate, spruce-wood, and graphite/epoxy.

INTERACTIVE COMPUTER PROGRAM: Experiment with Orthorhombic Anisotropy / update results

Figure 19 shows the solution of equation (43) for an orthorhombic fourth order stiffness tensor of Calcium-Formate, Ref.[13]. This particular tensor has a rare anisotropy that requires all three wave surfaces to connect into a single contiguous surface, see "combined" wave surface of Figure 19. A simple movie demonstrates both the connected surfaces and mode transitions previously described. The only other known material that exhibits this anisotropy (geometry) is Spruce wood. References [13-15] identify four possible geometries not previously known. With these geometries it is now possible to derive a subclassification scheme for orthorhombic anisotropy.

By studying the wave surface geometries it is also possible to study and compare elastic property anisotropy from isotropic to orthorhombic symmetry. Below is a table listing elastic material properties for each of these symmetries.

Table 3. Representative symmetries: isotropic, cubic, hexagonal, tetragonal, orthorhombic

Material: Stiffnesses GPa Density
gm/cm3
C11 C22 C33 C44 C55 C66 C12 C13 C23
Stainless Steel (isotropic) / Update 7.88 261.0 261.0 261.0 77.5 77.5 77.5 106.0 106.0 106.0
Aluminum, (cubic) / Update 2.7 108.0 108.0 108.0 28.0 28.0 28.0 62.0 62.0 62.0
Zinc (hexagonal) / Update 7.10 143.0 143.0 50.0 40.0 40.0 63.0 17.0 33.0 33.0
Aligned Graphite/Epoxy (hexagonal) / Update 1.61 14.5 14.5 160.0 7.07 7.07 3.62 7.26 5.53 5.53
Tellurium Oxide (tetragonal) / Update 4.99 55.7 55.7 106.0 26.5 26.5 65.9 51.2 21.8 21.8
Benzene (orthorhomibic Type-1) / Update
( C11, C22, C33 ) > ( C44, C55, C66 )
1.06 6.14 6.56 5.83 1.97 3.78 1.53 3.52 4.01 3.90
Calcium Formate (orthorhombic Type-2) / Update
( C11> C33> C 66> C22> C55> C 44 )
2.02 49.2 24.4 35.4 10.5 12.2 28.2 24.8 24.5 14.5
Hypothetical (orthorhombic Type-3) / Update
( C11> C66> C55> C22> C44> C33 )
6.00 180.0 45.2 10.3 20.4 45.5 90.6 0.100 0.200 0.300
Hypothetical (orthorhombic Type-4) / Update
( C11> C22> C66> C44> C55> C33 )
3.00 80.1 60.2 10.3 30.4 20.5 40.6 0.100 0.200 0.300

Figures below correspond to the last four orthorhombic symmetries and their corresponding elastic property inequality.


Figure 20. Examples of the four basic types of geometries and their corresponding
inequalities that can be used to subclassify within orthorhombic symmetry

The inequalities for each of the four types of orthorhombic symmetries are the necessary and sufficient condition that defines each unique geometry. These geometries and their corresponding inequalities can be used to create a subclassification scheme within orthorhombic symmetry. For those interested we have provided a graphical proof.

Other highly anisotropic materials such as graphite/epoxy can also significantly influence dynamic mechanical behavior. The section on Simulation and visualization of wave propagation in unidirectional graphite/epoxy composites Ref.s [16], [17], and [18], summarizes how unidirectional graphite/epoxy composites like Calcium Formate exhibit mode transitions but without resulting in a connected wave surface. Unlike Calcium Formate, which is a crystal, graphite/epoxy derives its anisotropy from its' constituent fiber and matrix properties (microscale), hence elastic properties at the microscale can signficantly influence anisotropy and dynamic mechanical behavior at the macroscale. For graphite/epoxy it can be shown that not only can anisotropy have a significant effect on dynamic mechanical behavior, see Figure 21, but degradation of either the fiber or matrix phase can preferentially alter the direction of the wave propagation from the initial undegraded state, see Figure 22. A section on simulation-visualization in unidirectional graphite/epoxy composites allows the reader to explore the influence of anisotropy on these and other possible dynamic mechanical behaviors using a small grid (30x60) or large grid (45x180) finite element grid model. For confirmation with an exact solution a simple isotropic one-dimensional model is provided here for the reader. Experimental confirmation of the accuracy of these physical based simulation models is given in Ref.[10] Figure 14, see Figure 23.


Figure 21. Schematic of simulation and animation summarizing
       wave properties of a full field simulation.


Figure 22. Influence of preferential microscale elastic component
          degragation on macroscale mechanical behavior.


Figure 23. Experimental confirmation of beam steering, Ref.[10] Figure 14.

The discussion above provides an excellent example of how elastic properties scale from the nanoscale (fiber transversely-isotropic nano-structure) to the microscale (transversley-isotropic unidirectional composite) and ultimately influence the dynamic mechanical behavior at the macroscale (beam steering and mode transitions). The diagram below illustrates how the wave length of a stress wave scales from the macro to the nano structures. Indeed the commponents of the elastic stiffness tensor are derived from the nanostructure of the atomic spacing.


Figure 24. Wave length and elastic properties scale from nano to macro.

Other sections on: 1) Visualization of Stress Waves in Elastic Solids and 2) Simulation and Visualization of Energy Waves Propagating in 2-D and 3-D Anisotropic Media, provide more information on how anisotropy controls dynamic mechanical behavior.

We have demonstrated that anisotropy in crystals like Calcium-Formate and undirectional composites like graphite/epoxy exhibit unusual dynamic and static mechanical behavior. As designers we can take advantage of this anisotropy. Because the elastic properties of graphite/epoxy can be tailored, engineers can design the material for the application instead of designing an application around fixed material properties. One popular method of exploiting the anisotropy of graphite/epoxy, in the design of a desired mechanical behavior, is to use laminate plate theory.



4.2 Introduction to Mechanical Behavior of Laminates (Laminate Plate Theory): Analytic Model

With laminated plate theory we bridge the scale of modeling mechanical behavior from the micro "lamina" scale to the macro "laminate" scale. Figure 25 shows the construction of a laminated plate with individual lamina layers. In practice laminated plates consist of hundreds stacked lamina but here we study simple laminated structures.


Figure 25. Lamina-Laminate Definitions

In the previous section, 4.1 Introduction to Mechanical Behavior of Anisotropic Materials, mechanical behavior of individual anisotropic layers, "lamina", of the "laminated" plate can now be used to model the mechanical behavior of the laminate. The purpose of this section is to introduce the reader to the basic ideas of laminated plate theory so that the reader gains an understanding of the relationship between the two scales. To gain a working knowledge of laminated plate theory, we encourage students to take ESM4044 MECHANICS OF COMPOSITE MATERIALS and ESM6104 MECHANICS OF COMPOSITE STRENGTH AND LIFE . With a working knowledge of laminate plate theory, designers can taylor the elastic properties and orientation of each lamina, hence optimize their design by controlling the properties and orientations of each lamina. At the end of this section we will outline several examples on how the mechanical response of a laminated plate can be controlled ("designed") by lamina properties, but first we briefly outline the derivation of the analytic model needed to understand how properties at the microscale control mechanical behavior at the macroscale.

Analytic Model:

Consider an orthorhombic lamina in plane stress and let the 1-2 axes shown in Figure 26 be aligned with the fiber axis. In plane stress we assume the lamina is very thin and the stress through the lamina thickness is very small.


Figure 26. Lamina axes defined

(44)

The elements of [Q] are related to the engineering constants E,G, and ν as follows

(45)

(46)

Transformation of lamina stress and strain from the 1-2 axes to the rotated x-y axes is shown here in matrix notation. This transformation is different from the transformation for all components of the second order tensor because here we have only three components of stress and strain because of the plane-stress condition.

(47)

where

(48)

The lamina stress-strain relationship in the rotated x-y coordinates

(49)

In reduced matrix notation

(50)

where

(51)

The inverse stress-stain relationship relates stresses(known) to strains(unknown) with the compliance, [S], matrix

(52)

where

(53)

Hence there is an inverse relationship between the [S] and [Q] matrices.

(54)

The equations above apply to each individual laminate layer ("lamina"). Next we stack each lamina into a laminate where we assume each lamina, denoted with a subscript k, is bonded ("glued") to the adjacent lamina. The composite structure is referred to as a laminated plate, see Figure 27.


Figure 27. Lamina-Laminate Configuration Geometry Defined

The laminated plate shown in Figure 28 is loaded, in the Z-direction. Deflections are assumed small and the Kirchoff thin plate approximations can be used here to derive relationships between in-plane and out-of-plane displacements, curvatures, twists, and strains.


Figure 28. Thin Plate Deflections

From Figure 28 it is not difficult to see how in-plane uo and vo displacements are augmented by out-of-plane (Z-direction) displacements and rotations.

(55)

From equation (7) the strain-displacement relationship are rewritten here for components of displacement, u and v, confined to the x-y plane ("in-plane").

(56)

From the expression of radius of curvature, ρ, we derive expressions for curvature, kx and ky assuming small angular deflections.

(57)

A similar expression for twist, kxy, can be derived from Figure 29 where twist is defined as an "angular change" per "unit length".


Figure 29. Exaggerated Thin Plate Deflections and Twists Defined

The coefficient of 2 is introduced by convention for engineering shear strain.

(58)

Combining equations (57) and (58), into (56) yields

(59)

Substituting equation (59) into equation (49) yields a more general stress-strain relationship for the k-th lamina which includes out-of-plane displacement implicitly in the curvature terms.

(60)

With stresses defined in terms of the Z-coordinate, through the laminate thickness, we integrate these expressions to get the laminate resultant in-plane loads (Nx, Ny, Nxy), and moments (Mx, My, Mxy).


Figure 30. Laminate Plate In-Plane Loads and Moments Defined

(61)

Substitute equation (60) into equation (61) and simplify.

(62)

Equation (62) can be simplified

(63)

where

(64)

Figure 31 shows the relative position of each k-th lamina used by equation (61) to integrate through the thickness to calculate the laminate loads, (Nx, Ny, Nxy), and moments (Mx, My, Mxy).


Figure 31. Lamina Positions Defined in the Laminate for Integration of Resultant Loads

Substituting the [A], [B], and [D] matrices defined in equation (64) into equation (63) results in the final relationship that defines the mechanical behavior of laminated fiber-reinforced plate. Equation (63) shows the relationship between known in-plane strains ( εox, εoy, γoxy ) and curvatures ( kx, ky, kxy ) and unknown in-plane loads and moments.

(65)

And the inverse relates known in-plane loads to unknown in-plane strains and curvatures.

(66)

The [A], [B], and [D] matrices for the laminate are similar to the [Q] matrix for the individual lamina. If you followed the basic ideas in this derivation it is not difficult to understand how microscale lamina properties control the macroscale laminate properties. In deed we now have a predictable way to control, and hence design, the mechanical behavior of a laminate at the macroscale.

Like the lamina at the microscale, anisotropy of the laminate has a significant effect on the mechanical behavior of the laminate. This can best be observed by inspection of equations (65) and (66). If the individual lamina are stacked and glued together such that the [B] matrix is zero then the mechanical behavior exhibits a decoupling of moments and in-plane loads. For example curvatures only affect moments and inplane strains only affect inplane load. On the other hand it is not difficult to see that if [B] is not zero then there is coupling between curvatures and inplane loads, and coupling also exists between inplane strains and moments. This is rather peculiar. For example if we apply an inplane strain to a specially configured laminate then the laminate can bend more then it would extend ("stretch"). Similarly if we apply an inplane load the laminate will twist. The later is actually used to an advantage by designers of aircraft structures where aircraft wings, that are constructed by fiber-reinforced laminates, are designed to twist into a predetermined aerodynamic angle-of-attack as the wing bends and extends. This design methodology is called PASSIVE aero-elastic tailoring. With "smart"-"intelligent" structures it is possible to insert active piezoelectric layers into the laminate stacking sequence and ACTIVELY control the mechanical behavior of the laminate. The degree of coupling is controlled by the [A], [B], and [D] matrices.

In this section a few simple examples are presented that illustrate how simple laminate configurations can provide predictable macroscopic mechanical behavior by inspection of equations (65) and (66).

Case I: Symmetric Laminates, [B]=0

By inspection when [B] is zero, inplane loads and moments decouple in equations (65) and (66).
Inplane loads {N} only produce inplane strains {ε} and inplane moments {M} only produce
curvature {k}

Case II: Specially Orthotropic Symmetric laminates, A16=0 and A26=0

Bi-directional laminates consist only of 0/90 degree lamina. It is not difficult to show Q16 and Q26
are zero. Hence A16 and A26 are also zero and consequently the inplane shear strain {γ} decouples
from the inplane normal strains. Therefore shearing loads can only produce shearing strains and
inplane normal loads can only produce inplane strains

Case III: Specially Orthotropic Bending Behavior, D16 & D26 -> 0 , no. of lamina -> large

For a large number of symmetric pairs, +θ at +h and -θ at -h, D16 and D26 go to zero. Hence
moments {Mx} & {My} and torque {Mxy} decouple so that bending caused by moments do not
introduce twisting and torque that causes twisting does not introduce bending.

Accept for a few simple cases outlined above, equations (65) and (66) are sufficiently complex where the designer can not deduce a solution by observation. These equations must be solved numerically, hence the designer must make educated guesses and experiment by submitting multiple jobs and hopefully converge on a particular design. To control the mechanical behavior of the laminate requires attention on controlling the [A], [B], and [C] matrices. For complex stacking sequences this becomes difficult. Many references on designing laminated structures exist. As an introduction we suggest the reader look at Refs.[19,20]. To that end we provide an interactive tool below that will allow the reader to experiment and design the desired laminate mechanical behavior.

INTERACTIVE COMPUTER PROGRAM: Experiment with Laminate Plate Analysis Program

In summary we have shown the significant role anisotropy has on the mechanical behavior of laminates both at the microscale and macroscale. The bridge between the microscale ("lamina") and the macroscale ("laminate") is accomplished with analytic models of laminated plate theory. In the theory of laminated plates we have ignored the stresses and strains in the z-direction (thin-plate assumption). For thicker laminates at the macroscale these stresses can not be ignored indeed these stresses are significant and cause the formation of ply-cracks, see Figures 1,3, and 4. In the next section we provide another model that predict stress singularities at the laminate free-edge. Singularities at the laminated free-edge are also largely influenced by anisotropy at the lamina microscale.


5.0 Laminate Singularities Caused by Anisotropy: "Free-Edge Problem": FEM Model

Mechanical behavior of thick laminated plates can be accurately predicted by simple laminated plate theory, but the stresses within these thick laminates are not accurately predicted by classical laminated plate theory. For thick laminates stress singularities result in controlling the sequence of failure events that lead to the ultimate fracture of the laminate. Figure 32 shows a typical sequence of damage events that lead to ultimate fracture at the macroscale.


Figure 32. Typical sequence of damage events leading to laminate fracture

To accurately predict stress singularities in thick laminates we introduce the laminate "Free-Edge" problem. In this section we provide:

Below is a description of the laminate geometry and finite element grid representation used to define the stress "Free Edge" problem. The figure below shows the laminate stacking sequence and FEM grid for a specific problem in Ref.[21].


Figure 33. Laminate geometry and FEM grid defined for "Free-Edge" problem, Ref.[21].

NOTE: the enlarged area showing the smallest element of the F.E.Model at the free edge is reaching the conceptual limit of the assumed homogeneous approximation ("continuum") for material properties in each element at the free edge. The finite element mesh superposed over a SEM photograph shows the heterogeneity of each graphite fiber diameter, 5 microns, is observed as a fuzzy white region. Here we are pushing the limit and assuming each element still behaves like a continuum with three fiber diameters in the z-direction and ten fiber diameters in the horizontal y-direction. For FEM elements smaller than this a continuum model can no longer be used. Instead a different model needs to be created at Y/B > 0.998 where mechanical behavior is now controlled by the individual fiber and matrix components. Hence the continuum approximation breaks down at this point.

Example of results prior to First-Ply-Failure (FPF) in the WET ("water saturated") condition, shown below, demonstrate the stress singularity of the stress distribution of the normal stress in the z-direction along each laminate interface. This singularity can only be accurately modeled up to Y/B=0.998 as a continuum.


Figure 34. Sample results of laminate edge-stress singularities, Ref.[21].


Figure 35. Sample results of stresses at laminate stress free edge, Ref.[21].

You might want to reproduce the results shown in this Figures 34 and 35 by going to the interactive computer program. You will also need to use the Laminated Plate Analysis program which is used to convert the Nx=1900 lb./in. to a strain in the x-direction. The strain in the x-direction, calculated by LELPA, is then used as the out-of-plane strain normal to the H-Z plane in the form below. Notice that the σz stress reverses and becomes negative ("compressive") along the +45/-45 interface and 0/+45 interface most likely because of the combination of residual stresses caused by the thermal contraction at low temperatures, ΔT=-180F, and swelling caused by moisture absorption, ΔM=1.2% moisture.

All stresses at the laminate stress free edge and interior are summarized below in Table 4. Elastic properties used to calculate these stresses are listed in Table 5.

Table 4. Edge and Interior Stress for Type-I Gr/Ep Laminate, Ref.[21]

Table 5. Elastic Properties for Fiber, Matrix, and Unidirectional Composite, Ref.[21]

In section 5.0 above a FEM model predicted stress free edge singularities in type-I quasi-isotropic laminate, Ref.[21]. In section 6.0 below a FEM model predicts stress singularities near a 90o ply crack terminating at the 0/90 interface in a type-II quasi-isotropic laminate, Ref.[1].


6.0 Laminate Singularities Caused by Ply Cracks: FEM Model

After the formation of a ply crack, see Figure 25, new singularities occur near the ply crack tip. The growth of this crack can continue into the adjacent ply, see Figures 1, 3, 4, and 5, or debond along the interface. Below we provide models that demonstrate how anisotropy and residual stress can augment singularities near ply cracks and influence how ply cracks grow and control laminate fracture strength. In this section we provide:

Stress Analysis Summary: Details of the stress analysis is given in Reference [1]. Here we show a schematic of the laminate and finite element mesh that was used to model a 90o ply crack in a [0/90/+45/-45]s laminate and the stress gradients ("concentrations") in the load bearing 0o ply near the ply crack tip, see also Figures 3 and 4.


Figure 36. Schematic and FEM Grid of Ply Crack in a [0/90/+45/-45]s laminate, Ref.[1]


Figure 37. Mode I stresses in the 0 degree load bearing ply above the ply crack, Ref.[1]

Summary: With the interactive computer program the reader can vary properties at the microscale and observe the effect these changes have on the stress distributions near the ply crack for two different types of quasi-isotropic laminates. Studying different laminates would require creating different FEM meshes. Future examples will target interactive FEM mesh generation. With this interactive computer program the reader can reproduce these results and study how variations on elastic properties and other properties at the microscale can influence the macroscopic stress state and influence laminate fracture.

Similar to the stress "Free-Edge" FEM problem in section 5.0, the FEM model above predicts stress singularities near the crack tip of a 90o ply crack in a eight-layered type-II quasi-isotropic laminate. Below we introduce a FEM model with a 90o ply crack in a [0/90]s laminate with and without a woven macro-structure, Kriz Ref.[23]. These FEM models are complete where stresses plotted in Ref.[23] are reproduced here, but graphical plots of these stress gradients are not shown here because the interactive NPIB modules are under construction.

---- UNDER CONSTRUCTION ----


7.0 Cracks Near Interfaces Between Dissimilar Anisotropic Materials

In previous sections we studied stresses near laminate stress free edges and ply cracks that exist at an interface between two dissimilar anisotropic materials. These results were predicted by a Finite Element Model (FEM) that did not include singularities. Instead stresses gradients near the free edge and ply crack tips were approximated with a high density FEM mesh near the singurality. In this section we introduce an analytic model that predicts these singularites: 7.1) at a bimaterial interface intersecting a stress free edge, see Figure 38, and, 7.2) at a ply crack tip intersecting an interface seperating two dissimilar anisotropic media, see Figure 40. Both models are solved using Stroh's method, Ref.[24]. These singularites can be introduced into enriched FEMs that included a variety of different boundary conditions (BCs) without using high density FEM meshes. With these enriched FEMs, stress intensity factors can be solved for as unknowns along with the nodal displacements for each unique geometry and BCs. Section 8 will outline how to create a finite element model with elements "enriched" with these singularities, Benzley Ref.[29] and Heyliger & Kriz Ref.[30].

For both models, shown in Figs. 38 and 40, it is important to note that the singularity coincides with the origin of the coordinates. Hence each model has a unique coordinate system. Since the objective here is to reproduce published results of these two models the same coordinate system established by Ting et al. will be used. This choice will require rederivation of an arbitray rotation in the 1-2 plane for orthorhombic symmetry outlined in section 4.1. The derivations below will follow the notation and procedures outlined by Ting et al, Refs. [25] and [26]. The notes below will expand on the solution introduced by Ting et al. by developing the numerical approach to solve for these singularities. These singularities were also calculated by using Green's functions at dissimilar anisotropic material interfaces terminating at a stress free edge, Tewary & Kriz Refs. [27] and [28].

7.1 Stress free edge singularities:

Figure 38. Coordinates defined for stress free edge singularity Ref.[25].

T.C.T. Ting and S.C. Chou, Ref.[25], used Stroh's method to solve for the singularity where the interface seperates two dissimilar anisotropic regions and terminates at a stress free edge, see Figure 38. Here we only expand on the development of the numerical solution and use the same notation in Ref.[25] except the exponent on the radius, r, is changed from, δ, to -k similar to the exponent in Ref.[26].

The solution is organized into three parts: 1) transformation of the fourth order stiffness tensor, Cijkl, for arbitrary rotation in the 1-3 plane, see eqn 35 Ref.[25], 2) calcuate eigenvalues, pk, and eigen vectors, νi and τij, from equations (9) and (10) Ref.[25], and 3) calcuate exponent (singularity) on radius, r, from the determinate of equations that are created when equations (30) and (31) of Ref.[25] are combined with appropriate boundary conditions.

7.1.1 Tensor Transformation Cijkl in 1-3 plane:

Hand-written notes are provided for the transformation of the fourth order stiffness tensor for arbitray rotations in the 1-3 plane.

7.1.2 Calculate Eigenvalue/Eigenvectors, pk / νi and τij:

We start with equations (9) and (10) in Ref.[25].

(9) Ting
(10) Ting
where
(11) Ting
(12) Ting

Below we expanded and simplify eqn. (10) into an eigenvalue problem that we can solve by using IMSL scientific subroutines. Combining equations (11) and (12) from Ref.[25] along with the rotated transformed stiffnesses, yields

(67)
where
(68)

Substituting eqn (68) into eqn (10) of Ref.[25] yields

(69)

Note that PL exists in every term of the coefficient matrix in eqn. (69) who's determinate yields a cubic characteristic equation in PL2. Although it is possible solve for the roots of this characteristic equation, numerical solutions to this type of problem yield more accurate eigenvalues. Hence eqn.(69) is further developed into an eigenvalue problem where PL are complex eigenvalues and their corresponding complex eigenvectors, ν i,L, where the ,L subscript associates the eigenvector with the appropriate eigenvalue (L=1,2,3). The L is not a free indice and the comma does not imply differentiation

The eigenvalue format is

(70)

It is possible to expand equation (69) into the form of equation (70). First we rewrite equation (69).

(71)

where we can write equation (71) in abbreviated notation.

(72)

Multiply equation (72) by the inverse of the coefficient matrix, [A'], yeilds.

(73)

where [A']-1[A'] becomes the indentity matrix [I], and matrices [B'] and [C'] are premultiplied by [A']-1 and become [B] and [C] respectively.

(74)

The array [A] in eqn.(70) can be constructed such that

(75)

where the nonzero components in array [A] of eqn.(70) are summarized below.

(76)

To show eqn.(70) and eqn.(73) are equivalent requires the following expansion. First substitute eqn.(75) into eqn.(70) where the {Z} eigenvectors are partioned into Xi and Yi components.

(77)

where

(78)

Equation (77) can be rewritten interms of Xi and Yi components.

(79)

Expanding eqn.(79) yields two equations

(80)

These two equations can be combined into one equation by substituting for Xi, hence eliminating Xi, and rearranging terms.

(81)

where comparing equations (81) and (73) yield the necessary eigenvalues and eigenvectors

(82)

Comparing equations (81) and (73) verifies that equation (70) can be used to solve for PL as eigenvalues instead of roots to the characteristic equation derived from the determinate of the coefficient matrix in equation (69). Both eigenvalues, PL, and eigenvectors from equation (81) will be used in the next section to calculate the singular exponent. IMSL routines EVCRG are used to solve equation (70) and EPIRG is used to obtain a performance index for the solutions.

7.1.3 Calculate Singular Exponent, -k:

Using Stroh's method Ting et al., Ref.[25], obtained the functional form for displacements and stresses.

(30) Ting

(31) Ting

where again we note that ,L does not imply differentiation and L is not a free indice but is used as a subscript to identify eigenvectors associated with the L-th eigenvalue. Here L is assigned 1,2,3 corresponding to the three eigenvalues of eqn.(69).

(26) Ting

The equations above have been slightly modified. The exponent has been changed from δ to -k and the coefficients, ML and NL have been renamed aL and a'L.

To solve for the exponent singularity, k, we create a systems of 12 linear equations from Ting's equations (30) and (31) and the boundary conditions listed below. Coefficients of these equations are written in terms of k and solved for k as a root to the determinate of the coefficient matrix of these equations. Since the solution is to find only one of the roots between 0.0 and 1.0, this root is first approximated by calculating the determinate as a function ("function-subroutine") of k where a value for k is incremented from 0.0 to 1.0 until the determinate changes sign (+/-). Within this interval a root exists. With an approximation for the root interval other IMSL routines can be used to obtain a more accurate solution for k. Next we construct the twelve equations for the boundary conditions defined in Fig. 38.

The boundary conditions are:

(83)

From these boundary conditions 12 equations are generated whose coefficients are functions of elastic properties and the exponent on the radius vector whose origin is located at the crack tip where the singularity exists. Below we generate this coefficient matrix.

Boundary condition 1:   

(84)

Boundary condition 2:   

(85)

Boundary condition 3:   

(86)

Boundary condition 4:   

(87)

The blocks defined above are substituted into the schematic shown below.


Figure 39. Block schematic of system of equations for determining singularity
where bimaterial interface terminates at a stress free laminate edge.

IMSL routine ZBREN is used to find a solution (root) to the determinate of the coefficient matrix shown in Fig. 39. The ZBREN requires an initial guess where the determinate changes sign in the interval from 0 to 1.0. Hence a function subroutine, F, is setup which calculates the determinate using IMSL routines LFTRG and LFDRG.

Equations summarized above are compiled into a Fortran computer program (edgesing.f) which is provided below as an interactive module (Module15a). This same program can be accessed and downloaded from Module15a listed in the Microstructure Interactive Computer Program Modules Table section of the CRCD pages.

With the interactive computer program provided above the student can reproduce the singularities tabulated below. In particular the default values set in the NPIB form of the computer program is configured to calculate the singularity of 2.9388 x 10-2 for θ'=0 and θ=75 with the Pipes-Pagano approximation for elastic properties:

E1=E2=+2.100E+06,    E1=+2.000E+07
G12=G13=G23=+0.850E+06
ν21= ν31= ν32=+0.210E+00

Table 6. Stress singularity at the free edge interface of an angle-ply graphite/epoxy laminate, Ref.[25]
Exponent
of r-k
θ'=0 θ'=90 θ'=-θ
θ=0 - 3.3388 x 10-2 -
θ=15 1.3528 x 10-4 3.2814 x 10-2 6.4322 x 10-4
θ=30 2.6286 x 10-3 2.8682 x 10-2 1.1658 x 10-2
θ=45 9.6461 x 10-3 2.0575 x 10-2 2.5575 x 10-2
θ=60 1.9866 x 10-2 1.0519 x 10-2 2.3346 x 10-2
θ=75 2.9388 x 10-2 2.6785 x 10-3 8.9444 x 10-3
θ=90 3.3388 x 10-2 - -


7.2 Ply Crack singularities:

Figure 40. Coordinates defined for ply crack singularity Ref.[26].

T.C.T. Ting and Hoang, Ref.[26], used Stroh's method to solve for the singularity at a ply crack tip acting normal to an interface seperating two dissimilar anisotropic regions, see Figure 40.

It is important to note that singularities coincide with the coordinate system origin. Hence each interface singularity has a unique coordinate system. Since the objective here is to reproduce published results the same coordinate system established by Ting et al. will be used. This choice will require rederivation of an arbitray rotation in the 1-2 plane for orthorhombic symmetry outlined in section 4.1. The derivations below will follow the notation and procedures outlined by Ting et al, Ref.[26]. For this module we have expanded on the solution introduced by Ting et al. and developed the numerical approach to solve for these singularities.

The solution is organized into three parts: 1) transformation of fourth order stiffness tensor, Cijkl, for arbitrary rotation in the 1-3 plane, see eqn 35 Ref.[25], 2) calcuate eigenvalues, pk, and eigen vectors, νi and τij, from equations (9) and (10) Ref.[25], and 3) calcuate exponent, k (singularity) on radius, r, from determinate of equations that are created when equations (21) and (22) of Ref.[26] are combined with the appropriate boundary conditions equation (36) of Ref.[26]. NOTE: The geometry shown in Fig.40 results in a set of boundary conditions different then those used in the previous section 7.1.

(88)

From these boundary conditions 18 equations are generated. Coefficients are functions of elastic properties and the exponent on the radius vector whose origin is located at the crack tip where the singularity exists. Unlike boundary conditions for the free edge problem of section 7.1, there is an additional notation where + or - appear as superscripts on variables that indicates these variables exist in regions x3 > 0 and x3 < 0 respectively.

Partial results from Ting et al. Ref.[26] are listed below in Table 7 with the Pipes-Pagano approximation for elastic properties:

E1=E2=+2.100E+06,    E1=+2.000E+07
G12=G13=G23=+0.850E+06
ν21= ν31= ν32=+0.210E+00

Table 7. Stress singularity at an interface with a ply crack acting normal to the interface, Ref.[26]
Exponent
of r-k
θ=0 θ=15 θ=30 θ=45 θ=60 θ=75 θ=90
θ'=90 .676635
.500212
.5
.676012
.502156
.500222
.672441
.505068
.500236
.660488
.504813
.500227
.630471
.502297
.500314
.567537
.500314
.500048
.5
θ'=75 .675051
.500189
.433918
.665411
.500305
.446234
.646234
.500401
.466560
.616354
.500401
.485037
.571159
.500269
.496300
.5 .499859
.499688
.4322561
θ'=60 .665263
.499839
.379767
.642118
.500352
.404849
.605628
.500735
.441810
.557636
.500673
.477316
.5 .503704
.499085
.430024
.498608
.497668
.371402
θ'=45 .643963
.499012
.369508
.606474
.500092
.407075
.556388
.500639
.456342
.5 .522801
.498293
.443947
.515052
.496531
.388842
.495691
.495008
.345965
θ'=30 .609394
.498562
.391694
.557954
.499986
.442491
.5 .543655
.498192
.445399
.558434
.495380
.400299
.533668
.493226
.364378
.494600
.492252
.339747
θ'=15 .560499
.499279
.438704
.5 .557552
.498554
.444039
.592813
.495627
.400073
.595139
.492603
.369244
.553905
.490548
.349787
.497659
.489828
.340358
θ'=0 .5 .562178
.498904
.441309
.609785
.496201
.397243
.631434
.493110
.367864
.620457
.490626
.350305
.566078
.489298
.342401
.5
.489038
.341108

We provide an interactive computer program module below where you can reproduce these results. The default values in the NPIB form will generate results highlighted in red in the Table above for θ=30 and θ'=60 which uses the Pipes-Pagano approximation for elastic properties. Unfortunately there are unresolved errors in the Fortran computer program (plysing.f) in module15b that predicts singularities listed in Table 7 above.

Singularites predicted for the stress "Free Edge" problem in section 7.1 have been reproduced, listed in Table 6, and can be used in enriched finite elements in Section 8.0 below.


8.0 Cracks in Homogeneous Isotropic Materials

Outline how enriched finite elements are used to augment existing nonsingular finite elements with singularities that exist at crack tips, Benzley Ref.[29]. Use this enriched FEM to model crack tip singularities in a homogeneous isotropic material and compare numerical results with analytic results that predict Stress Intensity factors for Mode-I and Mode-II, Heyliger & Kriz Ref.[30].

Benzely's Model reproduced here


9.0 Stress Concentrations in Homogeneous Anisotropic Plates

Predict stress gradient near holes in an anisotropic and isotropic plate and compare results with Mizra & Olson Ref.[31].

Mizra & Olson's results reproduced here


10.0 Wave Propagation in Homogeneous Isotropic/Anisotropic Materials

Introduction to 1-D isotropic wave propagation: Analytic and Numerical Models

Introduction to 2-D anisotropic wave propagation: Numerical models



11.0 References

  1. Kriz, R.D., "Influence of Ply Cracks on Fracture Strength of Graphite/Epoxy Laminates at 76K", Effects of Defects in Composite Materials, ASTM STP 836, American Society for Testing and Materials, pp. 250-265, 1984.

  2. Broeck, D., Elementary Engineering Fracture Mechanics, Kluwer Academic Publishers Group, 1982.

  3. Eftis, J., Subramonian, N., and Liebowitz, H., "Crack Border Stress and Displacement Equations Revisited," Engineering Fracture Mechanics, Vol. 9, pp. 189-210, 1977.

  4. Williams, M. L., "The Stresses Around a Fault of Crack in Dissimilar Media," Bulletin of the Seismological Society of America, Vol. 49, No. 2, pp. 199-204, April, 1959

  5. Zak, A. R., and Williams, M. L., "Crack Point Stress Singularities at a Bi-Material Interface," GALCIT SM 62-1, California Institute of Technology, January, 1962.

  6. Nye, J. F., Physical Properties of Crystals, Oxford University Press, London, 1957.

  7. Lekhnitskii, S. G., Theory of Elasticity of an Anisotropic Body, Holden-Day, San Francisco, 1963.

  8. Frederick, D. and Chang, T. S., Continuum Mechanics, Scientific Publishers, Inc. Boston, 1972.

  9. Kriz, R.D. and Ledbetter, H.M., "Elastic Representation Surfaces of Unidirectional Graphite/Epoxy Composites," Recent Advances in Composites in the United States and Japan, ASTM STP 864, J.R. Vinson and M. Taya, Eds. American Society for Testing and Materials, Philadelphia, pp. 661-675, 1985.

  10. Kriz, R.D. and Stinchcomb, W.S., "Elastic Moduli of Transversely Isotropic Fibers and Their Composites," Experimental Mechanics, Vol. 19, pp. 41-49, 1979.

  11. Golfein, S., "Investigation of Structure of Graphite FIbers with the Goal of Improving Their Mechanical Properties," Internal Report AD-A013 524, Army Mobility Equipment Research and Development Center, Fort Belvoir, Virginia, 1974.

  12. Barnet, F.R. and Norr, M.K., "A Three-dimensional Structural Model for a High Modulus Pan-based Carbon Fibre," Composites, pp. 93-99, 1976.

  13. Ledbetter, H.M. and Kriz, R.D., "Elastic Wave Surfaces in Solids," Phys. Stat. Sol. (b) 114, p. 475, 1982.

  14. Musgrave, M.J.P., "On an Elastodynamic Classification of Orthorhombic Media," Proc. R. Soc. Lond. A374, p. 401, 1981.

  15. Kriz, R.D. and Ledbetter, M.M., "Elastic-wave Surfaces in Anisotropic Media," Rheology of Anisotropic Materials, C. Huet, D. Fourgoin, S. Richmond, Eds. CR 19 Coll. GFR Paris Nov. 1984 Editions CEPADUES, Toulouse, pp. 79-92, 1986.

  16. Kriz, R.D. and Heyliger, P.R., "Finite Element Model of Stress Wave Topology in Unidirectional Graphite/Epoxy: Wave Velocities and Flux Deviations,"Review of Progress in Quantitative Nondestructive Evaluation, Vol.8A, Plenum Publishing Corp., pg 141-148, 1989.

  17. Kriz, R.D. and Gary J.M., "Numerical Simulation and Visualization Models of Stress Wave propagation in Graphite/Epoxy Composites,"Review of Progress in Quantitative Nondestructive Evaluation, Vol.9, Plenum Publishing Corp., pg 125-132, 1990.

  18. Vandenbossche, B. and Kriz, R.D.,and Oshima, T." Stress Wave Displacement Polarizations and Attenuation in Unidirectional Composites: Theory and Experiment, "J. Research In Nondestructive Evaluation, Vol. 8, No. 2, pp. 101-123, 1996.

  19. Jones, B., Mechanics of Composite Materials, McGraw Hill, 1975.

  20. Ashton, J.E. and Whitney, J.M., "Theory of Laminated Plates", Technomic Publishing Co., Inc., Stamford Conn., 1969.

  21. Kriz, R.D. and Stinchcomb, W.W., "Effects of Moisture, Residual Thermal Curing Stresses, and Mechanical Load on the Damage Development in Quasi-Isotropic Laminates," Damage in Composite Materials, ASTM STP 775, American Society for Testing and Materials, 1982, pp.63-80.

  22. Kriz, R.D., "Edge Stresses in Woven Laminates at Low Temperatures," Composite Materials: Fatique and Fracture, Second Volume, ASTM STP 1012, Paul A Lagace, Ed., American Society for Testing and Materials, 1989, pp. 150-161.

  23. R. B. Pipes and N.J. Pagano, "Interlaminar Stresses in Composite Laminates Under Uniform Axial Extension", J. Composite Materials, Vol. 4, pp.538-548, 1970.

  24. Stroh, A.N., "Steady State Problems in Anisotropic Elasticity,", J. Mathematics and Physics, Vol. 41, pp. 77-103, 1962.

  25. Ting, T.C.T and Chou, S.C., "Stress Singularities in Laminated Composites," Fracture of Composite Materials, eds. G.C. Sih and V.P. Tamuzs, 1982, Martinus Nijhoff, Proceedins of the 2nd USA-USSR Symposium, Lehigh Univ., Bethlehem, PA, March 9-12, 1981, pp. 265-277, 1982.

  26. Ting, T.C.T. and Hoang, P.H., "Singularities at the Tip of a Crack Normal to the Interface of an Anisotropic Layered Composite," Int. J. Solids and Structures, Vol. 20, No. 5, pp. 439-454, 1984.

  27. Tewary, V. K. and Kriz, R. D., "Generalized Plane strain Analysis of Bimaterial Composite a Free Surface Normal to the Interface", Journal of Material Research, Vol. 13, pp. 2609-2622, 1991.

  28. Tewary, V.K. and R. D. Kriz, "Effect of a Free Surface on Stress Distribution in a Bimaterial Composite", NIST Special Publication 802, National Institute of Standards and Technology, U.S. Gov. Printing Office, Washington D.C., March 1991.

  29. Benzley, S.E., "Representation of Singularities with Isoparametric Finite Elements," Int. J. Numerical Methods in Engineering, Vol. 8, pp. 537-545, 1974.

  30. Mirza, F.A. and Olson, M.D., "The Mixed Finite Element Method in Plane Elasticity," Int. J. Numerical Methods in Engineering, Vol. 15, pp. 273-289, 1980.

  31. Heyliger, P.R. and Kriz, R.D., "Stress Intensity Factors by Enriched Mixed Finite Elements," Int. J. Numerical Methods in Engineering, Vol. 28, 1461-1473, 1989.


Send comments to:
rkriz@vt.edu
Ronald D. Kriz, Short Bio
Engineering Science and Mechanics
College of Engineering
Virginia Tech
Blacksburg, Virginia 24061

Web Content Chronology:
Created ESM Dept On-Campus Web-server November 1997 / Modified December 2000
    http://www.jwave.vt.edu/~rkriz/crcd/lectures/OnePageLect.html
Moved to Off-Campus Web-server November 2014
    http://www.jwave.rkriz.net/crcd/kriz/lectures/OnePageLect.html

DVD: ESM Content Archive:
DVD-ESM-Archive/esm/ESM-DVD/classes/ESM5344/ESM5344_NoteBook/crcd/lectures/OnePageLect.html