Finite Element Analysis of Composite Layered Structures

Finite Element Analysis of Composite Layered Structures

___________________________________

ABSTRACT

A complete formulation of three-dimensional finite element analysis for multi-layered

structures was examined using MATLAB. Throughout the investigation, standard weak

form expressions, model discretization, local node-numbering, global matrix

partitioning, and solution procedures were followed in modeling a uniaxially-loaded,

symmetric laminate. Modeling each discretized element as a tri-linear element, the

effect of mesh refinement was observed to be significant, and agreed well with the

previously hypothesized possibility of a stress singularity. Furthermore, the results

obtained from the three-dimensional analysis were compared with that of Classical

Laminated Plate Theory, as well as ANSYS. For code validation purposes, a forcereaction balance was performed a posteriori to verify that it was implemented correctly.

___________________________________

AUTHORS

Connor Kaufmann

Neola Putnam

Ethan Seo

Ju Hwan (Jay) Shin

___________________________________

Sibley School of Mechanical & Aerospace Engineering

Cornell University, Ithaca, NY, 14853

Submitted on May 1, 2014 to Prof. N. Zabaras

Shin et al.

1. Introduction

more detail in the following paragraphs, the nodal

displacements, strain field, and the stress field were outputted.

Composite materials are known to be more favorable for

optimizing weight of the structure than monolithic materials.

This is due to the anisotropic nature of the fiber-reinforced

laminate, which allows engineers to utilize the composite

materials without overdesigning in certain orientations. While

these multidirectional laminates were shown to be effective at

reducing the overall payload of the structure, they introduced a

new level of complexity in the proper understanding of their

mechanics and led to investigative research in this field.

In the pre-processing stage, the model was discretized

into finite elements. In doing so, bias-factor was considered to

allow for a greater mesh density (Fig. 2) near the free-edge.

This optional step is helpful for obtaining more accurate results

without having to increase the number of elements.

A tri-linear element was employed throughout the project.

A first-order element was considered the easiest to implement

among the general p-order Lagrangian elements. The element

used was an eight-node hexahedron element (Fig. 3) with the

local node-numbering scheme as given, and each node

contained three translational degrees-of-freedom: ?, ?, and ?.

The appropriateness of the node-numbering convention was

confirmed by checking that the determinants of the elemental

Jacobian matrices were positive. This ensured that the mapping

between physical and natural coordinates was correct.

One particular area of such research is devoted to the

study of failure by delamination (Fig. 1). Of the many

mechanisms by which a layered structure fails due to

delamination, one notable cause of this is the development of

interlaminar stresses. Interlaminar stresses are out-of-plane

stresses, such as ??? , ??? , and ??? , that are induced at the ply

interface near the free-edges of a laminated structure. Though

completely neglected in the framework of 2-D CLPT, or Classical

Laminated Plate Theory, these through-thickness stresses can

be quite severe in some cases, as there can exist a state of stress

singularity from the abrupt mismatch of the Poisson¡¯s effects

along the layers. The present study sought to confirm this

behavior numerically by performing a complete threedimensional finite element analysis, as the Kirchhoff-Love¡¯s

assumptions used in CLPT was limited to estimating the

membrane stresses only, or in-plane stresses. It is noted,

however, that to model the physical mechanics with even

greater accuracy, one should take into account the material

nonlinearity in the analysis. Material nonlinearity plays an

important role in failure by delamination because failure is

often not governed only by the magnitude of the interlaminar

stress (which is infinite in linear elastic theory). Nonlinearity of

material properties causes the crack tip, formed at the onset of

progressive fracture to experience a blunting effect at the

material¡¯s yield point. Nevertheless, it is worth acknowledging

that this singular stress is important for understanding the

failure mechanism of a laminated structure.

Moreover, the elasticity, or the stiffness matrix, which

relates the strain to stress, was also defined for each layer of the

? ]{?}, is also

