Microstructure Lectures

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
   • 4.1 Introduction to Mechanical Behavior of Anisotropic Materials
   • 4.2 Introduction to Mechanical Behavior of 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, 10⁰m. 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 10⁰m

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 0^o 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 0^o 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.

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.

(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 -a 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 C_I, C_II, and C_III are functions of Stress Intensity and radius "r", and the variables F_ij^I, F_ij^III, and F_ij^III are functions of geometry only, i.e. location of boundaries and the angle, q, 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.

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, l on the radius "r". The exponent, l 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, l

---

Mechanical behavior of anisotropic materials was already introduced in the continuum (macroscale) section on Material Characterization where elastic fourth order stiffness tensor, C_ijkl, was contracted into a six by six stiffness matrix, C_ab, relating six components, s_a, of a symmetric stress tensor to the six components, e_b, 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 C_ijkl, where the stresses, s_ij, and Lagrangian strains, l_ij, 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, l_ij.

(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, r is the density, e_ij is the Eulerian strain associated with Eulerian displacement gradients, c is the radiative heat transfer per unit time per unit mass, and b_i 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, l_ij, 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)

(11)

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

(12)

If this relationship is true for any arbitrary nonzero dl_ij, 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 a 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 b_ij represent the residual stresses when the strains, l_ij, go to zero, and 2) the coefficient d_ijklmn is the nonlinear elastic term that gives rise to nonlinear 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 C_ijkl 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 C_ijkl transforms as a fourth order tensor.

(18)

(19)

(20)

---

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:

(a₁₁, a₁₂, a₁₃) = direction cosines of the y^'₁-axis with respect to the (y₁,y₂,y₃)-axes.

(a₂₁, a₂₂, a₂₃) = direction cosines of the y^'₂-axis with respect to the (y₁,y₂,y₃)-axes.

(a₃₁, a₃₂, a₃₃) = direction cosines of the y^'₃-axis with respect to the (y₁,y₂,y₃)-axes.

Figure 9. Direction Cosine Matrix With a Simple Rotation About the y₃ 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, a_ij, 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, C_ijkl .

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 S_ij and C_ij.

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, a_ij, would be the identity matrix but with the term a₃₃ 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 s₁ with s₁^' requires C₁₄ and C₁₅ 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 n, 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.

Another linear constitutive relationship, similar to equation (24), exists for compliances, S_ijkl, 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, S_ij and C_ij 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), S₃₃ 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, a_ij, are expanded into sines and cosines of the rotation angle, q, starting at zero along the z axis and rotating 90^o 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.

Material	C33	C11	C12	C13	C44	n32	n12
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

Material	Poisson's Ratio(n)    Modulus: Bulk(K), Shear(G), Young's(E) Subscripts: Matrix(m), Fiber-Axis(L), Transverse-Plane(T)
Isotropic: two independent properties->	Em	nm
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 E₃ 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 G₂₃ 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 n ₃₂ 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, n, 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 x₃ 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 x₃ axis because for hexagonal symmetry the x₃ axis is the conventional "unique" symmetry axis. However in engineering design the x₁ axis is the 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.

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, n, although interesting, can only reveal variations of just that one property. The modulus E₃ = S₃₃^-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.

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, n_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, n_l, and particle vibration direction, a_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, k_i, are terms used to define a plane (this will perhaps require a review of basic descriptive geometry in freshman calculus) and the circular frequency, w, 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, n_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, rv², are the eigen-values and the displacement direction cosines, a_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, n_i. Hence we conclude that if we solve for eigen-values (wave speeds) in all possible propagation directions, n_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.

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.

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.

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 a section of materials for each of these symmetries.

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 (15 min.) or large (4.5 hrs.) 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 25 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 25. Lamina axes defined

(44)

The elements of [Q] are related to the engineering constants E,G, and n 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 26.

Figure 26. Lamina-Laminate Configuration Geometry Defined

The laminated plate shown in Figure 27 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 27. Thin Plate Deflections

From Figure 27 it is not difficult to see how in-plane u_o and v_o 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, r, we derive expressions for curvature, k_x and k_y assuming small angular deflections.

(57)

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

Figure 28. 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 (N_x, N_y, N_xy), and moments (M_x, M_y, M_xy).

Figure 29. 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 30 shows the relative position of each k-th lamina used by equation (61) to integrate through the thickness to calculate the laminate loads, (N_x, N_y, N_xy), and moments (M_x, M_y, M_xy).

