### PPT

```Finite Element Analysis
Of Composite Layered Structures
Connor Kaufmann – B. Sc. ‘14
Neola Putnam – M. Eng. ‘14
Ethan Seo – M. Eng. ‘14
Ju Hwan (Jay) Shin – B. Sc. ‘14
Cornell University
Sibley School of Mechanical & Aerospace Engineering
Spring 2014 – Professor N. Zabaras
Objective
C. Kaufmann, N. Putnam, E. Seo, J. Shin
 Develop a linear 3D finite element analysis from scratch using MATLAB.
 Verify the results with expected results, namely the state of out-of-plane stresses near the free edges.
 Observe the effect of h-refinement (convergence of results through mesh refinements).
 Perform a simple sanity check by doing a force reaction balance with the applied traction.
 Compare the numerical results with a commercial FE software, or Ansys Composite PrepPost (ACP).
Stress_yz for [0 90]s laminate
14000
12000

yz [Pa]
10000
0
8000
6000
4000
2000
0
0
20
40
60
Along y-direction [m]
80
100
2
Some stress contours!
C. Kaufmann, N. Putnam, E. Seo, J. Shin

These results refer to a 0°/90°/90°/0° cross-ply laminate.
3
Overview of composite materials
C. Kaufmann, N. Putnam, E. Seo, J. Shin
 Composite materials are commonly used in aerospace structures to minimize mass.
 We considered symmetric, unidirectional, fiber-reinforced composites.
 Composite lamina (sheets) can be stacked to form high strength laminates.
 Laminate stack-ups are characterized by the orientation angles of the fibers, and the materials used.
 Anisotropy of the laminates allows one to tailor designs for stiffness and strength in specific directions.
i.e. 0°/+45°/−45°/90°
s
4
Complexities of layered structures
C. Kaufmann, N. Putnam, E. Seo, J. Shin
 Composite materials can give a much more complicated mechanical response than monolithic materials.
 Stress equilibrium must be satisfied in the laminate by way of interlaminar stresses.
 Special care must be taken to consider the free-edge and free-corner effect in composite samples!
 Non-intuitive effects, such as normal-shear coupling, can occur in anisotropic materials.
 As a result, finite element analysis is often found useful for predicting the behavior of complicated composites.
Normal-shear
coupling
Delamination
5
FEM Formulation
C. Kaufmann, N. Putnam, E. Seo, J. Shin
 Pre-processing:
 Define the required size dimensions, material properties, and the laminate configuration.
 Discretize the model into finite elements.
 Consider a tri-linear hexahedron element.
 Apply any bias factor when discretizing (optional).
 Calculate the 3D elasticity matrix,  .  i.e.   =
 Compute the elemental matrices necessary in developing the stiffness equation.
 Specify the boundary conditions (includes the external load).
 Processing:
 “Globalize” and assemble the local stiffness matrices and the local load vectors.
 Partition and rearrange the global stiffness equation.
 Solve for the nodal displacement field!
 Post-processing:
 Compute the strain field by applying the kinematic equation (displacement ↔ strain).
 Compute the stress field by applying the constitutive equation (Hooke’s Law).
6
Pre-processing (1/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a 8-noded hexahedron, or tri-linear element (3 translational DOF per node).
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 Define the coordinate axes, based on fiber orientations.
 x-axis: longitudinal direction
 y-axis: transverse direction
 z-axis: normal direction (or thickness direction)

 Specify the size dimensions.  , , and
 Specify the material properties.  1 , 2 , 3 , 12 , 13 , 23 , 12 , 13 , and 23
 Specify the fiber orientations of the off-axis plies.
7
Pre-processing (2/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a 8-noded hexahedron, or tri-linear element (3 translational DOF per node).
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 Store the global nodes as local nodes for each element.
 Adhere to the given node-numbering scheme to ensure that the determinant of the Jacobian matrix is positive.
 Guarantee invertible mapping to natural coordinate system.
8
Pre-processing (3/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a tri-linear hexahedron element.
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 Use an eight-noded hexahedron element.
 Each node has three translational degrees-of-freedom.
  ,  , and
9
Pre-processing (4/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a tri-linear hexahedron element.
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 Bias Factor: Allows us to have more concentrated mesh density near a particular region of interest.
 The spacing between the node becomes a geometrical series.
Node Profile
100
80
Finer mesh
y-axis
60
40
20
0
0
100
200
300
x-axis
400
500
10
Pre-processing (5/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a tri-linear hexahedron element.
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 The three-dimensional elasticity (stiffness) matrix is defined by applying the generalized Hooke’s Law.
 Take into account the anisotropy, assuming a transversely isotropic layer.
 Material nonlinearity (plasticity) is neglected!
= 1 −

=

2

=

1 − 23 32
2 3 Δ
21 + 23 31
2 3 Δ
= 31 + 21 32
2 3 Δ
0
0
0
Δ≡
21 + 23 31
2 3 Δ
1 − 13 31
1 3 Δ
32 + 12 31
1 3 Δ
0
0
0
31 + 21 32
2 3 Δ
32 + 12 31
1 3 Δ
1 − 12 21
1 2 Δ
0
0
0
0
0
0
0
0
0
0
0
0
23
0
0
0
13
0
0
0
12
1 − 12 21 − 23 32 − 13 31 − 221 32 13
1 2 3
11
Pre-processing (6/7)







Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a tri-linear hexahedron element.
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).

1
= 0
0
1

0

C. Kaufmann, N. Putnam, E. Seo, J. Shin

≡ s
0
1
0
0
1

0
0
0
1

=
1

1

0
1

0
0
1
0
0
1

1

1

0
2
0
0
0
2
0
2

0
0
0
2
0
2

0
0
0
2

2

2

0
2

⋯ nen
⋯
0
⋯
0
0
0
2

2

2

0
0

nen
0

nen
⋯

⋯
⋯
⋯
0
0
0

nen
0

nen

0
0
0

nen

nen
⋯

nen
⋯

0

nen

=
+1

=
≠
0
0
−

−
0

=
≠
+1
−

−

=
≠
+1
d
=
d
1

− ℎ
+1
−

−
≠ℎ ∧ ≠
ℎ≠

nen

nen

nen

+1
+1
d
=
d
1

− ℎ
+1
−
≠ℎ ∧ ≠

−
ℎ≠
+1
d
=
d
1

− ℎ
ℎ≠
+1
≠ℎ ∧ ≠
−

−
12
−

−
Pre-processing (7/7)







C. Kaufmann, N. Putnam, E. Seo, J. Shin
Define the required size dimensions, material properties, and the laminate configuration.
Discretize the model into finite elements.
Consider a tri-linear hexahedron element.
Apply any bias factor when discretizing (optional).
Calculate the 3D elasticity matrix,  .  i.e.   =
Compute the elemental matrices necessary in developing the stiffness equation.
Specify the boundary conditions (includes the external load).
 Essential Boundary Condition
   = 0, ,  = 0
  ,  = 0,  = 0
  , ,  = 0 = 0
 Natural Boundary Condition
 Pressure-based load, 0 , at  =
13
Processing (1/3)



C. Kaufmann, N. Putnam, E. Seo, J. Shin
“Globalize” and assemble the local stiffness matrices and the local load vectors.
Partition and rearrange the global stiffness equation.
Solve for the nodal displacement field!
 The weak form of our finite element formulation is given below.
 Use Gauss Quadrature rule to numerically evaluate the local integration.
 Apply the transformation rule to  and   using the connectivity matrix.
   glb =  ⊤
   glb =  ⊤
 Sum individual matrices and vectors for global assembly.
  =  glb
  =   glb
nel
nel
+1 +1 +1

⊤
⊤
d d d   =

−1 −1 −1
=1
=1

⊤
⊤
dΓ
Γt

14
Processing (2/3)



C. Kaufmann, N. Putnam, E. Seo, J. Shin
“Globalize” and assemble the local stiffness matrices and the local load vectors.
Partition and rearrange the global stiffness equation.
Solve for the nodal displacement field!
 Partition the global stiffness equation into the known and unknown components.
 Apply the transformation rule to rearrange them as shown below.
=

⊤

=

15
Processing (3/3)



C. Kaufmann, N. Putnam, E. Seo, J. Shin
“Globalize” and assemble the local stiffness matrices and the local load vectors.
Partition and rearrange the global stiffness equation.
Solve for the nodal displacement field!
 Solve the system of equations efficiently by using the Gaussian elimination method.
 In MATLAB, a built-in function, d=K\f can be employed.
⊤
+   =
= −1
, ∀  = 0
=   =  −1

16
Post-processing (1/2)


C. Kaufmann, N. Putnam, E. Seo, J. Shin
Compute the strain field by applying the kinematic equation (displacement ↔ strain).
Compute the stress field by applying the constitutive equation (Hooke’s Law).
 The elemental strain vector can be computed as shown below.
 Weighted average of the strain values, evaluated the Gauss points.
= s

0

0

=
0

0
0

0
0

0

0
17
Post-processing (2/2)


C. Kaufmann, N. Putnam, E. Seo, J. Shin
Compute the strain field by applying the kinematic equation (displacement ↔ strain).
Compute the stress field by applying the constitutive equation (Hooke’s Law).
 The elemental stress vector can be computed as shown below.
 Weighted average of the stress values, evaluated the Gauss points.
=   = 1 −
vm =
1
2
−
2
+  −
2
2
+  −
2

2
2
2
+ 6
+
+
det  −  = 0
3 − 1 2 − 2  − 3 = 0
1 =  +  +
2
2
2
2 =   +   +   −
−
−
2
2
2
3 =    −
−
−
+ 2
1 = p1
max = MAX
2 = p2
3 = p3
p2 − p3 p1 − p3 p1 − p2
,
,
2
2
2
18
Pathwise-results!
C. Kaufmann, N. Putnam, E. Seo, J. Shin
These results refer to a −35°/+35°/+35°/−35° angle-ply laminate.
Stress_yz for [-35 35]s laminate
6
x 10
Stress_xz for [-35 35]s laminate
7
x 10
Stress_yz for [-35 35]s laminate
1
-2000
-0.5
-4000
0.5
0
-1.5
-2
-0.5
-1
-6000
yz [Pa]
0
xz [Pa]
yz [Pa]
-1

5
-3
10
15
20
25
Along z-direction [m]
30
35
0

5
-10000
-12000
Stress
singularity!
-2.5
Ply interface
-8000
-14000
-16000
-18000
10
15
20
25
Along z-direction [m]
30
35
0

20
40
60
Along y-direction [m]
80
100
16800 seconds ≈ 4.5 hours!
19
Contour Plots!
C. Kaufmann, N. Putnam, E. Seo, J. Shin

These results refer to a −35°/+35° s angle-ply laminate.
20
More contours…
C. Kaufmann, N. Putnam, E. Seo, J. Shin

These results refer to a −35°/+35°
s
angle-ply laminate.
21
Sanity Check!
C. Kaufmann, N. Putnam, E. Seo, J. Shin
 A force reaction balance check would indicate that our analysis was modeled correctly!
• Reaction load is extracted from nodes (DOF to be more specific), where an essential BC is specified.
Nodal
Reaction
Forces,
External
Pressure
+   =
=
Input
Laminate width,
Laminate thickness,

.
×
= ext + r = 0

22
Error Analysis
C. Kaufmann, N. Putnam, E. Seo, J. Shin
• L2 and Energy error norms are considered.
• In the below formula,  is equal to 1, since a linear element is considered.

en
L2
=   − h
=   − h

L2
L2
en
=
2
=
1
1
2
2
1
1
− h
2
d
1
−  h
≤ ℎ+1

ℎ=
2
en
2
2
d
≤ ℎ
ℎ2 + ℎ2 + ℎ2
23
Comparison to Ansys?
C. Kaufmann, N. Putnam, E. Seo, J. Shin
Averaged over an element
5.3539e7 Pa
(Ansys ACP)

5.358e7 Pa
(MATLAB)
24
```