laminate. This constitutive relationship, {?} = [?

referred to as the generalized Hooke¡¯s Law (App. A-1).

The rest of the components required to formulate the

elemental stiffness equation were calculated in the preprocessing stage. These included the [?? ], [?? ], [?? ], and [?? ]

matrices evaluated at the Gauss integration points. Here it is

noted that [?? ] is the shape function matrix, [?? ] = [?s ?? ] is

the symmetric gradient of [?? ], and [?? ] and [?? ] are referred

to as the Jacobian matrix and the connectivity matrix,

respectively (App. A-2).

The last step of the pre-processing stage involved

specifying the boundary conditions. In the case of uniaxiallyloaded laminate, three symmetry planes within the model were

observed, such that only one octant of the model was analyzed.

The essential boundary conditions were that the appropriate

degrees-of-freedom on each symmetry plane were set to zero

directional displacements. In addition, the applied pressure

load boundary condition was specified as follows.

Ply interface

?(? = 0, ?, ?) = 0

?(?, ? = 0, ?) = 0

?(?, ?, ? = 0) = 0

??(? = ?, ?, ?) = ?0

Fig. 1: Failure due to delamination

2. FEM Formulation

In developing the FEM code, the analysis was divided into

three main stages: pre-processing, processing, and postprocessing. Using MATLAB, any necessary inputs, including the

anisotropic material properties, the applied load, the laminate

configuration, and the laminate dimensions, were specified.

Following the standard solution procedures to be discussed in

Fig. 2: Mesh profile of a laminated structure

-2-

Shin et al.

?

Similarly, the stress field was calculated by applying

Hooke¡¯s Law using the [?? ] matrices that were previously

computed.

?

{?} = [?? ]{?}

While there are other ways of displaying the strain/stress

fields, such as an interpolated field that would result in a smooth

variation over the model, as opposed to a discontinuous set of

elemental results, the latter was chosen for its slightly greater

efficiency. Furthermore, a number of useful features were built

into the code, including the options to show the deformed

model, display the elemental boundaries, and generate a

comprehensive solution file.

?

Fig. 3: 8-node hexahedron element with 3 degrees-of-freedom per node

Once the pre-processing stage was completed, the

elemental stiffness equations were assembled to construct the

global stiffness equation. Recalling that the weak form of the

governing equation is as shown below, the stiffness matrix, [?? ],

and the load vector, {? ? }, were simply summed for all elements.

3. Numerical Results

Several salient path-wise and contour plots (App. B) are

illustrated in this section. A number of observations with

emphasis on the free-edge effects are also briefly mentioned.

nel

+1 +1 +1

¡Æ (?? ?

??=1

¡Ò ¡Ò ¡Ò ?? ? ?? ?? |?? | d? d? d? ?? ) ?

?1 ?1 ?1

?

nel

= ¡Æ (?? ? ¡Ò ?? ? ?? ? d¦£)

?

¦£?t

?=1

?

To solve for the unknown components of the global

degrees-of-freedom, {?? } , the global stiffness equation was

rearranged and partitioned by applying the appropriate

transformation, as shown below.

???

?? = ?

??

[ ?

???

(Free-edge)

??? ??

?

] [ ] = [ ?]

??

?? ??

Fig. 4: Path-wise distribution of ??? along y

?

?? = ??1

? ?? , ? ?? = 0

It is noted that the method of Gaussian elimination

implemented in MATLAB was used to solve the above set of

equations; this is done by setting d_F = K_F \ f_F for greater

efficiency than direct matrix inversion. The known components

of the degrees-of-freedom, {?? } were appended to the solved

solutions, or {?? }, in order to construct the full set of solutions.

Possibility

of stress

singularity

Finally, the strain field and the stress field were calculated

in the post-processing stage. Using the kinematic relationship,

the element-averaged strain was calculated, as follows.

{?} = [?? ]{?},

where [?? ] was computed as the weighted average of the

values obtained at each Gauss integration point.

(Ply interface)

???

Fig. 5: Path-wise distribution of ??? along z

-3-

Shin et al.

???

Fig. 9: ??? distribution along ? for various mesh profiles

Fig. 6: ??? distribution along ? for various mesh profiles

???

Fig. 7: ??? distribution along ? for various mesh profiles

Fig. 10: ??? distribution along ? for various mesh profiles

???

Fig. 8: ??? distribution along ? for various mesh profiles

Fig. 11: ??? distribution along ? for various mesh profiles

-4-

(Free-edge)

(Free-edge)

???

(Free-edge)

(Free-edge)

???

(Free-edge)

(Free-edge)

???

Shin et al.

It is observed (Fig. 4-11) that the magnitude of ???

begins to decay rapidly near the free-edge, which is consistent

with the behavior for a traction-free surface. Similarly, ??? ,

which is assumed to be zero for a cross-ply laminate under the

framework of CLPT, displayed a gradual decay toward the freeedge. On the other hand, ??? exhibited a possible state of stress

singularity. It is noted that this observation was consistent with

findings from prior research.

orthogonal to the direction of the applied load. The table above

summarizes a sample computation. Since the applied load is

equal and opposite to the corresponding reaction traction, their

sum is zero and the force balance is satisfied.

A simple error analysis is shown in the following

paragraph. Due to the difficulty in solving for the exact solution

for a layered structure, the theoretical L2 and energy norm

errors were computed to demonstrate the error convergence

rate with mesh refinement.

Furthermore, the validity of CLPT for various laminates

was investigated. The table below indicated that CLPT was an

accurate assumption only for a laminate with small mismatch in

Poisson effect between the neighboring plies. Moreover, it was

noted that the mesh refinement had a significant effect on the

axial stress results, especially for a laminate with a large

mismatch in fiber orientation.

2

||?||L = (¡Æ ¡Ò (?(?) ? ?(? )) d¦¸)

2

||?||en

(%)

nel

¦¸?

1?

2

¡Ü ???+1

1

= (¡Æ ¡Ò ? ? (?(?) ? ? ? (?)) d¦¸)

2 ¦¸?

1?

2

¡Ü ???

nel

Here, it is noted that ? is a constant that is dependent on

the interpolation function used for the FE solution, ? refers to

the order of the element (? = 1 for a linear element), and ? ¡Ô

¡Ì??2 + ??2 + ??2 is the largest diameter (Fig. 13) that can

inscribed in each element.

4. Further Verification

Fig. 13: Schematic illustrating element size, ?

After the model was implemented in MATLAB, a forcereaction balance test was conducted as a sanity check to confirm

that the code had solved the mathematical model correctly. The

configuration being modeled was in static equilibrium, and as

such, summation of forces in each coordinate direction should

be equal to zero. The following figure (Fig. 12) shows how the

reaction forces at nodes with essential boundary conditions

specified should equal the externally applied traction forces.

fr

Thus by reducing the element parameter, ?, by a factor of

2, the L2 and energy errors are expected to be reduced to ? and

? the original values, respectively. As successively refined

meshes were tested with the developed code, mesh

convergence was observed, and it could be reasonably be

assumed that the FE solution approached the exact solution.

Lastly, the effect of varying the aspect ratio of an element

was also investigated in the present study. Here, the aspect

ratio, ? ¡Ô ?max /?min , is defined as the ratio of the largest

dimension that can be inscribed in the element to the smallest

dimension. While the reasoning behind how the aspect ratio

has an effect on the solution will not be explained in detail, it is

noted that this is due to the shear locking phenomena. The firstorder element¡¯s inability to capture the effect of large shear

deformation leads to an inaccurate result as the stiffness will be

too high. Thus it is recommended that the element be

discretized appropriately so as to attain a moderate aspect ratio.

fext

Fig. 12: Free-body-diagram of the force-reaction balance

Laminate width, ?

Laminate thickness, ?

Applied load, ?0

Total reaction force

Normalized traction

10 m

3.6 m

5 ¡Á 107 Pa

-1.8 ¡Á 109 N

-5 ¡Á 107 Pa

¡Æ ?? = ?ext + ??

0

5. Concluding Remarks

This paper discusses the formulation methodology for

implementing a 3-D finite element analysis of a layered

structure using MATLAB. Following the standard procedures of

FEM for an eight-node tri-linear element, several path-wise and

contour plots were examined. The results exhibited our

expected behavior for a uniaxially-loaded, symmetric laminate.

The reaction traction was computed by dividing the

external force by the cross-sectional area of the plane

-5-

 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download