Figure 30. 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 ( e^o_x, e^o_y, g^o_xy ) and curvatures ( k_x, k_y, k_xy ) 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 {e} and inplane moments {M} only produce
curvature {k}

Case II: Specially Orthotropic Symmetric laminates, A₁₆=0 and A₂₆=0

Bi-directional laminates consist only of 0/90 degree lamina. It is not difficult to show Q₁₆ and Q₂₆
are zero. Hence A₁₆ and A₂₆ are also zero and consequently the inplane shear strain {g} 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, D₁₆ & D₂₆ -> 0 , no. of lamina -> large

For a large number of symmetric pairs, +q at +h and -q at -h, D₁₆ and D₂₆ go to zero. Hence
moments {M_x} & {M_y} and torque {M_xy} 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.

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 31 shows a typical sequence of damage events that lead to ultimate fracture at the macroscale.

Figure 31. 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:

• Theory and formulation of the "Free-Edge" problem.
• Comprehensive description, development, and user's guide of a finite element model (FEM) ( Computer Program) of the "Free-Edge" problem.

• Interactive computer program: Experiment with Free-Edge FEM program: [0/90/45/-45]_s
• Interactive computer Program: Experiment with edge stresses in a [0/90/]_s laminate:
                                                      Compare with Pipes and Pagano (1970) Finite Difference Model

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

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

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 33. Sample results of laminate edge-stress singularities, Ref.[1].

Figure 34. Sample results of stresses at laminate stress free edge, Ref.[1].

You might want to reproduce the results shown in this Figures 33 and 34 by going to the interactive computer program. You will also need to use the Laminated Plate Analysis program which is used to convert the N_x=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 sigma-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, Delta-T=-180F, and swelling caused by moisture absorption, Delta-M=1.2% moisture.

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

In the section above the model predicts edge singularities in a eight-layered quasi-isotropic laminate. For comparison we provide models that predict edge stress states for simple [0/90]_s laminates with and without woven macro-structures, For the nonwoven [0/90]_s laminate a classic solution by Pipes and Pagano, Ref. [23], (Finite Difference Model) is compared with this solution to check the accuracy of our model, but here we are interested in comparing edge stresses with and without woven geometry.

• Interactive Computer Program: Predict stresses in a [0/90/45/-45]_s thin laminate
• Interactive Computer Program: Predict stresses in a [0/90/+45/-45]_s thick laminate prior to 90 ply crack
• Interactive Computer Program: Predict stresses in a [0/90/+45/-45]_s thick laminate near a 90 ply crack

• Case study on influence of ply-cracks on [0/90/+45/-45]s laminate fracture strength:
  Compares stress analysis (FEM program) with fractographs and statistical analysis of fracture strengths.

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 ply crack in a [0/90/+45/-45]_s laminate and the stress gradients ("concentrations") in the load bearing 0 degree ply near the ply crack tip, see also Figures 3 and 4.

Figure 28. Schematic and FEM Grid of Ply Crack in a [0/90/+45/-45]_s laminate

Figure 29. Mode I stresses in the 0 degree load bearing ply above the ply crack.

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 "Free-Edge" problem the model above predicts stress singularities of a ply crack in a eight-layered quasi-isotropic laminate. Below we introduce an FEM model for a simpler [0/90]_s laminate with and without a woven macro-structure.

• Interactive Computer Program: Experiment with lamina stress in a [0/90/]s laminate prior to 90 ply crack
• Interactive Computer Program: Experiment with singularities in a [0/90]s laminate near a 90 ply cracks

• Interactive Computer Program: Experiment with stresses in a 0/90 woven laminate prior to fill crack
• Interactive Computer Program: Experiment with stresses in a 0/90 woven laminate near a fill crack

• Interactive Computer Program: Experiment with singularities at Interfaces Between Dissimilar Anisotropic Materials: Analytic Model (Stroh's Method)

• Interactive Computer Program: Experiment with singularities at cracks in homogeneous isotropic materials: Enriched Finite Element Model (Benzley)

• Interactive Computer Program: Experiment with stress concentrations in homogeneous anisotropic materials: Finite Element Model (Mizra & Olson)

• Interactive computer program: Experiment with 2-D grid force step function.

Introduction to 2-D anisotropic wave propagation: Numerical models

• Interactive computer program: Experiment with small grid FEM.

• Interactive computer program: Experiment with large grid FEM.