Report

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. Consider a uniaxial loading of a symmetric laminate. 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! External load ≟ Reaction load • Reaction load is extracted from nodes (DOF to be more specific), where an essential BC is specified. Nodal Reaction Forces, External Pressure Load, 0 + = = Input Laminate width, Laminate thickness, Pressure Load, 0 . × = 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