Scott E. Stapleton and Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

The Analysis of Adhesively Bonded Advanced

Composite Joints Using Joint Finite Elements

NASA/CR—2012-217606

April 2012

Notice for Copyrighted Information 

This document contains material copyrighted by the party submitting it to NASA—see the copyright notice affi xed thereto. 

It may be reproduced, used to prepare derivative works, displayed, or distributed only by or on behalf of the Government 

and not for private purposes. All other rights are reserved by the United States Government. 
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to help@

sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443–757–5803

• Telephone the NASA STI Help Desk at

 443–757–5802

• Write to:

 NASA Center for AeroSpace Information (CASI)

 7115 Standard Drive

 Hanover, MD 21076–1320
Scott E. Stapleton and Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

The Analysis of Adhesively Bonded Advanced

Composite Joints Using Joint Finite Elements

NASA/CR—2012-217606

April 2012

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared under Contract NNX09AM29H

Prepared under Grant NCC3–989

Notice for Copyrighted Information 

This document contains material copyrighted by the party submitting it to NASA—see the copyright notice affi xed thereto. 

It may be reproduced, used to prepare derivative works, displayed, or distributed only by or on behalf of the Government 

and not for private purposes. All other rights are reserved by the United States Government. 
Acknowledgments

Portions of this work were fi nancially supported by the Space Vehicle Technology Institute under grant NCC3–989 jointly funded 

by NASA and the Department of Defense. The bulk of the fi nancial support was provided by NASA Glenn Research Center 

through the GSRP Fellowship, with Dr. Brett Bednarcyk as mentor.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Trade names and trademarks are used in this report for identifi cation 

only. Their usage does not constitute an offi cial endorsement, 

either expressed or implied, by the National Aeronautics and 

Space Administration.

Level of Review: This material has been technically reviewed by NASA expert reviewer(s).

Copyright © 2012 by Scott E. Stapleton
NASA/CR—2012-217606 iii 

Contents 

Chapter 1 Introduction .................................................................................................................................. 1 

1.1 General Comments About This Thesis ............................................................................................. 1 

1.2 Introduction ...................................................................................................................................... 1 

1.3 Thesis Organization .......................................................................................................................... 1 

1.4 Significant Contributions.................................................................................................................. 2 

1.5 Publications ...................................................................................................................................... 2 

Chapter 2 Linear Elastic Modeling of Adhesively Bonded Joints Using an Enhanced Joint Finite 

Element ................................................................................................................................................... 3 

2.1 Introduction ...................................................................................................................................... 3 

2.2 Analytical Formulation ..................................................................................................................... 4 

2.2.1 Adherend Models .................................................................................................................... 4 

2.2.2 Adhesive Models .................................................................................................................... 6 

2.2.3 General Formulation ............................................................................................................. 10 

2.3 Validation and Results .................................................................................................................... 14 

2.3.1 Adhesive Model Comparison ............................................................................................... 14 

2.3.2 Composite Adherends ........................................................................................................... 22 

2.3.3 Tapers and Steps ................................................................................................................... 23 

2.4 Conclusions .................................................................................................................................... 29 

Chapter 3 Functionally Graded Adhesives for Composite Joints ............................................................... 31

3.1 Introduction .................................................................................................................................... 31 

3.2 Method ............................................................................................................................................ 33 

3.2.1 Formulation ........................................................................................................................... 33 

3.3 Results and Discussion ................................................................................................................... 36 

3.3.1 Baseline Configuration ......................................................................................................... 36 

3.3.2 Stress Comparison ................................................................................................................ 37 

3.3.3 Sensitivity Study ................................................................................................................... 40 

3.3.4 Effect of Loading .................................................................................................................. 42 

3.3.5 Proof-of-Concept Testing ..................................................................................................... 44 

3.4 Concluding Remarks ...................................................................................................................... 48 

Chapter 4 Progressive Failure of Adhesively Bonded Composite Joints .................................................... 51 

4.1 Introduction .................................................................................................................................... 51 

4.2 Formulation .................................................................................................................................... 51 

4.2.1 Co-rotational Formulation .................................................................................................... 51 

4.2.2 Material Nonlinearities ......................................................................................................... 59 

4.2.3 Crack Growth ........................................................................................................................ 60 

4.2.4 Adaptive Shape Functions .................................................................................................... 62 

4.2.5 Adhesive Model Characterization ......................................................................................... 63 

4.3 Results and Validation .................................................................................................................... 68 

4.3.1 Geometric Nonlinearities ...................................................................................................... 68 

4.3.2 Material Nonlinearities ......................................................................................................... 72 

4.3.3 Crack Growth ........................................................................................................................ 74 

4.3.4 Adaptive Shape Functions .................................................................................................... 75 

4.3.5 Experimental Comparison .................................................................................................... 76 

Chapter 5 Conclusion .................................................................................................................................. 85 

5.1 Brief Summary ............................................................................................................................... 85 

5.2 Future Work ................................................................................................................................... 86 

5.2.1 3-D Plate Joint Element ........................................................................................................ 86 

5.2.2 Joint Element Extensions ...................................................................................................... 87 
NASA/CR—2012-217606 iv 

5.2.3 Functionally Graded Adhesives ............................................................................................ 87 

5.2.4 Adherend Yielding ................................................................................................................ 87 

References ................................................................................................................................................... 89 

List of Tables 

Table 2-1. Parameters of in-depth study cases. ........................................................................................... 15 

Table 2-2. Material properties used in Hypersizer analysis [22]. ............................................................... 22 

Table 2-3. Layup and boundary conditions for Hypersizer analysis [22]. .................................................. 22 

Table 3-1. Parameters of the baseline single strap joint configuration used for the theoretical study. ....... 36 

Table 3-2. Geometric parameters for proof-of-concept testing. ................................................................. 46 

Table 4-1. For an unbalanced joint, the peel to stress ratio is different on each side and produces a 

different strength prediction. .......................................................................................................... 75 

Table 4-2. Material properties and geometric parameters for Song/Waas [70] DCB specimens. .............. 77 

Table 4-3. Material properties and geometric parameters of DCB specimens. .......................................... 79 

Table 4-4. Material properties of the single lap joint adherends and adhesive [55]. .................................. 81 

Table 4-5. Experimental and predicted strengths of the single lap joint. .................................................... 82 

List of Figures 

Figure 2-1. Geometric parameters for overlap region of an adhesively bonded joint with multiple bonded 

layers and the width of the joint in the y-direction is b. ................................................................... 4 

Figure 2-2. Boundary conditions for adherend i: prescribed nodal displacements and rotations at x=0 and 

x=l. ................................................................................................................................................. 11 

Figure 2-3. The building block approach facilitates modeling complex joints with simple joint element 

building blocks. .............................................................................................................................. 13 

Figure 2-4. Joints connected together with an offset between centerlines are related using Eq. (21). ....... 14 

Figure 2-5. Boundary conditions, loading scenario and geometric parameters for the single lap joint. 

Material properties include: E1=E2= E, Ea, and Ga. ....................................................................... 15 

Figure 2-6. Relative difference in maximum adhesive centerline (za1=0) shear and peel stress between 

Models 3 and 2. .............................................................................................................................. 16 

Figure 2-7. Comparison of the maximum shear and peel stress predicted by Models 2 and 1 for different 

joint configurations, along with points indicating the cases for comparison with 2-D FEMs. ...... 16 

Figure 2-8. Reaction force predicted by the three adhesive models compared to the 2-D element dense 

mesh model for single lap joins with different geometries. ........................................................... 17 

Figure 2-9. Adhesive centerline stresses predicted by different models for Case 1 parameters. ................ 17 

Figure 2-10. Adhesive stress contours predicted by different models for Case 1 parameters. ................... 18 

Figure 2-11. Adhesive centerline stresses predicted by different models for Case 2 parameters. .............. 19 

Figure 2-12. Adhesive stress contours predicted by different models for Case 2 parameters. ................... 19 

Figure 2-13. Adhesive centerline stresses predicted by different models for Case 3 parameters. .............. 20 

Figure 2-14. Adhesive stress contours predicted by different models for Case 3 parameters. ................... 20 

Figure 2-15. Adhesive centerline stresses predicted by different models for Case 4 parameters. .............. 21 

Figure 2-16. Adhesive stress contours predicted by different models for Case 4 parameters. ................... 21 

Figure 2-17. Geometric parameters for composite stiffened panel comparison with HyperSizer [22]. ..... 22 

Figure 2-18. Comparison of normal stresses in the x-direction of the mid-joint cross section of adherend 1 

of a stiffened panel found using HyperSizer [22] and the joint element........................................ 23 

Figure 2-19. Geometric and material parameters of the stepped double-strap joint. .................................. 23 
NASA/CR—2012-217606 v 

Figure 2-20. Two models compared in the study, (a) joint element model using 26 elements and (b) 

Abaqus 2-D plane stress model using 99,000 elements. ................................................................ 24 

Figure 2-21. Stresses in the (a) upper and (b) lower adhesive centerlines, and (c) left and (d) right 

adherends in the double-strap joint. ............................................................................................... 25 

Figure 2-22. Stresses in the (a) lower and (b) upper surfaces of the top doubler, and (c) lower and (d) 

upper surface of the bottom doubler in the double-strap joint. ...................................................... 26 

Figure 2-23. Geometric and material parameters of the tapered joint. ....................................................... 27

Figure 2-24. Models compared in the study: joint element using (a) 7, (b) 12, (c) and 22 elements, and (d) 

Abaqus 2-D plane stress model using 120,000 elements. .............................................................. 27 

Figure 2-25. Stresses in the (a) left and (b) right adherends and (c) peel and (d) shear stress in the 

adhesive of the tapered joint. ......................................................................................................... 28 

Figure 3-1. Method of linear segments involves breaking up a function into small segments and assuming 

the function is constant within the segment. .................................................................................. 33 

Figure 3-2. A comparison of solutions to the homogeneous, linear, non-constant coefficient differential 

equation 3

, ( sin 4 ) xx y x e xy   with initial conditions y(0) 1  and , (0) 0 x y  (a) using the method 

of constant segments with two, four, and six segments and (b) the relative error as a function of 

the number of segments used. ........................................................................................................ 35 

Figure 3-3. Single strap joint (a) geometric and material parameters and (b) joint finite element 

representation assuming symmetry. ............................................................................................... 36 

Figure 3-4. Single strap joints with (a) constant modulus adhesives were compared with joints with 

functionally graded adhesives, including (b) step-wise graded, (c) linearly graded, and (d) 

exponentially graded. ..................................................................................................................... 37 

Figure 3-5. Optimized configurations for the single strap joints for different functions of graded adhesive 

compared: (a) modulus across the adhesive and (b) centerline peel stress across the adhesive. ... 38 

Figure 3-6. Comparison of the maximum stresses in the joints with different adhesives, where all 

maximum stress values are normalized by the maximum stress value found in the stiffer constant 

modulus (single adhesive) joint, and adherend and doubler stresses refer to those in the xdirection. ........................................................................................................................................ 39 

Figure 3-7. The maximum stresses in joints with linear FGAs and different lower bound modulus, El, 

where all maximum stress values are normalized by the maximum stress value found in the stiffer 

constant modulus (single adhesive) joint, and adherend and doubler stresses refer to those in the 

x-direction. ..................................................................................................................................... 39 

Figure 3-8. Increasing the number of adhesives used for a discretely graded adhesive and optimizing the 

adhesive spacing lowers the maximum adhesive peel stress; (a) stress vs. number of adhesives, 

and (b) modulus at the end of the adhesive for 20 steps compared with the other functions. ....... 40 

Figure 3-9. Observing the maximum adhesive (a) peel and (b) shear stress as a function of l as defined in 

Figure 3-4 shows how sensitive the maximum stress is to the shape of the grading. .................... 41 

Figure 3-10. The sensitivity of the maximum adhesive peel and shear stress to decreasing the lower bound 

modulus, El, for the discrete (a and b), linear (c and d), and exponential (e and f) FGAs. ............ 42 

Figure 3-11. For different loading scenarios at the end of the joint ( (a) axial load, (c) moment, and (e) 

vertical load) the stress gradients of the adhesive peel stress (b, d, and f) remains the same 

direction; high at the ends and low in the middle. ......................................................................... 43

Figure 3-12. Maximum peel and shear adhesive stress and its sensitivity to grading parameter, l, for 

different amounts of end moment (a, b) and end shear force (c, d) for joints with a linear FGA. . 44 

Figure 3-13. For different loading scenarios for a (a) single lap joint, (b) bonded doubler, and (c) double 

strap joint, the stress gradients of the adhesive peel stress (b, d, and f) remains the same direction; 

high at the ends and low in the middle. ......................................................................................... 45 

Figure 3-14. Diagram of single strap joint with placed glass beads to test functionally graded adhesive 

concept. .......................................................................................................................................... 45 
NASA/CR—2012-217606 vi 

Figure 3-15. Photographs of the adhesive layer before placement of the doubler for the single strap joint 

with (a) no beads, (b) uniform beads, and (c) graded beads. ......................................................... 46 

Figure 3-16. Single strap joint test: (a) photograph of typical specimen prior to failure and (b) joint 

strengths of the specimens tested with different adhesive systems. ............................................... 47 

Figure 3-17. Post mortem photographs of the failed adhesive layer for the single strap joints with (a) no 

beads, (b) uniform beads, and (c) graded beads. ............................................................................ 47

Figure 4-1. The nodal displacements can be broken up into two parts: a) rigid translation and rotations and 

b) local deformations. .................................................................................................................... 52 

Figure 4-2. Displacements, ir qrot , of node ir due to a rigid body rotation of the joint element................... 53 

Figure 4-3. The initial and current lengths of the 1st adherend are used to determine the rotation angle. . 55 

Figure 4-4. Diagram showing a) an uncracked joint element, b) a partially cracked element, and c) a fully 

cracked joint element. .................................................................................................................... 60 

Figure 4-5. Flow chart showing how cracked element Sub-Assembly is incorporated into joint element 

solution procedure. ......................................................................................................................... 61 

Figure 4-6. For a single lap joint with (a) nonlinear material properties, the (b) tangent modulus of the 

adhesive layer becomes similar to a functionally graded adhesive in advanced stages of loading.62 

Figure 4-7. Adhesive may be characterized by (a) experimental bulk adhesive tensile tests, then (b) fitting 

a curve to the stress-strain plot. ...................................................................................................... 63 

Figure 4-8. Assuming that the adhesive is perfectly bonded to the adherends, the adhesive can be 

considered a constrained body under triaxial stress. ...................................................................... 64 

Figure 4-9. The peel to shear ratio for adhesive i can be approximated by dividing the maximum peel 

stress by the maximum shear stress for the linear elastic adhesive case. ....................................... 66 

Figure 4-10. Fracture mechanics properties such as critical stress and fracture toughness can be used to 

form an adhesive stress-strain law for the joint element. ............................................................... 67 

Figure 4-11. Geometric dimensions for aluminum beam problem for validation of the co-rotational 

formulation. .................................................................................................................................... 68 

Figure 4-12. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the 

cantilever beam, along with (b) a convergence study for the joint element model. ....................... 69 

Figure 4-13. Load vs displacement plots for a beam loaded in compression with a slight (a) moment and 

(b) imperfection angle. ................................................................................................................... 69 

Figure 4-14. Single lap joint used to validate joint element co-rotational formulation. ............................. 70 

Figure 4-15. Comparison of (a) joint element and 2-D dense mesh finite element representation of the 

joint in Figure 4-14 and (b) the resulting load vs displacement plot. ............................................ 70 

Figure 4-16. (a) Geometric parameters and (b) buckled shape of layered beam under compressive loading.70 

Figure 4-17. Load vs. displacement plot comparing a dense, 2-D finite element mesh with the joint 

element using various adhesive and adherend models. .................................................................. 71 

Figure 4-18. Load displacement plots for the layered beam under compressive loading showing the effects 

of varying the adhesive (a) Young's modulus and (b) shear modulus illustrating the importance of 

the shear modulus on the buckling load. ........................................................................................ 72 

Figure 4-19. Stress vs. strain approximation for aluminum beam depicted in Figure 4-11. ....................... 72 

Figure 4-20. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the 

cantilever beam with nonlinear materials, along with (b) a convergence study for the joint 

element model. ............................................................................................................................... 73 

Figure 4-21. Plots showing the effect of element size on the load-displacement response of the single lap 

joint featured in Figure 4-14with nonlinear adherends for the (a) joint element, and (b) dense 2-d 

finite element mesh. ....................................................................................................................... 73 

Figure 4-22 Example of the peeling of a (a) layered beam where the adhesive is modeled as (b) linear 

until failure to demonstrate the joint element crack growth ability. .............................................. 74 

Figure 4-23. Load displacement plots for the peeling of a layered beam with different numbers of joint 

elements using (a) no re-meshing and (b) re-meshing. .................................................................. 75 
Figure 4-24. Load displacement plots for the joint depicted in Figure 4-14 with an elastic perfectly plastic 

adhesive with yield stress of 40 MPa. Plot shows (a) results of basing the constitutive properties 

on the peel to shear ratio of the highest stressed side and the lower side, elemental convergence 

for the (b) 2-D dense mesh model, (c) joint element, and (d) joint element with adaptive shape 

functions. ........................................................................................................................................ 76 

Figure 4-25. A typical DCB specimen, (a) the geometric parameters and boundary conditions for the DCB 

specimens and (b) the mesh for the joint element DCB model. ..................................................... 77 

Figure 4-26. The stress-strain relation of the adhesive was defined based on the Mode I energy release 

rate and modulus. ........................................................................................................................... 78 

Figure 4-27. Load vs. displacement curves for a) E719/IM7 and b) E7T1/G40 DCB specimens tested by 

Song and Waas [70] along with the present joint element model. ................................................. 78 

Figure 4-28. Using a) cylindrical tensile specimens and digital image correlation, the b) stress-strain 

relations of the adhesives EA 9394 and EA 9309.3NA could be defined by curve-fitting 

experimental data. .......................................................................................................................... 79 

Figure 4-29. Two adherends of an EA 9394 DCB specimen after complete failure. Adhesive found on 

both adherends indicates that failure occurred within the adhesive layer as desired. .................... 80 

Figure 4-30. Load vs. displacement curves for DCB specimens with aluminum adherends and EA 9394 

and EA 9309.3NA adhesive, along with the joint finite element model prediction. ..................... 81 

Figure 4-31. Geometric parameters for single lap joint tested by Harris and Adams [55]. ........................ 81 

Figure 4-32. (a) Stress-strain relation for bulk adhesive, along with peel and shear components for a single 

lap joint with ψ=1.4, and (b) corresponding load-displacement plots. .......................................... 82 

Figure 4-33. Peel to shear stress ratio in adhesive layer of the single lap joint as a function of end 

displacement. ................................................................................................................................. 83 

NASA/CR—2012-217606 vii 

NASA/CR—2012-217606 1

Chapter 1

Introduction

1.1 General Comments About This Thesis

Each chapter of this thesis contains work that has already being published as a conference or journal 

paper, and/or is being considered for journal publication. Because of this, each chapter has an introduction 

that motivates the particular aspects of the work relevant to that chapter. Therefore, in this introductory 

chapter, a general but brief overview that is relevant to the thesis topic is presented.

1.2 Introduction

Joining dissimilar materials is an important and contemporary subject receiving great attention from the 

aerospace engineering sector due to the large number of structural components that are being produced with 

a variety of materials, including metals, polymers, fiber reinforced composites (FRCs) and structural foams 

due to increasing demands on weight reduction. Traditional methods of joining, which have been largely 

based on bolted joint technology, are not very suitable for joining FRCs to each other or FRCs to other 

materials [1].

As bonded joints increase in popularity and use, the demand for modeling techniques increases also. In 

the past, analytical models have been favored as the preferred method of predicting stresses and strength, but 

finite element (FE) methods have emerged as the new standard in preliminary design due to necessity of 

analyzing and designing components that contain multiple joints where analytical techniques become 

intractable. FE based methods have been proven to be extremely powerful, but the small scale of the 

adhesive thickness when compared to the dimensions of the surrounding structure has kept joint FE analysis

largely out of global vehicle models. A fine mesh is needed to correctly model the adhesive layer producing 

an incompatibility in simultaneously analyzing the joint stresses accurately in conjunction with a very coarse 

model of an entire vehicle. Therefore, the actual design and sizing of joints is often put off until a later time, 

when small sub-models are used to look into the details of a vehicle.

To further complicate things, the eccentricity of the load path in most aerospace structural joints (such 

as the popular single lap joint [2] featured in Figure 2-5) causes significant bending when axially loaded.

This bending-introduced rotation makes nonlinear geometric effects significant, even in early stages of 

loading. Furthermore, most modern advanced polymeric adhesives show considerable nonlinear material 

behavior, which causes the joint to remain intact even when the yield stress of the adhesive is reached.

Therefore, it is imperative that that these factors be included in models in order to fully utilize the superior 

capability of bonded joints.

Motivated by these reasons, the overarching objective of this thesis was to develop a single finite 

element that can capture accurately the stress and strain states of a bonded joint while still facilitating its 

merger with the surrounding structure without having an incompatibility in finite element mesh densities.

The bonded joint element (simply referred to as "joint element") uses analytical structural models to find 

shape functions for a joint, allowing the joint region to be modeled with one element. That way, with a 

single element modeling the joint, it can be inserted into a larger finite element model efficiently.

Furthermore, large rotations, nonlinear material properties, and crack growth capabilities are included and 

analysis guidelines are developed to keep the number of elements required to a minimum.

1.3 Thesis Organization

The chapters are based on previously published papers, and can be read as stand-alone pieces of work.

For this reason, the introduction and motivation for each chapter might seem repetitive. However, in the 

actual body of the chapter, many aspects have been removed from their previous state to shorten up the 

work and avoid repetition, especially in the formulation. The notation has been largely unified, with one 
NASA/CR—2012-217606 2

difference being that in the earlier chapters, there is only one coordinate system. In the final chapter before 

the conclusion, the co-rotational formulation necessitates a local and global coordinate system. This local 

coordinate system correlates with the coordinate system in the earlier two chapters.

1.4 Significant Contributions

This section provides an executive summary of the significant contributions of this dissertation to the 

scientific community.

A linear elastic joint finite element used to find the stresses, strains, and displacements in an 

adhesively boned joint.

Method of modeling complex, modern joints with tapers, ply drops and angle changes using 

discrete joint element “building blocks.” 

A joint element that can model the performance of a joint with a functionally graded adhesive with 

just one element.

Demonstrated the stress reduction, sensitivity, and universal applicability of functionally graded 

adhesives.

A technique for applying a co-rotational formulation to layered beam models to capture large 

rotation effects.

A technique of internally re-meshing an element to introduce and grow a crack within the element.

A technique of adapting the shape functions in an element as the loading increases to account for 

material softening and improve elemental mesh convergence.

A method of using bulk adhesive tensile data in an adhesive spring model with uncoupled shear 

and peel responses.

1.5 Publications

The following related publications were available at the time of the dissertation defense:

[1] S.E. Stapleton, A.M. Waas, S.M. Arnold, Functionally Graded Adhesives for Composite Joints, 

International Journal of Adhesion and Adhesives. 35 (2012) 36-49.

[2] S.E. Stapleton, A.M. Waas, B.A. Bednarcyk, Modeling Progressive Failure of Bonded Joints Using 

a Single Joint Finite Element, AIAA Journal. 49 (2011) 1740-9.

[3] Ahn J, Stapleton SE, and Waas AM. Advanced Modeling of the Behavior of Bonded Composite 

Joints in Aerospace Applications. In: P.P. Camanho, L. Tong, Composite joints and connections: 

Principles, modelling and testing, Woodhead Publishing, 2011.

[4] S.E. Stapleton, A.M. Waas, B.A. Bednarcyk, Bonded Joint Elements for Structural Modeling and 

Failure Prediction, in: Proceedings of the 52nd AIAA/ASME/ASCE/AHS/ASC SDM Conference, 

Denver, CO, 2011.

[5] S.E. Stapleton, A.M. Waas, Reduced-order Modeling Of Adhesively Bonded Joints Using An 

Enhanced Joint Finite Element, in: Proceedings of the 52nd International SAMPE Symposium, Salt 

Lake City, UT, 2010.

[6] S. Stapleton, A. Waas, B. Bednarcyk, Modeling Progressive Failure of Bonded Joints Using a 

Single Joint Finite Element, in: Proceedings of the 18th AIAA/ASME/AHS Adaptive Structures 

Conference, Orlando, Florida, 2010.

[7] S.E. Stapleton, A. Waas, Macroscopic Finite Element for a Single Lap Joint, in: Proceedings of the 

AIAA/ASME/ASCE/AHS/ASC 50th SDM Conference, Palm Springs, California, 2009.
NASA/CR—2012-217606 3

Chapter 2

Linear Elastic Modeling of Adhesively Bonded Joints Using an

Enhanced Joint Finite Element

The design and sizing of adhesively bonded joints has always been a major bottleneck in the design of 

composite vehicles. Dense meshes are required to capture the full behavior of a joint, but these dense 

meshes are impractical in vehicle-scale models where a course mesh is more desirable to make quick 

assessments and comparisons of different joint geometries. Analytical models are often helpful in sizing, but 

difficulties arise in coupling these models with full-vehicle finite element (FE) models. Therefore, a reduced 

order joint finite element was created that can be used within structural FE models to make quick 

assessments of bonded composite joints. The shape functions of the joint finite element were found by 

solving the governing equations for a joint where the adherends were modeled as beams connected by 

adhesive modeled using various assumptions. By analytically determining the shape functions of the joint 

element, the complex joint behavior can be captured with very few elements. Analyses of joint stresses for 

different joints using the enhanced joint finite elements were found to agree well with analyses using 

standard, 2-D plane stress elements.

2.1 Introduction

With the increasing demand for composites in lightweight aerospace structures, adhesively bonded 

joints are becoming more critical than ever. Bolts and rivets introduce holes which cause significant stress 

concentrations and premature failure in composite materials, while adhesives spread the load more evenly 

over the composite while facilitating a lighter overall structure.

Traditionally, analytical models have been used to assess the performance of joints [2–7]. However, FE 

modeling has emerged as a popular and robust method for structural analysis. In order to properly predict 

the stresses in a joint using FE analysis, a dense mesh must be used, particularly near the adhesive ends.

While such a model is ideal for detailed analysis, it can be crippling when it comes to joint design and sizing

[8,9]. Furthermore, such a model does not couple easily with coarse, vehicle-scale models used for sizing of 

the vehicle. Therefore, a need exists to develop predictive tools for bonded joints that can be seamlessly 

coupled with large scale structural analyses without adding major computational demands. Such tools can 

be used to make quick mesh-independent assessments of bonded composite joints. Furthermore, they fit in 

into the computational hierarchy of virtual testing of aircraft structures [10], an area that is getting increased 

attention in the aerospace industry with the aim of lowering design cycle and certification costs. 

A solution to this problem involves merging analytical models with finite elements. Simplified 

structural models can be used to obtain shape functions that are exact for the assumptions of the model.

These shape functions can be used to formulate stiffness matrix for the problem at hand. As long as the 

assumptions remain valid, such an element would give the exact solution regardless of the number of 

elements used. 

This method has been used to calculate an stiffness matrix for different beam on elastic foundation 

problems [11,12]. More recently, Waas and Gustafson [13] have created an element to capture the behavior 

of a double overlap joint subjected to mechanical and thermal loads. The current authors have continued this 

effort by creating an stiffness matrix for a single lap joint [14].

However, joints in application rarely resemble the simple joints that can be solved easily to obtain 

enhanced shape functions. Adherend tapers and steps do not conform well to a constant thickness adherend 

model. The current study lays out the formulation for the stiffness matrix for a joint with an arbitrary 

number of adherends. Additionally, an approach is introduced to model realistic, complex joints using 

simple joint elements as building blocks. This approach is demonstrated using two different joint types, and 

compared with dense mesh 2-D element models to show: 1) how well the model predicts stress near an 

adherend step and 2) how well a smooth taper can be modeled using stepped joint elements.
NASA/CR—2012-217606 4

Figure 2-1. Geometric parameters for overlap region of an adhesively bonded joint

with multiple bonded layers and the width of the joint in the y-direction is b.

2.2 Analytical Formulation

The enhanced joint finite element uses an analytical formulation to obtain the stiffness matrix for N

number of adherends joined by N-1 adhesive layers. The geometric parameters around the ith adhesive layer 

are shown in Figure 2-1. The adhesives and adherends were assumed to be linearly elastic, but not 

necessarily isotropic. It is also assumed that the adherend stiffness is much greater than that of the adhesive, 

which is typical of standard aerospace applications. The strains and rotations are considered small. The 

subscript i refers to values associated with adherend i, and ai refers to adhesive layer i. Each adhesive and 

adherend layer has its own z-coordinate starting at the centerline of the layer, marked as zai and zi

respectively. The variables ui(x) and wi(x) refer to the x and z-direction centerline displacements of

adherend i.

The first section below shows different models considered for the adherends, namely an Euler-Bernoulli 

and a Timoshenko model. The next section lays out various models utilized to represent the adhesive layers.

Using various combinations of adherend and adhesive models, a system of governing equations can be 

formulated. Regardless of the models chosen, this system of governing equations can be solved following 

the method outlined in the third section, and the stiffness matrix can be obtained. Subsequent discussion will 

consider the different models and why or under what circumstances their use is favorable.

2.2.1 Adherend Models

The models for the adherends utilize assumptions for the stress, strain, and displacement in order to 

simplify the problem to make the possibility of obtaining an analytical solution within reach. For all models, 

the adherends are assumed to be linearly elastic, and layered transversely isotropic. It is assumed that the 

beam is very long in the y-direction, which puts the adherends in a state of cylindrical bending. The first 

model assumes the adherends to behave as an Euler-Bernoulli beam (with the assumption of plane strain 

rather than plane stress), while the second model includes shear deformations. Rather than presenting a 

detailed derivation for the three beam theories, three vectors/matrices needed to use these models will be 

defined for each model: k Di , Gi , and ui . The matrix k Di is defined with the relation

k k σ =D ε i ii

2-1

where

adhesive i

adherend i+1

l

adherend i

ηai

ti

ti+1

x

zi

zi+1

zai

ui(x), wi(x)

ui+1(x), wi+1(x)
NASA/CR—2012-217606 5

T k kk

i ii é ù s t

ë û σ =

2-2

and

[ ]

T

i ii ε = e g .

2-3

The normal stress and strain are in the x-direction of the kth layer of the ith adherend, while the shear stress 

and strain are in the xz plane. Similarly, the matrix Gi relates the strain to the centerline displacements with 

the equation

i ii ε =Gu 2-4

where ui , the centerline displacements of the ith adherend, will be defined for each model. Now k Di , Gi , 

and ui can be defined for each model, enabling a unified notation and simplifying a change between the 

beam theories.

2.2.1.1 Euler-Bernoulli Model

The Euler-Bernoulli beam model assumes that the axial stress due to extension and bending is the most 

dominant stress, and that all others can be ignored. The stress in layer k is related to the strain by k Di , which 

is given by:

11 0

0 0

k

k é ù Q

ê ú

ê ú ë û

D = i

.

2-5

where 11

k

Q is the 1,1 component of the transformed lamina stiffness matrix, 

k

Q , as defined by Classical 

Lamination Theory [11,12]. The adherend strain is related to the adherend centerline displacements by Gi , 

given as:

0100 0

0000 0 0

i

i

é ù z

ê ú

ë û

G =

.

2-6

Furthermore, the centerline displacement vector is defined as

, ,, , () () () () () ()

T

i i i x i i x i xx i xxx = é ù ux ux wx wx wx wx ë û u

2-7

where ,x denotes the derivative with respect to x. This form of defining the centerline displacements might 

not be conventional, but it is used to lead into our solution strategy of the governing equations. Using state 

variables with higher order derivatives as is done here allows the governing equations to be reduced to a 

series of first order differential equations.
NASA/CR—2012-217606 6

2.2.1.2 Shear Deformation Model

The shear deformation model assumes that plane sections remain plane, but not necessarily 

perpendicular to the centerline. The stress in layer k is related to the strain by k Di , which is given by:

11

66

0

0

k

k

k

Q

kQ

é ù

ê ú

ê ú

ë û

D = i

.

2-8

where the variable k is a correction factor used to offset the error caused by the simplifying assumption that 

the shear stress is uniform throughout a cross section of the beam. For a rectangular cross-section, the value 

of 5/6 is most commonly used [15]. The adherend strain is related to the adherend centerline displacements 

byGi , given as:

0100 0

0001 10

i

i

é ù z

ê ú

ë û

G =

.

2-9

Furthermore, a new variable is introduced to define the centerline displacements. The rotation due to 

bending only is defined by ( ) i x , making the centerline displacements:

,,, () () () () () ()

T

i i ix i ix i ix = é ù ux ux wx wx x x y y ë û u .

2-10

2.2.2 Adhesive Models

The adhesive models considered differ in the amount of detail included in the stress-strain relationship 

and the strain-displacement relationship. However, all models shown here will start out with the assumption 

that the displacements vary linearly in the z-direction. This allows the formulation to be strictly in terms of 

functions of x, as with the adherends. The adhesive is assumed to be perfectly bonded to the adherends. The 

difference in the models will be defined by two matrices: Dai , and Gai . Similar to the adherends, the 

matrix Dai relates the adhesive stresses to strains through

σ =D ε ai ai ai

2-11

where

[ ]

T

σ = ai ai ai xai sts 2-12

and

[ ]

T

ai ai ai xai ε = ege .

2-13

The normal stress and strain, ai and ai , are the stresses and strain in the z-direction, also commonly 

referred to as the peel stress and strain. The shear stress and strain, ai and ai , are in the xz plane, while the 
NASA/CR—2012-217606 7

axial stress and strain, xai and xai occur in the x-direction. Although this might not be the conventional 

way of ordering stresses and strains in elasticity, this is the chosen scheme that adheres closer to conventions 

when dealing with adhesive layers. The peel and shear stresses are the components of most concern, and the 

axial stress and strain are often ignored [3,4,6,5]. Therefore, it is placed at the end to allow for easy deletion 

if desired.

The matrix ˆ

Gai relates the strain of the adhesive to the adhesive displacements with the equation

ˆ ai ai ai ε =G u 2-14

where the adhesive displacements and derivatives,uai , are defined as 

,

,

(, )

(, )

(, )

(, )

ai ai

ai ai x

ai

ai ai

ai ai x

u xz

u xz

w xz

w xz

é ù

ê ú

ê ú

ê ú

ê ú

ë û

u =

.

2-15

The adhesive displacements are related to the adherend centerline displacements of the adherends above and 

below the adhesive layer through the equation

1

ˆ i

ai ai

i+

é ù

ê ú

ë û

u

u =H

u

2-16

where 

1 2

1 2

1 2

1 2

0 0 0 00 0 0 0 00

0 0 0 00 0 0 0 00 ˆ

0 0 0 00 0 0 0 00

0 0 0 00 0 0 0 00

ai

C C

C C

C C

C C

é ù

ê ú

=

ë û

H

2-17

and

1

1

2

ai

ai

z C = +

, 

2

1

2

ai

ai

z C = -

.

2-18

This relation was found by assuming that the displacements in the adhesive layer vary linearly in the xdirection, and that the adhesive and adherends are bonded perfectly. Using the above relations, the adhesive 

strain can be related to the adherend centerline displacements by

1

i

ai ai

i+

é ù

ê ú

ë û

u

ε =G

u

2-19

where
NASA/CR—2012-217606 8

ˆ ˆ G =G H ai ai ai .

2-20

Gai will be defined for each adhesive model in the following sections.

2.2.2.1 Model 1: Winkler Foundation

In this model, it is assumed that the adhesive is a bed of uncoupled linear shear and normal springs. The 

adhesive layer is assumed to be so thin that the stress in the adhesive layer is independent of the zcoordinate. The stress is related to the strain by the matrix 

0 0

0 0

0 00

ai

ai ai

E

G

é ù

ê ú

ê ú

ê ú

ë û

D =

2-21

where Eai and Gai are the normal and shear moduli of the ith adhesive layer, respectively. As can be seen, 

the axial stress is neglected which is a common assumption for adhesive layer models [4–6,16]. The normal 

strain is defined in the standard small strain elasticity manner, but the shear strain is simplified to

(, ) ai ai

ai

ai

u xz

z = 2-22

which makes it a true Winkler foundation, or bed of uncoupled linear shear springs. From this 

simplification, the matrix relating the adhesive strains with the adherend centerline displacements, Gai is 

given as

1 1

1 1 1

2 2

0 0 0 00 0 0 0 00

0 0 00 0 0 00

0 0 0 0 00 0 0 0 0 00

ai ai

t t i i

ai ai ai ai ai

h h

h hh h

+

é ù

ê ú

ë û

G =

.

2-23

This model is the simplest model, and will be utilized heavily in the study. The main reason for this is not 

necessarily its simplicity, but lies in the fact that it resembles a Cohesive Zone Model (CZM), with peel and 

shear components uncoupled. Although a significant simplification, the model is “tuned” using experimental 

tests to characterize the peel and shear components (Mode I and Mode II) and is capable of predicting 

mixed-mode failure events.

2.2.2.2 Model 2: Coupled Springs

The second model resembles the first model, except in the strain-displacement relation. The full shear 

strain equation is used:

(, ) (, ) ai ai ai ai

ai

ai ai

u xz w xz

z x

¶ ¶ + ¶ ¶ =

.

2-24
NASA/CR—2012-217606 9

This is the equivalent of having a non-Winkler foundation as the adhesive, where the normal springs are 

“tied together” and add to the shear rigidity. This model makes the adhesive shear strain/stress a function of 

z, as can be seen in

1 1

1 1 1

2 2 1 2

0 0 0 00 0 0 0 00

0 0 00 0 0 00

0 0 0 0 00 0 0 0 0 00

ai ai

t t i i

ai ai ai ai ai

C C

h h

hh h h

+

é ù

ê ú + +

ë û

G =

.

2-25

2.2.2.3 Model 3: Plane Stress

Model 3 uses the full stress-strain and strain/displacement relations for a body in plane stress. This 

model is appropriate when the depth of the joint (b) is much smaller than the length of the adherends. The 

stress for the plane stress problem is related to the strain by the matrix 

3 3

3 3

0

0 0

0

ai

ai ai

ai

C C

G

C C

é ù

ê ú

ë û

D =

.

2-26

where ai is the Poisson’s Ratio of the ith adhesive layer, and 

( ) 3 2 1

ai

ai

E C =

.

2-27

The adhesive strain is related to the adherend centerline displacements by Gi , given as:

1 1

1 1 1

2 2 1 2

1

1 2 2 1 2 2

0 0 0 0 00 0 0 0 0

00 00 00 0 0

0 0 0 00 0 0 0

ai ai

t t i i

ai ai ai ai ai

t t i i

C C

C CC C

h h

hh h h

+

+

é ù

ê ú + +

ë û

G =

.

2-28

2.2.2.4 Model 4: Plane Strain

Model 4 is similar to model 3, except that the stress-strain relation is for a body in plane strain.

Although this model was not used in the results, it is important for modeling the joining of wide panels. The 

stress for the plane strain problem is related to the strain by the matrix 

4 4

4 4

(1 ) 0

0 0

0 (1 )

ai ai

ai ai

ai ai

C C

G

C C

n n

n n

é ù

ê ú

ë û

D =

.

2-29

where ai is the Poisson’s Ratio of the ith adhesive layer, and 
NASA/CR—2012-217606 10

( )( ) 4 1 12

ai

ai ai

E C n n = + - .

2-30

2.2.3 General Formulation

Regardless of the adherend and adhesive models chosen, the strain energy of the joint, UJoint, can be

written as:

1

0

1 1 1

1

2

N M N l Tk T

joint k i i ai ai A Aai i i k i

U dA dA dx = = =

é ù = + ê ú

ë û òò ò åå å σε σ ε

2-31

where k Ai is the area (in the yz plane) of the kth layer of the ith adherend, Aai is the area of the ith adhesive 

layer, and M is the number of transversely isotropic layers of the ith adherend, and all other variables are 

previously defined. Using the principle of stationarity of potential energy, 2N fully coupled governing 

equilibrium differential equations are obtained from the energy expression. Of the 2N governing equations, 

N equations correspond to the axial equilibrium, while N equations correspond to the transverse equilibrium.

The axial displacement equilibrium equations contain second order derivatives, while the transverse 

displacement equations have fourth order derivatives. The order of these equations can be reduced and 

assembled into a system of first order constant coefficient homogeneous ordinary differential equations of 

the form

u Au ,x = .

2-32

where 

1

T TT T

i N = é ù

ë û uu u u K K 2-33

or a vector containing the centerline displacements of all N adherends in the overlap region.

Inspecting the matrix A can be helpful in determining the nature of the solution and determining the 

solution method. There are 6N eigenvalues of A : N real eigenvalues, 2N complex eigenvalues, and 3N

repeating eigenvalues. Therefore, the solution is made up of N exponential terms, 2N exponential terms 

multiplied by a sine or cosine, and the 3N repeating eigenvalues correspond to a third order polynomial 

found in a standard beam solution. Such a complex solution shows that merely employing standard beam 

shape functions to the joint problem would be inadequate in capturing the nature of the whole solution.

The system in can be solved using various methods, but calculating the matrix exponential was the 

chosen method because numerical boundary conditions are not required to obtain a solution. The solution of 

the system can be written in terms of the matrix exponential, Ax e , and a vector of unknown constants, C , as

x = A u eC .

2-34

The matrix exponential can be expressed as the infinite series [17]

0 !

k

x k

k

x

= k = Ae A .

2-35
NASA/CR—2012-217606 11

In order to obtain faster convergence, a method of scaling and squaring [18] was employed, and the 

series was calculated up to a value of k which yields an acceptable error, . The error can be defined many 

ways, but the current study defined the error as the difference between the 1-norms of Ax e for k-1 and k. The 

value of the acceptable error was set at = 0.0001.

The next step was to solve for the vector of constants, C , using the boundary conditions. This is where 

the analytical formulation is discretized, and the displacements are obtained in terms of the nodal

displacements as defined in Figure 2-2. For adherend i, the boundary conditions on the left side of the joint 

(x=0) can be expressed in the following equation:

(0) il q bu = c i

. 2-36

where (0) ui is ui evaluated at x=0, il q is a vector containing the prescribed nodal degrees of freedom of 

adherend i at x=0 (or the left side), and

100000

001000

000100

c

é ù

ê ú =

ë û

b .

2-37

Equation 2-36 for all N adherends can be assembled together, and a relation between the nodal degrees of 

freedom at x=0 and the vector of constants can be found using:

0

0

l = NA q Be C 2-38

Figure 2-2. Boundary conditions for adherend i: prescribed nodal displacements and rotations at x=0 and x=l.

adherendi 1

il

q

2

il

q

3

il

q 1

ir

q

2

ir

q

3

ir

q
NASA/CR—2012-217606 12

where

1l

l il

Nl

é ù

ê ú

ê ú

= ê ú

ê ú

ê ú

ê ú

ë û

q

q q

q

M

M .

2-39

and

c

N

c

é ù

ê ú = ê ú

ê ú

ë û

b

B

b

O .

2-40

where the subscript N denotes the number ofmatrices on the diagonal. After performing the same operations 

at x=l, all of the boundary conditions can be combined in the form

0

l

r

ì ü = = í ý

î þ

q zC q

q .

2-41

where

0

2N l

ì ü = í ý

î þ

A

A

e z B

e .

2-42

Using this relation, one can obtain an expression for the vector of unknown constants:

1 C zq 0 = . 2-43

This relation can be inserted into Equation 2-34 to get the adherend centerline displacements in terms of the 

nodal degrees of freedom,

u Nq = 2-44

where the shape functions, N are defined as:

x 1 = A N ez . 2-45

Next, Equation 2-31 is rewritten in terms of the centerline displacements using the stress-strain and straindisplacement relations and put into matrix form:
NASA/CR—2012-217606 13

Figure 2-3. The building block approach facilitates modeling complex joints with

simple joint element building blocks.

0

l T U joint = u Uudx . 2-46

The usefulness of this form is that the integral has been performed analytically over the area since the 

variation in the z-direction and y-direction is known, and numerical integration must only be performed in 

the x-direction. Inserting Equation 2-24 into Equation 2-46 and minimizing the energy yields the stiffness 

matrix of the joint:

0

l T

joint k N UN = dx . 2-47

This formulation gives the stiffness matrix for a simple region of constant thickness adherend overlap.

However, many joints in application contain complicated geometries, including ply steps and tapers. To use 

the simple joint element for efficient modeling of complex joints, a building block approach was 

implemented. This approach involves combining simple, constant-thickness joint sections to create 

complicated joints with very few elements. This concept is illustrated in Figure 2-3, where single, double 

and triple adherend joint building blocks are combined to make complicated joints such as a pi joint, tapered 

single lap joint, and spliced sandwich joint.

Adherends joined together at the adherend centerlines are related to each other with the equation

iL iR q q = , i =1..3

2-48

where the subscript L is for the adherend on the left, R is for the adherend on the right, and the numerical 

subscripts 1, 2, and 3 refer to axial, transverse, and rotational degrees of freedom respectively. To model a 

ply step or taper, the transverse and rotational degrees of freedom are equal, but the axial degree of freedom 

of the left adherend is related to that of the right adherend through the following equation: 

11 3 L R offset R q q tq = - 2-49

where offset t is the vertical (z-direction) offset distance between the two nodes (Figure 2-4). A discussion of 

the accuracy of this approximation will be discussed subsequently.

Beam

Single Lap

Double Lap
NASA/CR—2012-217606 14

Figure 2-4. Joints connected together with an offset between centerlines are related using Eq. (21).

2.3 Validation and Results

2.3.1 Adhesive Model Comparison

Two studies were done to validate the joint element and discover the simplest adhesive model required 

to accurately predict the shear and peel stress in a single lap joint. First, the adhesive models were compared 

with each other over a broad range of joint geometric parameters to show which models have an impact on 

the predicted stress for different parameter values. Second, a 2-D solution, based on the finite element 

method (FEM) was generated for four parameter cases and the predicted adhesive stresses were compared 

with the three models. This was done to illustrate the accuracy of the joint element for different geometric 

parameters.

The three adhesive models were compared over a range of parameters to determine the difference in 

peak shear and peel stress along the adhesive centerline (za1=0) predicted by the models as a function of the 

parameters. This is useful to show when assumptions about the adhesive stresses are valid and when one 

should be careful using them. Figure 2-5 shows the dimensions and material properties of the single lap joint 

considered. The adherends were aluminum, with a Young’s modulus of 70 GPa (E) and Poisson’s ratio of 

0.33 (υ1 and υ2). They were both 5 mm thick (t), 2 mm wide (b), and extended 5000 mm long past the joint 

overlap (ladherend). The adhesive was FM300 and had a Young’s modulus of 2.17 GPa (Ea) with a shear 

modulus of 0.89 GPa (Ga). The adhesive thickness (η) was varied from 0.005 to 5 mm, and the overlap 

length (l) was varied from 50 to 5000 mm. Although these parameters are not necessarily typical for a joint, 

it was necessary to test a wide range of parameters to demonstrate the working range of the joint element 

and the limits of the assumptions.

On either side of the joint element, 50 beam elements were used to model the adherends outside of the 

overlap region. The left end was clamped, or restrained from displacement. The right end was extended in 

the x-direction by 10 mm (Δ) and restrained from rotating and from displacing in the z-direction as shown in 

Figure 2-5.

Along with comparing the three models to each other, they were also compared with the results from a 

2-D FEM created using the commercial package Abaqus [19] for four different geometric parameter cases.

The same geometric parameters and boundary conditions shown in Figure 2-5 were used for the cases, and 

Table 2-1 shows the values of t/l and η/t used for each case. Cases were chosen at four extreme corners of 

the parameters tested for the comparison study.

q1R

q3R

q2L

toffset

q2R

q3L

q1L

q3Rtoffset
NASA/CR—2012-217606 15

Figure 2-5. Boundary conditions, loading scenario and geometric parameters for the single lap joint.

Material properties include: E1=E2= E, Ea, and Ga.

Table 2-1. Parameters of in-depth study cases.

t/l η/t

Case 1 0.001 1

Case 2 0.001 0.1

Case 3 0.1 1

Case 4 0.1 0.1

The finite element models were constructed using 70,000 - 100,000 2-D plane stress quadrilateral 

elements, with an element bias towards the corners of the adhesive. One issue that deserves mention is that 

the joint element models do not fulfill the traction free (σa and τa ) boundary condition at the free edges of 

the adhesive boundary. Therefore, the FEMs and the joint element models are not expected to predict 

similar stresses at the free edges of the adhesive. Moreover, the inside corners of the adhesive in the FEM 

causes a stress singularity, making the model mesh dependant in the corner singularity region. In 

application, spew fillets are present, which eliminate this stress concentration and cause the edges of the 

overlap to not have zero stress [20]. Therefore it is of no major concern that the joint element does not 

reflect the stresses predicted by the FEA model at the ends of the adhesive. Since the model is meant to 

serve as an initial vehicle-scale model element for initial sizing and not necessarily for detailed analysis, the 

goal of the validation is to show that the overall behavior of the joint is reflected by the joint element.

The plots comparing the relative difference between the maximum centerline peel and shear stresses 

predicted by Models 3 and 2 for different η/t and t/l values can be found in Figure 2-6. The relative 

difference between the maximum peel stress predicted by Models 3 and 2 is between 8% and 5% for the 

whole range of adhesive thicknesses and overlap lengths considered. The difference is greatest for the very 

thin adhesive layers, and least for the thickest adhesive layers. Additionally, the difference between the 

maximum adhesive centerline shear stress predicted by the two models is below 1% for all of the parameters 

considered in this study. Considering the uncertainty of the stress state at the corner of the adhesive due to 

the previously mentioned issues related to stress concentrations and violation of the traction free condition, 

the difference between these models is very small. If one is only concerned with the peel and shear stress in 

the adhesive layer, it would be advantageous to use Model 2 over Model 3 due to its increased simplicity.

However, Tsai and Morton [21] note that the magnitude of σxa near the end of the overlap is often 

comparable to the magnitude of the peel and shear stresses. Therefore, it should be kept in mind that Model 

3 is the only one of the three models to yield the extensional stress, σaxx.

Figure 2-7 contains a comparison of the relative difference between the maximum centerline peel and 

shear stresses predicted by Models 2 and 1 for different η/t and t/l values. It can be seen that these models 

predict very different values of maximum stress. The shear stress difference can reach up to 40%, while the 

peel stress difference can be almost 30%. It appears that the differences between the two models are 

especially large for thicker adhesive layers and shorter overlaps. Interestingly, the sign change between t/l

values of .001 and .01 in both plots indicates that there may be a t/l value which causes the difference to 

hover around zero.

l ladherend

t Δ

b

ladherend

η
NASA/CR—2012-217606 16

Figure 2-6. Relative difference in maximum adhesive centerline (za1=0) shear and peel stress between

Models 3 and 2.

Figure 2-7. Comparison of the maximum shear and peel stress predicted by Models 2 and 1 for different

joint configurations, along with points indicating the cases for comparison with 2-D FEMs.

2-D plane stress FEMs were constructed for the four cases, each case with parameters as described in

Table 2-1 and the points labeled in Figure 2-7. The results of the comparison between the FEM models and 

the joint element models are found in Figure 2-8 up to Figure 2-16. Figure 2-8 compares the reaction force 

predicted by the three adhesive models, normalized and compared with that of the 2-D element dense mesh 

model. Since the joints were given a prescribed displacement, this plot of the reaction force also represents 

how well the stiffness was captured by the joint models. This corresponds with the shear stress comparison 

discussed earlier, since the reaction force for a linearly elastic joint can be equated to the integral of the 

shear stress over the are in the area of the adhesive in the xy-plane. Therefore, model 1 deviated from the 

other models and the Abaqus model the most for cases 1 and 3, where the adhesive was the thickest. An 

investigation of the adhesive centerline stresses along with 2-D contour plots showing the shear and normal 

stress distribution in the adhesive as predicted by the FEM and three models for each case will be presented 

in the following section for each case for a more detailed discussion.

2.3.1.1 Case 1

Case 1 had a thick, long adhesive layer with η/t = 1 and t/l = 0.001. The length of the adhesive displayed 

in Figure 2-9 and Figure 2-10 is only 0.5% of the length because the adhesive in the middle region is 

effectively stress free, which makes it this region unimportant for the current study. Obviously, this joint 

would be very inefficient because almost all of the stress is held by less than 1% of the adhesive. All three 

models predict the peel stress fairly well, although Models 3 and 2 appeared to be slightly more accurate. 

Looking at the adhesive peel stress distribution in Figure 2-10, Models 2 and 1 predict constant peel stress in 

the z-direction, while the peel stress distribution of the 2-D FEM varies nonlinearly in the z-direction, even 

past the free end.

0.0

0.2

0.4

0.6

0.8

1.0

0 0.2 0.4 0.6 0.8 1

η/t

max τa difference (%)

-8.5

-8.0

-7.5

-7.0

-6.5

-6.0

-5.5

-5.0

0 0.2 0.4 0.6 0.8 1

η/t

max σa difference (%)

-10

0

10

20

30

40

0 0.2 0.4 0.6 0.8 1

η/t

max τa difference (%)

-30

-20

-10

0

10

20

0 0.2 0.4 0.6 0.8 1

η/t

max 

σa difference (%)

Case 3

Case 1

Case 4

Case 2

Case 3

Case 1

Case 4

Case 2
NASA/CR—2012-217606 17

Figure 2-8. Reaction force predicted by the three adhesive models compared to the

2-D element dense mesh model for single lap joins with different geometries.

Figure 2-9. Adhesive centerline stresses predicted by different models for Case 1 parameters.

0

0.5

1

1.5

2

Case 1 Case 2 Case 3 Case 4

Abaqus 2-D

Model 1

Model 2

Model 3

F/FAbaqus

-2

0

2

4

6

8

0.495 0.496 0.497 0.498 0.499 0.500

 

0

1

2

3

4

5

0.495 0.496 0.497 0.498 0.499 0.500

 

Abaqus 2-D elements

Model 1

Model 2

Model 3

 

 

 

 

 

 

τa (MPa)

σa (MPa)

x/l x/l
NASA/CR—2012-217606 18

Figure 2-10. Adhesive stress contours predicted by different models for Case 1 parameters.

Looking at the centerline shear stress, Model 1 does not seem to match up with the FEM very well. 

While Models 3 and 2 are low, they appear to capture the general trend. Obviously, none of the models 

reflect the free end condition of zero traction. The shear stress contour plots in Figure 2-10 show that Model 

1 has shear stress constant in the z-direction, but the FEM model shows that the stress varies significantly in 

the z-direction, which might explain the inaccuracy of the model. The distribution of shear stress in Models 

3 and 2 is a closer match to the FEM model, although the slopes of the lines of constant stress appear to be 

incorrect past the free end.

2.3.1.2 Case 2

Case 2 was composed of a joint with a thin, long adhesive layer (η/t = 0.1 and t/l = 0.001). The length of 

the adhesive displayed in Figure 2-11 and Figure 2-12 is only 0.5% of the length because only the ends 

display behavior worth comparison. All three models predict the peel and shear stress very well. The reason 

can be seen in the contour plots in Figure 2-12. Although it is difficult to see, the effect of the stress 

singularity and traction free boundary at the free end of the adhesive predicted by the 2-D FEM dies out 

very quickly, and the remainder of the adhesive distribution in the z-direction similar to that predicted by the 

joint element models. Since the thickness of the adhesive is so small, the stress can be effectively modeled 

as constant in the z-direction. Therefore, using the simplified Model 1 would still yield very accurate results 

for this case.

7

 

 

0

 

 

-2

 

Model 2

Model 3

Abaqus 2-D elements

Model 2

Model 3

Model 1

.495 0.496 0.497 0.498 0.499 0.500

 

 

Abaqus 2-D elements

Model 1

.495 0.496 0.497 0.498 0.499 0.500

 

 

7

τa (MPa)

σa (MPa)

x/l x/l
NASA/CR—2012-217606 19

Figure 2-11. Adhesive centerline stresses predicted by different models for Case 2 parameters.

Figure 2-12. Adhesive stress contours predicted by different models for Case 2 parameters.

2.3.1.3 Case 3

Case 3 had a thick, short adhesive layer with η/t = 1 and t/l = 0.1. The right half of the adhesive is 

shown in the plots. The three model predictions of the shear stress were not very consistent with the 2-D 

FEM. Model 1 vastly over predicts the stress levels, and Models 3 and 2 under predict the shear stress in 

most of the adhesive (Figure 2-13). Looking at the shear stress distribution of the FEM in Figure 2-14, the

effects of the free end boundary do not die out quickly as with the thin adhesive cases. The constant zdirection stress distribution of Model 1 is insufficient to correctly model the adhesive, especially when it 

comes to the shear stress.

The peel stress, on the other hand, was predicted reasonably well. All three models predicted similar 

trends for the adhesive centerline peel stress, even though the prediction appears to lag behind the stress 

predicted by the 2-D FEM. Generally, Model 1 would be a poor choice for a joint of this type. Models 2 and 

3 are closer, but still do not accurately predict the stress levels of the single lap joint for Case 3.

-5

0

5

10

15

20

0.495 0.496 0.497 0.498 0.499 0.500

 

0

3

6

9

12

15

0.495 0.496 0.497 0.498 0.499 0.500

 

Abaqus 2-D elements

Model 1

Model 2

Model 3

 

 

 

 

 

 

τa (MPa)

σa (MPa)

x/l x/l

Model 2

Model 3

Model 1

.495 0.496 0.497 0.498 0.499 0.500

 

 

Abaqus 2-D elements

.495 0.496 0.497 0.498 0.499 0.500

 

 

Model 2

Model 3

Model 1

Abaqus 2-D elements

0 τa (MPa) 15 -4 σ 17 a (MPa)

x/l x/l
NASA/CR—2012-217606 20

Figure 2-13. Adhesive centerline stresses predicted by different models for Case 3 parameters.

Figure 2-14. Adhesive stress contours predicted by different models for Case 3 parameters.

-2

0

2

4

6

8

10

0.0 0.1 0.2 0.3 0.4 0.5

 

 p p

0

1

2

3

4

5

6

7

0.0 0.1 0.2 0.3 0.4 0.5

 

Abaqus 2-D elements

Model 1

Model 2

Model 3

 

 

 

 

 

 

τa (MPa)

σa (MPa)

x/l x/l

7

 

 

0

 

 

8

 

-2

 

Model 2

Model 3

Abaqus 2-D elements

Model 2

Model 3

Model 1

Abaqus 2-D elements

Model 1

0.0 0.1 0.2 0.3 0.4 0.5

 

 

0.0 0.1 0.2 0.3 0.4 0.5

 

 

τa (MPa)

σa (MPa)

x/l x/l
NASA/CR—2012-217606 21

2.3.1.4 Case 4

Case 4 was composed of a joint with a thin, short adhesive layer (η/t = 0.1 and t/l = 0.1). All three 

models under predicted the shear stress slightly, and surprisingly, Model 1 is the closest to the 2-D FEM

(Figure 2-15). The contour plot of the FEM shear stress in Figure 2-16 shows that the effects of the free 

edge disappear quickly in the adhesive. The lines of constant shear stress appear linear, but Models 3 and 2 

have the wrong slopes, which probably accounts for the inaccuracy. Since Model 1 results in stresses 

constant in the z-direction, it is slightly more accurate than Models 3 and 2. The models predicted the FEM 

peel stress slightly more accurately, but there was still a lag similar to Cases 1 and 3. It is unclear what 

causes the lag; whether it is an effect created by the difference in free edge conditions, or simply due to 

inaccuracies in the model. 

The models presented here reflect the behavior of the joint in a global sense to be used for early sizing 

studies. When one keeps that in mind, the stress states predicted by the models for all of the cases appear to 

be adequate. The models are especially suited for thin adhesive layers, and appear to be even more accurate 

for long, thin adhesive layers. This tends to minimize edge and through-thickness effects, which are only 

approximated by the analytical models.

Figure 2-15. Adhesive centerline stresses predicted by different models for Case 4 parameters.

Figure 2-16. Adhesive stress contours predicted by different models for Case 4 parameters.

-10

-5

0

5

10

15

20

25

30

0.0 0.1 0.2 0.3 0.4 0.5

 

0

4

8

12

16

20

0.0 0.1 0.2 0.3 0.4 0.5

 

Abaqus 2-D elements

Model 1

Model 2

Model 3

 

 

 

 

 

 

τa (MPa)

σa (MPa)

x/l x/l

Model 2

Model 3

Model 1

Abaqus 2-D elements

Model 2

Model 3

Model 1

Abaqus 2-D elements

0 19 -6 26

0.0 0.1 0.2 0.3 0.4 0.5

 

 

0.0 0.1 0.2 0.3 0.4 0.5

 

 

τa (MPa) σa (MPa)

x/l x/l
NASA/CR—2012-217606 22

2.3.2 Composite Adherends

The composite adherend modeling capability of the joint element model was compared with the 

commercial sizing code HyperSizer [22], which uses an analytical formulation derived by Mortensen and 

Thompsen [6]. The two were compared for a bonded doubler under tension. The material properties are 

outlined in Table 2-2, the layup and boundary conditions can be found in Table 2-3, and the geometric 

parameters are shown in Figure 2-17. The joint element adherend and adhesive models chosen for the study 

were the Euler Bernoulli adherend model and adhesive model 1 (Winkler foundation), which are most 

consistent with the HyperSizer joint analysis method. Figure 2-18 shows a comparison of the normal 

stresses in the x-direction (σ1) through the thickness of the upper adherend at a distance of half of the length 

of the stiffener from the supported edge. As can be seen, the stresses in the layers are nearly identical.

Table 2-2. Material properties used in Hypersizer analysis [22].

E1 

(GPa)

E2 

(GPa)

E3 

(GPa)

G12 

(GPa)

G31 

(GPa)

G23 

(GPa) ν12 ν13 ν23

Boron/Epoxy 223 24.1 24.1 8.48 8.48 8.48 0.23 0.23 0.32

Epoxy 3.07 3.07 3.07 1.14 1.14 1.14 0.348 0.348 0.348

Table 2-3. Layup and boundary conditions for Hypersizer analysis [22].

Material Layup

Loads and boundary 

conditions

Adherend 1 Boron / Epoxy (±45/0/90/0/90/±45/0) Left 

1

1 3 0 l

q = 2

1,3 0 l

q =

Adherend 2 Boron / Epoxy (0/90/±45/90/0) Right F=1 N/mm

Adhesive Epoxy -

Figure 2-17. Geometric parameters for composite stiffened panel comparison with HyperSizer [22].

25.4 mm 30 mm

2.29 mm F

0.762 mm

0.102 mm

Adherend 1

Adherend 2
NASA/CR—2012-217606 23

Figure 2-18. Comparison of normal stresses in the x-direction of the mid-joint cross section

of adherend 1 of a stiffened panel found using HyperSizer [22] and the joint element.

2.3.3 Tapers and Steps

In order to check the validity of the model and its ability to accurately predict stresses in a joint with 

tapers and steps, two joints were modeled with the joint elements and a dense mesh model using the 

commercial software Abaqus with 2-D plane stress elements. The joint element adherend and adhesive 

models chosen for the study was the Euler Bernoulli adherend model and adhesive model 1 (Winkler 

foundation).The first joint modeled was a stepped double-strap joint, which was chosen specifically to 

demonstrate how well the joint element could predict stresses in the vicinity of the steps. Second, a tapered 

joint was modeled to explore how well a smooth taper could be modeled using discrete, stepped joint 

elements. Stresses obtained using both models were compared for both joints, and conclusions were drawn 

in order to guide joint designers in the correct use of the joint element.

2.3.3.1 Stepped Double-Strap Joint

A stepped double-strap joint was modeled using reduced-order joint elements and a dense mesh with

2-D plane strain elements. The geometric and material parameters used are found in Figure 2-19. The joint 

was modeled using 26 joint elements and around 99,000 2-D plane strain elements in Abaqus (Figure 2-20).

The left end was clamped, and the right end was restrained from rotation and transverse displacement. The 

stresses in the upper and lower adhesive layers are shown Figure 2-21a and b, while the stresses in the upper 

and lower surfaces of the left and right adherends are shown in Figure 2-21c and d.

Figure 2-19. Geometric and material parameters of the stepped double-strap joint.

-1.20

-0.80

-0.40

0.00

0.40

0.80

1.20

0 0.2 0.4 0.6 0.8 1

z1 (mm)

σ1 (MPa)

HyperSizer

JointElement

EAl=110 GPa, ETi=70 GPa, EEA 9394=4 GPa, and GEA 9394=1.79 Gpa

40 mm

50 mm 40 mm

Δ= 0.1 mm

5 mm

12 mm

40 mm

η = 0.5 mm

Aluminum Titanium

EA 9394

3 mm
NASA/CR—2012-217606 24

Figure 2-20. Two models compared in the study, (a) joint element model using 26 elements and (b) Abaqus 2-D 

plane stress model using 99,000 elements.

The shear and normal stresses in the centerline of the adhesive layers agreed very well between the two 

models, even in the regions near the steps. The only discrepancies were found near the edges of the 

adhesive. The 2-D Abaqus model treats the adhesive as a continuum, while the joint element treats the 

adhesive as normal and shear springs. Therefore, stress singularities and traction free boundary conditions 

are not captured in the joint element model. However, for long overlap regions and thin adhesive layers, the 

difference is minimal.

The adherend stresses are also in good agreement (Figure 2-21c and d) except at the reentrant corners 

and the loaded ends. As stated above, the stress singularity of the reentrant corners is not captured by the 

model. However, overall trends are captured (reaction force for the Abaqus 2-D model was 4.323 kN, while 

the joint element model predicted 4.326 kN), making this model a good tool for initial design and sizing, 

particularly considering that in reality, the spew fillet would mute the discrepancy at the reentrant corners.

It is the stress in the doublers which requires special attention. First, consider the surface of the doublers 

away from the surface of the steps (Figure 2-22a and d). The stress peaks in the joint model around the steps 

appear to correlate quite well with the stress predicted by the Abaqus 2-D model. Near the steps, the joint 

element model, which treats the adherends as beams, predicts stress oscillations. The FE model, which treats 

the adherends as continua, does not predict these oscillations at this surface, which is some distance from the 

steps.

The stress on the surface of the steps (Figure 2-22b and c) shows fairly good correlation between the 

models, except for the last step (the outer peaks in Figure 2-22b and C). It is believed that since the joint 

element model does not take the traction free surface of the step into account, the difference between the 

two models grows as a larger percentage of the thickness is a stress free surface. At the last step, half of the 

surface is traction free, which causes a discrepancy between the two models. However, despite the fact that 

the joint element model doesn’t capture stress singularities at reentrant corners and ignores traction free 

surfaces of the steps, the joint element model appears to correlate well with dense 2-D finite element 

meshes.

26 Elements

99,000 Elements

Joint Element

Abaqus 2-D

(a)

(b)
NASA/CR—2012-217606 25

Figure 2-21. Stresses in the (a) upper and (b) lower adhesive centerlines, and (c) left and

(d) right adherends in the double-strap joint.

x z

x z

τxz

σz

τxz

σz

x z x z

Joint Element Model Abaqus 2-D Model

(a) (b)

(c) (d)
NASA/CR—2012-217606 26

Figure 2-22. Stresses in the (a) lower and (b) upper surfaces of the top doubler, and (c) lower and (d) upper surface 

of the bottom doubler in the double-strap joint.

x z

x z

(a) (b)

(c) (d)
NASA/CR—2012-217606 27

2.3.3.2 Tapered Joint

A tapered joint was also modeled with both models to see how well the smooth taper could be 

approximated by a discrete stepped model using joint elements. The parameters for the joint are illustrated 

in Figure 2-23. The joint was modeled using 7, 12, and 22 joint elements (Figure 2-24a-c). This was to show 

the effect of mesh density in the vicinity of the joint. All three were compared with a 120,000 2-D plane 

stress element model in Abaqus (Figure 2-24d).

The stresses in the upper and lower surface of the left adherend near the taper are shown in

Figure 2-25a and b, while those of the right adherend can be found in Figure 2-25c and d. The normal and 

shear stresses in the adhesive layer are shown in Figure 2-25e and f, respectively. The stresses in the stepped 

joint model appear to oscillate around those of the smooth taper Abaqus 2-D model in the adherends

(Figure 2-25a-d). However, as the density of the joint elements increases, the stress predicted appears to 

converge on the stress in the Abaqus 2-D smooth taper model. 

Figure 2-23. Geometric and material parameters of the tapered joint.

Figure 2-24. Models compared in the study: joint element using (a) 7, (b) 12, (c) and 22 elements, and

(d) Abaqus 2-D plane stress model using 120,000 elements.

15 mm 50 mm

Δ= 0.1 mm

5 mm

12 mm

50 mm

η = 0.1 mm Aluminum EA 9394 Titanium

7 Elements

120,000 Elements

Joint Element

Abaqus 2-D

12 Elements

22 Elements

(a)

(b)

(c)

(d)
NASA/CR—2012-217606 28

Figure 2-25. Stresses in the (a) left and (b) right adherends and (c) peel and

(d) shear stress in the adhesive of the tapered joint.

-30

0

30

60

90

120

150

180

45 50 55 60 65

Abaqus 2-D

5 Steps

10 Steps

20 Steps

σx (MPa)

x (mm)

-30

0

30

60

90

120

150

180

45 50 55 60 65

Abaqus 2-D

5 Steps

10 Steps

20 Steps

σx (MPa)

x (mm)

x z

-50

0

50

100

150

200

250

0 10 20

Abaqus 2-D

5 Steps

10 Steps

20 Steps

σx (MPa)

x (mm)

-50

0

50

100

150

200

250

0 10 20

Abaqus 2-D

5 Steps

10 Steps

20 Steps

σx (MPa)

x (mm)

x z

-10

0

10

20

30

40

0 5 10 15

5 Steps

10 Steps

20 Steps

Abaqus 2-D

σz (MPa)

x (mm)

-10

0

10

20

30

40

0 5 10 15

Abaqus 2-D

5 Steps

10 Steps

20 Steps

τxz (MPa)

x (mm)

Abaqus 2-D

5 Steps

10 Steps

20 Steps

5 Steps 10 Steps 20 Steps

(a) (b)

(c) (d)

(e) (f)
NASA/CR—2012-217606 29

The peel stresses in the adhesive centerline appear to be quite different between the Abaqus 2-D model 

with a smooth taper and the stepped joint element models. Both have maximum stresses at the ends of the 

adhesive layer, but the joint element model shows a near zero stress away from the ends while the Abaqus 

2-D model has a nearly constant stress of about 8 MPa in the middle of the joint. This is due to the way the 

joint element represents the smooth taper. The adhesive layer is modeled as horizontal discrete pieces of 

adhesive only, while the vertical component is ignored. The vertical component would carry a near-constant 

peel stress, which is why the peel stress in the middle of the joint is not captured. Therefore, one must keep 

in mind that this component of stress will not be captured by the joint element model. However, the joint 

element model appears to converge to the Abaqus 2-D solution for the shear stress in the adhesive layer.

2.4 Conclusions

Eccentric load paths caused by single lap joints often require a detailed, dense mesh in structural finite 

element models of assemblies containing single-lap joints, which can be costly for global, vehicle scale 

models. To reduce computational time required for these large-scale models and aid joint sizing early on in 

the design phase, a joint element was created to model the behavior of a joint with a single element. This 

was accomplished by embedding an approximate analytical solution of the single lap joint into a finite 

element. The model discussed in this chapter assumes that all materials remain linearly elastic, but the 

extension of this model to include progressive failure discussed in Chapter 4 will require an iterative 

solution procedure. With this in mind, it is imperative that the simplest formulation be used to reduce the 

complexity and computing time required for each iteration. Therefore, three adhesive models of varying 

complexity were compared.

The first adhesive model, Model 1, treats the adhesive as a bed of uncoupled linear springs (Winkler 

foundations), making the stresses independent of the z-direction. Model 2 couples the normal springs 

together so that the variation of w displacements across the joint aids in shear strain accumulation. Finally, 

Model 3 assumes that the adhesive is in plane stress, and considers σaxx, σazz, and τaxz. Additionally, a model 

assuming plane strain was created, but not compared in this study. All models assumed the displacements in 

the adhesive layer to vary linearly in the z-direction. The maximum peel and shear stress predicted by these 

three models was compared to ascertain how much of a difference each additional assumption matters. It 

was found that there was little difference between Models 2 and 3. Since Model 2 is simpler than Model 3, it 

can be concluded that Model 3 is unnecessary, unless the value of σaxx is desired. On the other hand, Models 

1 and 2 predicted very different maximum stresses, especially for joints with thick adhesive layers and short 

joint overlap lengths.

Four joints with different combinations of thick and thin adhesive layers and long and short overlap 

lengths were studied in depth and compared with 2-D finite element models to 1) identify which models 

were more accurate, and under what parameter conditions and 2) to show that the joint finite element could 

be used to give a general prediction of joint behavior. It was found that all of the models were more accurate 

for joints with thin adhesive layers. Also, all of the models were slightly more accurate for longer joint 

overlaps. Therefore, for thin joints, Model 1 should be used because it can predict an answer just as accurate 

as the other models, with a much simpler formulation. For thick joints, Model 2 would be preferred because 

it is more accurate than Model 1, although the models are generally not as accurate for thick joints.

However, most joints in application have very thin adhesive layers, so Model 1 should be adequate for most 

real-life situations. It was also shown that for all cases, the joint elements were more than adequate at 

predicting the behavior of a joint for early design sizing purposes.

A second validating example compared the joint element with HyperSizer [22] sizing software for a 

stiffened composite panel. The comparison showed an almost perfect match of the stresses in the upper 

adherend layer.
NASA/CR—2012-217606 30

Finally, it was shown how well the joint element can model joints with tapers or steps. For discrete 

steps, the joint element appeared to agree with the 2-D dense mesh finite element solution. For a gradual 

taper, it was illustrated that the stresses in the adherends predicted by the joint element model appear to 

oscillate around the dense mesh solution, and that the stress converges to the solution by increasing the 

number of steps. Also, the lower bound of the adhesive shear stress is a good representation of the smooth 

taper solution. The average peel stress is not captured using the joint element model, but the maximum 

stress at the ends of the adhesive are still predicted, and these high stresses are of concern in predictions 

related to joint strength and durability.
NASA/CR—2012-217606 31

Chapter 3

Functionally Graded Adhesives for Composite Joints

Adhesives with functionally graded material properties are being considered for use in adhesively 

bonded joints to reduce the peel stress concentrations located near adherend discontinuities. Several 

practical concerns impede the actual use of such adhesives. These include increased manufacturing 

complications, alterations to the intended grading due to adhesive flow during manufacturing, and whether 

changing the loading conditions significantly impact the effectiveness of the grading. An analytical study is 

conducted to address these three concerns. An enhanced joint finite element, which uses an analytical 

formulation to obtain shape functions, is used to model the joint. Furthermore, proof–of-concept testing is 

conducted to show the potential advantages of functionally graded adhesives. In this study, grading is 

achieved by strategically placing glass beads within the adhesive layer at different densities along the joint.

3.1 Introduction

One major drawback of adhesively bonded joints is that the load path eccentricity causes the appearance 

of peel stress concentrations at the end of the adhesive layer. There has been a vast amount of research 

conducted in an attempt to reduce these stress concentrations, such as tapering the end of the adherend [23], 

increasing thickness of the adhesive at the end [24], fillets [20], novel joint geometries [25], and joint 

insertions [26], to name a few. All of these methods involve local details of adherend geometry (except the 

adhesive fillets), which typically increases part complexity.

Material grading occurs in nature at material interfaces to reduce stress concentrations [27]. Biological 

interfaces such as tendon to bone joints have been found to have graded material properties to distribute 

stress more evenly across the joint [28]. In this same spirit, material grading has been applied to 

adhesively bonded joints. Although grading the adherends has shown promise [29], many more researchers 

have investigated grading the adhesive properties. Some of the earliest grading of the adhesive was reported

by Patrick [30] and Raphael [31], where grading was discretely achieved using two adhesive materials (i.e., 

bi-adhesive). Recently many other researchers have investigated bi-adhesive joints with mixed results.

Sancaktar and Kumar [32] graded the adhesive by making rubber toughened regions, and found that 

the selectively toughened joints had the same strength as the fully toughened joints. Piresa et al. [33] used 

two adhesives to bond aluminum, and found up to a 22% increase in joint strength. Fitton and Broughton 

[34] bonded carbon reinforced plastic (CFRP) to steel, and found that it was crucial to optimize the amount 

of each adhesive and that some configurations did not benefit from grading. Da Silva and Lopes [35] looked 

at the influence of the ductility of adhesives on bi-adhesive joint strength experimentally and theoretically. It 

was found that the bi-adhesive joints out-performed joints with the more brittle adhesive alone, but did not 

always improve on the more compliant single adhesive joints. More recently, Kumar and Pandey [36]

modeled bi-adhesive joints using nonlinear 3D finite elements compared with a 2D finite element model, 

and found that the 2-D model could not fully capture the complex multi-axial stress state. Valleé et al. [37]

investigated bi-adhesive joints, among other stress reduction methods, and found that the adhesive stress 

was not linked to the joint strength of the configurations tested (which displayed adherend failure).

It appears that the first instance of grading the adhesive with a non-stepwise function was Kumar [38].

This purely theoretical investigation first compared a continuous (non step-wise) functionally graded 

adhesive (FGA), where the modulus was graded using a quadratic function, with a step-wise graded 

equivalent for different overlap lengths and adhesive thicknesses in a tubular joint. It was found that the 

continuous FGA reduced the shear and peel stresses in all cases. Second, four “arbitrarily chosen” functions 

were compared to show that the grading function can be manipulated to optimize joint performance. The 

current study aims to increase the understanding of both bi-adhesive and non-stepwise FGAs to make them 

a more viable, realistic, and advantageous choice for actual application in composite structures.
NASA/CR—2012-217606 32

Some potential drawbacks to FGAs have also been discussed in the literature. The first two were 

identified by Hart-Smith [24], where he pointed out that: 1) small gains over just using the ductile adhesive 

alone may be inadequate when considering the production difficulties and 2) during manufacturing, there is 

the “inevitable tendency for the stiff adhesive to squeeze out and displace the ductile adhesive,” making it 

probable that the resulting joint will be worse-off than using the ductile adhesive alone. A third concern was 

raised by Aboudi et al. [39] while investigating the response of metal matrix composites with tailored 

microstructures. They found that functionally grading the properties of a material may be detrimental when 

the loading is changed, such that the stress gradient in the material is reversed; this is a valid concern for 

practical situations where all loading cases are not necessarily known.

The current study aims to address these practical concerns to show that FGAs are a viable means of 

decreasing the peel stress in an adhesively bonded joint. An analytical model is constructed and used to 

compare the stresses in a butt-end joint configuration with four different functions of modulus graded 

adhesive: constant (single adhesive), discrete (bi-adhesive), linear, and exponential. First, the potential gains 

in stress reduction for FGA joints over joints with a single adhesive are shown. Along these lines, it is 

shown that additional stress reductions can be achieved by lowering the modulus of the more compliant 

adhesive. Since it is likely that step-wise grading will appeal from a manufacturing perspective, stress 

reduction of a step-wise graded adhesive with many steps is investigated (with single or bi-adhesive being a 

special case). Second, the study addresses the issue of adhesive flow during bonding by showing how 

sensitive the optimum for the three FGAs is to perturbation of the grading. Third, multiple load scenarios are 

examined to address the concern of changing loading conditions. Results indicate that the stress magnitude 

gradient in the adhesive remains unchanged under different loading conditions, making joints a perfect 

application for material grading. Additionally, it is shown that the grading can still be optimal under 

different loading cases. Addressing these three concerns provides significant impetus for the use of FGAs in 

industrial applications.

The model used to obtain the adhesive stress for different FGAs is a structural finite element made 

specifically for adhesively bonded joints. Motivated by the desire to create a computationally efficient tool 

for designing joints within a coarse, vehicle scale finite element model [13], we combined an analytical 

formulation with a finite element in the joint element. This concept has been often referred to as the exact 

stiffness matrix method, and has been previously applied to the beam on an elastic foundation problem 

[11,12]. The joint element is capable of capturing the stresses in a mesh-independent, efficient manner. Such 

an efficient method is pivotal to the current study, allowing the analysis of over 20,000 different joints for 

optimization and parametric studies on a desktop computer in a fraction of a second per joint. The 

formulation adopted here is altered from the previous section to account for a graded adhesive modulus and 

is presented below. A linear elastic material model is used for several reasons, simplicity being the most 

prominent. Also, since a controlled method and material system for manufacturing FGAs has not yet been 

identified, failure of the joint and post-failure response is not addressed. Finally, it should be noted that after 

initial departure from material linearity (due to damage or plasticity) and before crack formation, the 

adhesive modulus is effectively a continuous function across the joint, which causes more load to be 

transferred to the inner regions of the joint. However, the main idea of a FGA is that this effect can be 

achieved without taking on irrecoverable damage. Since the benefit of FGAs can be realized without 

material damage, this study will be limited to the linearly elastic regime of the adhesive. Geometric 

nonlinearity is also ignored for simplicity and because it is not expected to have a large effect on a stress 

comparison study between adhesive systems.

The analytical findings are complemented with an experimental “proof-of-concept” testing to illustrate 

the benefits of FGAs. The adhesive was graded by adding different volume percentages of glass beads, 

although no precise method was used to control the grading except the eye and hand of the person preparing 

the specimens. By showing that a joint can benefit from grading in such a rudimentary manner, the potential 

for more drastic gains through controlled and precise grading can be argued for.
NASA/CR—2012-217606 33

3.2 Method

3.2.1 Formulation

The formulation for the joint element with a functionally graded adhesive is almost identical to the 

formulation of the element with a non-graded adhesive as shown in Section 2.2. Therefore, only details that 

diverge from the previous derivation will be shown here.

The adhesive and adherend models utilized for this study are the Euler-Bernoulli adherend model 

(Section 2.2.1.1) and the simple Winkler Foundation adhesive model (Section 2.2.2.1), although any 

combination of adhesive and adherend models could be used. The adhesive moduli are now considered to be 

functions of x, causing the stresses in the adhesive to be

( ) aaa iii s e = E x 3-1

and

( ) aaa iii t g = G x .

3-2

This causes the system of governing equations to now have non-constant coefficients (compare with 

Equation 2-32)

, ( ) x u Au = x 3-3

where the coefficient matrix, A( ) x , is now a function of x. Since the coefficient matrix is non-constant, a 

new method is needed to solve the system of ordinary differential equations.

3.2.1.1 Method of constant segments for solving linear, homogeneous, non-constant coefficient system 

of ordinary differential equations 

In order to solve the system of ordinary differential equations, a semi-numerical method of solution was 

adopted. Traditional differential equation solving techniques employing numerical boundary conditions 

could not be employed because the boundaries (nodes) contain unknown, symbolic conditions. Therefore, 

the domain was split into segments in which the coefficient matrix, A(x), is considered constant and solved 

using the matrix exponential. First, consider segment n+1, with a local x-direction coordinate system x’

which originates at the left side of the segment, x=xn (Figure 3-1). The other end of the segment is at x’=Δx

and x=xn+1.

Figure 3-1. Method of linear segments involves breaking up a function into small segments

and assuming the function is constant within the segment.

u

x

u(x)

u

xn x

un

xn+1

un+1

u

x'

Δx

Δx
NASA/CR—2012-217606 34

It is assumed that Δx is sufficiently small so that A(x) can be considered constant within each such 

segment. The linearized coefficient matrix An+1, is taken to be the coefficient matrix evaluated at the 

midpoint of the segment:

1 2 ( )x

n n x A A + = + .

3-4

Within segment n+1, the system can now be expressed as a system of ordinary constant coefficient 

differential equations in the local coordinate system, x’, of the form:

u Au ,' 1 x n = + .

3-5

As before, there are 6N eigenvalues of An+1: N real eigenvalues, 2N complex eigenvalues, and 3N

repeating eigenvalues. Therefore, the solution is made up of N exponential terms, 2N exponential terms 

multiplied by a sine or cosine, and the 3N repeating eigenvalues correspond to a third order polynomial 

found in a standard beam solution. Such a complex solution shows that merely employing standard beam 

shape functions to the joint problem would be inadequate in capturing the nature of the complete solution.

The solution of the system in Equation 3-3 can be written in terms of the matrix exponential, 1 ' n x A + e , and 

a vector of unknown constants, Cn+1, as

' 1 1 ( ') x n n x + = +

A u eC .

3-6

The matrix exponential can be expressed as the infinite series [17]

' 1

0 !

k x n k

k

x

k

+

=

= A

e A

.

3-7

To get faster convergence, a method of scaling and squaring [18] is employed, and the series is calculated 

up to a value of k which yields an acceptable error, . The error can be defined many ways, but the current 

study defined the error as the difference between the 1-norms of 1 ' n x A + e for k-1 and k. The value of the 

acceptable error was set at = 0.0001.

In the local coordinate system, the solution un at x’=0 can be expressed as

01 1 1 n n nn + = = + +A ue C C ,

3-8

and the solution, un+1, at the end of the segment (x’=Δx) can be written in terms of the solution at the 

beginning of the segment,

1 1

x n n n + + = A ueu ,

3-9

to eliminate the vector of constants, Cn+1. Similarly, the solution at the end of the previous segment can be 

expressed as 

1

x n n n = A ueu 3-10

and so on, down to the first segment, which has the solution:
NASA/CR—2012-217606 35

0

0 0

x = A ueC . 3-11

Therefore, the solution at any segment, n+1, can be expressed in terms of the vector of constants from the 

first segment, C0, by the equation:

( ) 1 1 0

x xn n + + = A ue C 3-12

where

1

( ) 1

0

n x x x n m

m

+

+

=

= A A e e

.

3-13

To demonstate this method, Figure 3-2 shows a convergence study for the differential equation

2

, ( sin 4 ) xx y x e xy = - with initial conditions y(0) 1 = and , (0) 0 x y = . It can be seen that the percent 

error between the method of constant segments and the solution found using a 4th order Runge-Kutta [40]

solver with a step size of 0.001 diminishes with the number of segments. For the current study, 200 

segments were used in each element. 

It should be acknowledged that the discretization involved in this method is very similar to conventional 

finite element discretization. This may serve to reduce the advantages of using the joint element, but it still 

allows for the joint to be represented with one element, simplifying analysis steps like mesh generation and 

post-processing.

Figure 3-2. A comparison of solutions to the homogeneous, linear, non-constant coefficient differential equation

3

, ( sin 4 ) xx y x e xy = - with initial conditions y(0) 1 = and , (0) 0 x y = (a) using the method of constant segments 

with two, four, and six segments and (b) the relative error as a function of the number of segments used.

(a) (b)

0

20

40

60

80

1 3 5 7 9

% Error

Number of Segments 

0.4

0.6

0.8

1.0

0 0.2 0.4 0.6 0.8 1

y

x 

2 Segments

4 Segments

6 Segments

y(x)
NASA/CR—2012-217606 36

3.3 Results and Discussion

3.3.1 Baseline Configuration

To assess the performance of functionally graded adhesives (FGAs), a baseline joint configuration was 

identified based on the configuration in the proof-of-concept testing (Section 3.3.5). A single strap joint (or 

butt end joint) was chosen because of its ease of manufacturing, symmetry, and single dominant stress 

concentration in the middle of the joint. The geometric and material properties are defined in Figure 3-3a, 

and the finite element representation is shown in Figure 3-3b. Only half of the joint was modeled due to 

symmetry, and the overlap section was modeled by one or several joint elements (depending on the

number of discrete regions of continuous adhesive modulus) while the non-joint adherend section is 

modeled with one beam element. The loaded end is constrained from rotation and vertical translation, while 

the symmetric face of the doubler is constrained from horizontal translation and rotation.

The values of the material and geometric parameters used for the analytical modeling are found in

Table 3-1, and are based on realistic values that coincide with the proof-of-concept test specimens discussed 

in the next section. The FGAs were compared with two different single adhesive systems (Figure 3-4a).

These two adhesives provided upper and lower bounds for the grading functions, and will be referred to as 

Eu and El, respectively. The grading functions chosen for investigation included a step (bi-adhesive), linear, 

and exponential function. These functions were reduced to a single grading variable, l, as defined in

Figure 3-4b, c, and d. For the step function, l is the length of the more compliant adhesive. For the linear 

function adhesive, l is the length of the section in which the modulus decays linearly, and for the exponential 

function it is the length of adhesive which has a modulus less than 99% of Eu. Although these may not be 

the optimal grading functions, these functions were chosen because of their ability to be reduced to one 

variable, allowing for simple and clear sensitivity studies. It was assumed that regardless the grading, the 

relationship between the Young’s modulus and shear modulus remained constant or, in other words, that the 

Poisson’s ratio remained constant. A similar assumption was made by Apetre et al. [41] for functionally 

graded sandwich beam cores.

Figure 3-3. Single strap joint (a) geometric and material parameters and

(b) joint finite element representation assuming symmetry.

Table 3-1. Parameters of the baseline single strap joint configuration used for the theoretical study.

P

(kN)

la

(mm)

lo

(mm)

t

(mm)

b

(mm)

η

(mm)

E

(GPa)

Eu

(GPa)

El

(GPa)

ν τnet

(MPa)

4000 82.6 38.1 1.1 25.4 0.4 108.5 2.5 1.1 0.34 8.3

P

lo

la b

t η

P

(a)

(b)

E

Ga, Ea

x

z
NASA/CR—2012-217606 37

Figure 3-4. Single strap joints with (a) constant modulus adhesives were compared with joints with functionally 

graded adhesives, including (b) step-wise graded, (c) linearly graded, and (d) exponentially graded.

Most stress values reported in the theoretical study were normalized by the net shear stress, τnet, defined 

by:

net

o

P

l b = 3-14

3.3.2 Stress Comparison

The first goal of this study was to address the concern that functionally graded adhesives (FGAs) are not 

worth the increase in manufacturing difficulty. This section addresses this concern by showing the potential 

adhesive peel stress reduction due to FGAs over single adhesive joints. This will be shown for the three 

functions defined previously (Figure 3-4). Next, the result of changing the lower bound modulus, El, will 

show how further stress reductions can be achieved. Finally, since stepped adhesives will most likely be the 

easiest FGAs to manufacture, the effect of increasing the number of steps and spacing optimization will be 

investigated. All optimization was carried out using a genetic algorithm.

The grading parameter, l, was optimized to reduce the maximum peel stress in the adhesive for the three 

FGAs, and the resulting moduli are plotted in Figure 3-5a. As shown, the region of gradation is very small, 

about 2% of the overlap length for the step and linear function adhesives, and around 5% for the 

exponential. The peel stress in the adhesive for half of the symmetric joint is plotted in Figure 3-5b for each 

adhesive. The single adhesive joints are in blue, and the FGAs are in black. The step FGA has two stress 

peaks; one at the end of the adhesive and one at the interface between the two adhesives. The linear and 

exponential FGAs have a rounded stress peak, and appear to result in very similar stress distributions.

0.0

0.5

1.0

0 0.25 0.5

Ea / Eu

x / lo

Eu

El

(a)

0.0

0.5

1.0

0 0.25 0.5

Ea / Eu

x / lo

l/lo

(b)

0.0

0.5

1.0

0 0.25 0.5

Ea / Eu

x / lo

l/lo

(c) 

0.0

0.5

1.0

0 0.25 0.5

Ea / Eu

x / lo

l/lo

99% Eu

(d)
NASA/CR—2012-217606 38

Figure 3-5. Optimized configurations for the single strap joints for different functions of graded adhesive compared: 

(a) modulus across the adhesive and (b) centerline peel stress across the adhesive.

The maximum stress in the adherend, doubler, and adhesive for the single adhesive joints and the 

optimized FGA joints is found in Figure 3-6. All stress values are normalized by the maximum stress found 

in the stiffer single adhesive for viewing all stress components in the same plot. The stress reported for the 

adherend and the doubler is the normal stress in the x-direction, and the maximum value of the stress is 

found at the upper and lower surfaces of the adherend and doubler, respectively. With composite laminate 

adherends, the most important stress in the adherend is usually the peel stress (z-direction) between the plies 

because this is where failure often initiates. However, the current configuration contains only one ply, so 

there is no such interlaminar shear stress.

Some important aspects of using FGAs are illustrated in Figure 3-6. First, the FGAs in this study 

outperformed the stiffer single adhesive joint, El. Adhesive stresses were considerably lower and adherend 

and doubler stresses were not significantly impacted. This is important because the more compliant single 

adhesive joint had lower adhesive stresses, but higher adherend stress. The FGAs were able to lower the 

adhesive stress without affecting the adherend stress. Second, when compared with the El adhesive, the 

FGAs reduced the adhesive peel stress but did not necessarily reduce the shear stress. It will be shown later 

that the optimum value of the grading variable l is not the same for minimizing peel as it is for minimizing 

shear stress. Therefore, the relative levels of peel and shear must be considered when designing FGAs so 

that the dominant stress can be minimized. However, typical adhesives are more ductile under shear loading, 

so peel stresses will normally be the minimized variable.

-6

4

14

24

0.45 0.475 0.5

σz / τnet 

x / lo

-10

0

10

20

30

0 0.25 0.5

FM 300

AF 163-2k

Mix

Linear

Exponential

σz / 

τnet 

x / lo

0

0.5

1

0 0.25 0.5

Upper

Lower

Discrete

Linear

Exponential Ea / Eu

x / lo

0.4

0.6

0.8

1

0.45 0.475 0.5

Ea / Eu

x / lo

(a)

(b)

Constant (Eu )

Constant (El )

Step

Linear

Exponential

Constant (Eu )

Constant (El )

Step

Linear

Exponential
NASA/CR—2012-217606 39

Figure 3-6. Comparison of the maximum stresses in the joints with different adhesives, where all maximum stress 

values are normalized by the maximum stress value found in the stiffer constant modulus (single adhesive) joint,

and adherend and doubler stresses refer to those in the x-direction.

Figure 3-7. The maximum stresses in joints with linear FGAs and different lower bound modulus, El, where all 

maximum stress values are normalized by the maximum stress value found in the stiffer constant modulus (single 

adhesive) joint, and adherend and doubler stresses refer to those in the x-direction.

In order to further reduce the maximum adhesive peel stress, the lower bound modulus can be 

decreased. The affect of changing the lower bound modulus on the maximum joint stresses is shown in 

Figure 3-7 for the linear FGA only (Figure 3-4c). For each value of El, the grading parameter l has been 

selected to produce for the smallest maximum peel stress. By grading the adhesive from a very low modulus 

to a higher one, the peel stress can be lowered significantly, with peel stress reductions in this case of up to 

90%. At the same time, the shear stress is also significantly lowered. For low lower bound modulus values, 

the maximum adherend and doubler stresses are also impacted. The maximum adherend stress increases 

with decreasing El, while the maximum doubler stress decreases slightly. While it may appear that the 

0.0

0.5

1.0

1.5

Adherend Doubler Normal Shear

FM 300

AF 163-2k

Mix

Linear

Exponential

Adhesive

σmax/σEu

Constant (Eu )

Constant (El )

Step

Linear

Exponential

0

0.5

1

1.5

0 0.5 1

El/Eu

σmax/σE

u

Adherend

Doubler

Adhesive Normal

Adhesive Shear
NASA/CR—2012-217606 40

Figure 3-8. Increasing the number of adhesives used for a discretely graded adhesive and optimizing the

adhesive spacing lowers the maximum adhesive peel stress; (a) stress vs. number of adhesives, and

(b) modulus at the end of the adhesive for 20 steps compared with the other functions.

adherend stress is raised considerably while the doubler is only slightly decreased, it should be kept in mind 

that these values are normalized by the stresses found in the stiffer single adhesive system. The maximum 

doubler stress is around three times higher than the maximum adherend stress, so the stress is being taken 

from the highly loaded doubler to the less loaded adherend, which is advantageous when they are all made 

of the same material as in this case.

Since it is likely that grading of adhesives will be accomplished using a step function with different 

adhesives for each step, the effect of the number of steps used on the adhesive peel stress is plotted in

Figure 3-8a. For each point, it was assumed that the moduli of the steps are equispaced between El and Eu.

The length of each step was optimized using a genetic optimization routine to reduce the peel stress. Using 

more steps led to a lower maximum peel stress, eventually even lower than the linear and exponential 

functions. The optimized moduli for 20 steps is shown in Figure 3-8b, along with the optimal linear and 

exponential moduli. This sheds a little light on a better grading function, consisting of a constant modulus 

region followed by an exponentially higher order declining region.

This section has shown that there are significant stress reductions possible by using FGAs. Unlike 

simply using a lower modulus adhesive, FGAs do not have a significant effect on adherend and doubler 

stresses. Furthermore, lowering the lower bound modulus leads to even greater stress reductions. These 

stress reductions are possible using stepped modulus adhesives, although the length of each step needs to be 

optimized.

3.3.3 Sensitivity Study

The second goal of this study was to address the concern that during manufacturing, the adhesive is 

pressurized and heated, often causing the adhesive to flow and even squeeze out of the joint. If a 

functionally graded adhesive (FGA) is specifically designed for a certain joint, such squeeze-out could 

either change the shape of the grading, the lower bound modulus, or both. This could result in an FGA 

which has higher stress than using the more compliant adhesive alone. This section seeks to address this 

concern by presenting a grading sensitivity study. For the purposes of this study, sensitivity will refer to the 

narrowness of the range of l values that results in lower maximum stresses than those obtained using the 

lower bound adhesive, El. First, the effects of changing the grading parameter l on the maximum adhesive 

stress are quantified (Figure 3-9), followed by an investigation into changing the lower bound modulus El

and l (Figure 3-10). Through this study, insight is gained into the sensitivity of the grading and which 

functions are more tolerant to changes in the grading.

12

14

16

18

2 8 14 20

σmax / τnet 

Number of Adhesives (b)

0.4

0.6

0.8

1

0.47 0.48 0.49 0.5

Upper

Lower

Linear

Exponential

20

x/lo

Ea / Eu

Constant (Eu )

Constant (El )

Linear

Exponential

Step (20 Steps)

(a)

Constant (El )

Step 

Linear

Exponential
NASA/CR—2012-217606 41

Figure 3-9. Observing the maximum adhesive (a) peel and (b) shear stress as a function of l as defined in Figure 3-4

shows how sensitive the maximum stress is to the shape of the grading.

The effects of changing the grading parameter l on the maximum adhesive peel and shear stresses for 

the three FGAs are found in Figure 3-9a and b respectively. For all three FGAs, the maximum shear stress is 

less sensitive to the grading parameter than the peel stress. Also, the optimum value of l is always greater for 

minimizing the shear stress than the peel stress. This is most likely because the peel stress peak is much 

more concentrated than the shear stress peak, so a steeper gradation is needed to minimize the peak. Also, if 

l becomes too short, the stress goes above the stress that would be found in a joint with just the more 

compliant adhesive. Unfortunately, with the addition of pressure, a decreasing l is more likely. The step 

function adhesive was the most sensitive: only a very small range of values of l results in lower stresses than 

just using the more compliant adhesive, El. The linear function adhesive was not as sensitive, and it 

converges to a stress less than that of El when l is large. Finally, the exponential function adhesive had a 

broad range of l values resulting in low maximum stress, making it the most tolerant to changing the grading 

parameter.

When the joint is being manufactured, the adhesive is heated and pressure is applied. The most likely 

result will be that the adhesive will spew out and the modulus of the adhesive at where the adherend ends 

will be greater. At the same time, the grading parameter l may also change. Therefore, it is important to 

know how sensitive the stress is to increasing the modulus and how that affects the grading parameter 

sensitivity. Based on Figure 3-10, the sensitivity of the maximum adhesive stress for all three FGAs to l

increases with increasing lower bound modulus El. Additionally, the optimum grading parameter l decreases 

as the lower bound modulus increases. Therefore, it is recommended that while using FGAs in joints, one 

14

18

22

26

0 0.25 0.5

FM300

Af 163-2k

Dis

Lin

Exp

σmax / 

τnet 

l/lo

Constant (Eu )

Constant (El )

Step

Linear

Exponential

7

9

11

13

0 0.25 0.5

FM300

Af 163-2k

Dis

Lin

Exp

l/lo

-τmax / 

τnet 

Constant (Eu )

Constant (El )

Step

Linear

Exponential

(a)

(b)
NASA/CR—2012-217606 42

Figure 3-10. The sensitivity of the maximum adhesive peel and shear stress to decreasing the lower bound modulus, 

El, for the discrete (a and b), linear (c and d), and exponential (e and f) FGAs.

should either figure out how much flow will occur and plan accordingly, or use a grading system which will 

not move when adhesive flows. An example of such a system is the glass beads used for the proof-ofconcept testing which will be described in Section 3.3.5. Since beads are used for thickness control also, 

applying pressure causes the beads to be clamped between the adherend and doubler and unable to move.

Such a system of grading is nearly independent of adhesive flow.

3.3.4 Effect of Loading

The third goal of this study was to show the effect of changing the loading conditions when using FGAs 

in joints. In many applications of functionally graded materials, a change in loading can lessen or even 

reverse the stress gradient, rendering the material grading useless or even detrimental. To show that this is 

not the case with adhesively bonded joints, six different loading scenarios and their respective adhesive peel 

stress distributions are shown in Figure 3-11. The boundary conditions were altered from that of the baseline 

0

5

10

15

0 0.1 0.2 0.3 0.4 0.5

-τmax / τnet 

l/lo

0

10

20

30

0 0.1 0.2 0.3 0.4 0.5

σmax / τnet 

l/lo

0

5

10

15

0 0.1 0.2 0.3 0.4 0.5

l/lo

-τmax / τnet 

0

10

20

30

0 0.1 0.2 0.3 0.4 0.5

l/lo

σmax / τnet 

0

5

10

15

0 0.1 0.2 0.3 0.4 0.5

-τmax / τnet 

l/lo

0

10

20

30

0 0.1 0.2 0.3 0.4 0.5

σmax / τnet 

l/lo

1.0

0.8

0.6

0.4

0.2

0.0

El

Eu

1.0

0.8

0.6

0.4

0.2

0.0 El

Eu

1.0

0.8

0.6

0.4

0.2

El 0.0

Eu

1.0

0.8

0.6

0.4

0.2

0.0 El

Eu

1.0

0.8

0.6

0.4

0.2

0.0

El

Eu

1.0

0.8

0.6

0.4

0.2

0.0

El

Eu

(a)

(c)

(e)

(b)

(d)

(f)
NASA/CR—2012-217606 43

Figure 3-11. For different loading scenarios at the end of the joint ( (a) axial load, (c) moment, and

(e) vertical load) the stress gradients of the adhesive peel stress (b, d, and f) remains the same direction;

high at the ends and low in the middle. 

configuration (Figure 3-3) to represent an end moment and shear force. The end of the doubler was 

clamped, with no other constraints. The peel stress in the adhesive layer for a tensile and compressive axial 

load, positive and negative moment, and up and down vertical shear force is shown in Figure 3-11b, d, and f 

respectively, while the loading, deformed configuration, and boundary conditions are shown in

Figure 3-11a, c, and e.

For all of these loading cases, the gradient of stress goes from high stress (tensile or compressive) at the 

ends to low stress in the middle. Therefore, functionally grading the adhesive in the manner prescribed will 

always be beneficial under mechanical loading; thermal loading conditions will be addressed in the future. 

Compressive peel stresses are not considered to initiate or propagate cracks in an adhesive layer. 

Therefore, when the loading results in a compressive peel stress at the end of the joint, the subsequent stress 

“trough” (see Figure 3-5 for an example of a stress distribution and trough), which is comprised of tensile 

stress, becomes the critical stress. Although this trough is always significantly lower in magnitude than the 

peak stress at the end of the joint, the gradation should be designed to minimize the trough. By minimizing 

the trough rather than the peak stress at the end of the joint, the strength of a joint with a compressive peak 

stress at the end of the adhesive can be maximized.

To simplify the parametric studies, the FGAs investigated were only designed to minimize the stress 

concentration in the middle of the joint (x/lo=0.5) because this was the dominant stress concentration in the 

baseline configuration. However, Figure 3-11d and f shows that changing the loading can cause both ends to 

be highly stressed. Therefore, it is highly recommended that the grading be applied to every area of a joint 

where there will be a stress concentration; i.e. at the ends of the adherends and doubler. Applying grading to 

all ends of the joint will ensure that the FGA will still be optimal regardless of loading conditions.

Although Figure 3-11 does not explore all possible loading conditions, it is believed that these represent 

a large and varied enough sample to draw general conclusions. Different loading conditions away from the 

actual joint overlap region can always be resolved into a shear load, axial load, and moment near the joint 

overlap as has been done for many classical analytical formulations [21]. Therefore, changing the loading at 

the end of the adherend results in merely changing the relative magnitudes of the load components near the 

joint overlap region.

-1

0

1

0 0.25 0.5

σz / σmax 

x / lo

-1

0

1

0 0.25 0.5

σz / σmax 

x / lo

σz / σmax -1

0

1

0 0.25 0.5

σz / σmax 

x / lo

(a)

(c)

(e)

(b)

(d)

(f)

P

-P

M

-M

-V

V

P

-M

-V

-P

M

V
NASA/CR—2012-217606 44

Figure 3-12. Maximum peel and shear adhesive stress and its sensitivity to grading parameter, l, for different 

amounts of end moment (a, b) and end shear force (c, d) for joints with a linear FGA.

To study the effect of changing the loading on the optimal grading shape and grading sensitivity, the 

normalized maximum normal (Figure 3-12a, c) and shear (Figure 3-12b, d) stress in the adhesive for joints 

with the linear FGA with different grading shapes was plotted. Figure 3-12a and b show the effect of 

different ratios of an end vertical shear force to axial force, and Figure 3-12c and d show the effect of 

different ratios of end moment to axial force. The loading was chosen so that the maximum stress in the 

adhesive remained at x/lo=0.5. As can be seen, the loading investigated here did not affect either the 

optimum grading or the sensitivity of the grading shape. Therefore, it can be concluded that the same 

grading would be beneficial for multiple loading scenarios.

In addition to varying the loading scheme, the joint type was also altered to show that the adhesive peel 

stress magnitude gradient always goes from high stress near the adherend or doubler discontinuities to low 

stress elsewhere. The normalized adhesive peel stress for a single lap joint, bonded doubler, and double 

strap joint for different loadings is shown in Figure 3-13a, b, and c, respectively. This figure not only shows 

that the adherend peel stress magnitude gradient remains in the same direction for different joint types, but 

also suggests that the normalized stress gradient is very similar to that of the single strap joint. It might be 

the case that the same adhesive gradient could be utilized for not only different loading conditions, but 

different joint types. If it was the case that one grading type was beneficial for many different joint types and 

geometries, graded adhesives could be mass produced for general application rather than custom designed 

for each individual joint type and geometry.

3.3.5 Proof-of-Concept Testing

To illustrate the potential gains of FGAs, proof-of-concept testing was performed on single strap joints 

(the baseline case). The adherends and doubler consisted of a single ply of 0/±45 triaxially braided 

composites, with the axial tows consisting of around 80,000 fibers, biased tows with 12,000 fibers, a fiber 

volume fraction of 52%, and matrix of Epon 862 epoxy resin [42]. The axial 0° direction was aligned with 

the global x-direction as defined in Figure 3-3. The adhesive used was AF 163-2k [43] and the bond line 

was two layers thick in order for glass beads, normally used for thickness control, to be inserted in between 

(Figure 3-14).

10

20

30

40

0 0.1 0.2 0.3 0.4 0.5

σmax / τnet 

l/lo

5

10

15

20

0 0.1 0.2 0.3 0.4 0.5

-τmax / τnet 

l/l (a) (b) o

V

P

V

P

0.01 0.008

0.006

0.004

0.002

0.0

0.01

0.008

0.006

0.004

0.002

0.0

5

10

15

20

25

0 0.1 0.2 0.3 0.4 0.5

-τmax / τnet 

l/lo

10

20

30

40

50

0 0.1 0.2 0.3 0.4 0.5

σmax / τnet 

l/lo

1.0

0.8

0.6

0.4

0.2

0.0

M

P

1.0

0.8

0.6

0.4

0.2

0.0

M

P

(c) (d)
NASA/CR—2012-217606 45

Figure 3-13. For different loading scenarios for a (a) single lap joint, (b) bonded doubler, and

(c) double strap joint, the stress gradients of the adhesive peel stress (b, d, and f)

remains the same direction; high at the ends and low in the middle.

Figure 3-14. Diagram of single strap joint with placed glass beads to test functionally graded adhesive concept.

-1

0

1

0 0.5 1

σz / σmax -1

0

1

0 0.5 1

σz / σmax 

x / lo

-1

0

1

0 0.25 0.5

σz / σmax 

x / lo

σz / σmax 

x / lo

σz / σmax 

x / lo

(a)

(b)

P

M

P

M

M M

V

(c)

M

V

P

M M

P

M

Doubler

Adhesive

Adherend

Adherend Glass Beads

x / lo
NASA/CR—2012-217606 46

All specimens were manufactured together as a plate, with individual specimens cut out of the plate by 

means of a wet saw with a diamond-coated blade. The adherends and doubler were pre-infiltrated and cured 

prior to bonding. After bonding, the adhesive was cured in a heated press at 177° C under 0.43 MPa for 120 

minutes, and specimens were not handled for 24 hours. The average dimensions of the specimens along with 

the standard deviation are shown in Table 3-2. The adherend material properties used in Section 3.3.1 are

representative of the adherends in this study, while the upper bound modulus is representative of the pure 

adhesive modulus. The specimens were placed in grips and pulled in tension under displacement control at a 

rate of 0.005 mm/s. Load and displacement data was recorded, and the joint strength was defined as the 

maximum load held by the joint.

The grading was accomplished by means of inserting glass beads between the adhesive layers to change 

the properties of the adhesive. Bead placement was done purely by hand, and a photograph of the adhesive 

prior to doubler bonding of a representative specimen for each adhesive system is shown in Figure 3-15. 

Two specimens were made with no beads in the adhesive, two with a high density of beads uniformly 

spread throughout the adhesive region, and three with beads placed strategically by hand. Care was taken to 

ensure that the location of the stress concentration, where the two adherends meet, was devoid of beads with 

a gradual increase in beads going out from the center.

An image of a single strap joint prior to complete joint separation is shown in Figure 3-16a. The 

resulting strengths of the joints tested are plotted in Figure 3-16b. Each bar represents a specimen; with two 

specimens with no beads, two with a uniform distribution of beads, and three specimens with graded beads 

in the adhesive. The specimens with no beads performed the worst, with an average strength of 3.18 kN. 

The uniform beaded joints had an average strength of 4.29 kN, and the graded bead specimens had an 

average strength of 5.43 kN. With so few specimens, no statistically significant conclusions can be drawn, 

but trends suggest that grading the beads could have resulted in higher strengths over the uniform 

specimens.

Table 3-2. Geometric parameters for proof-of-concept testing.

la (mm) lo (mm) t (mm) b (mm) η (mm)

68.76±0.38 37.49±0.25 1.08 ±0.02 22.73±0.17 0.42±0.03

Figure 3-15. Photographs of the adhesive layer before placement of the doubler for the single strap joint with

(a) no beads, (b) uniform beads, and (c) graded beads.

1 mm 1 mm 1 mm (a) (b) (c)

border 

between 

adherends
NASA/CR—2012-217606 47

Figure 3-16. Single strap joint test: (a) photograph of typical specimen prior to failure and

(b) joint strengths of the specimens tested with different adhesive systems.

Figure 3-17. Post mortem photographs of the failed adhesive layer for the single

strap joints with (a) no beads, (b) uniform beads, and (c) graded beads.

A few words should be said about the effects of FGAs and failure modes. As was shown earlier in 

Figure 3-6, grading the adhesive has little or no effect on the adherends for the configuration studied.

Therefore, a difference in strengths was only manifested when the failure occurred in the adhesive, not in 

the adherend or doubler.

For these specimens, an improvement appears to exist using FGAs because the failures consisted mostly 

of cohesive failure with a few areas of adhesive failure (Figure 3-17). Another round of testing was 

performed with the same type of specimens cured at a lower temperature. These specimens all failed in the 

adherend, and there was no statistical difference in strengths between graded and non graded specimens.

0

2

4

6

No Beads Uniform Beads Graded Beads

Joint Strength (kN)

(a) (b)

5 mm

Doubler

Adherend

1 mm (a) (b) (c) 1 mm 1 mm
NASA/CR—2012-217606 48

Therefore, it should be noted that grading the adhesive can improve joint strength when failure occurs in or 

around the adhesive, and may even drive a joint to fail in the adherend when ungraded specimens fail in the 

adhesive, but FGAs most likely will not have a significant effect when failure already occurs in the 

adherends without grading. However, it is desirable for joints to fail in the adherend and adhesive at the 

same time so that the joint is maximized for strength, while minimized for weight.

3.4 Concluding Remarks

The peel stress concentration found in adhesively bonded joints can be reduced by grading the modulus 

of the adhesive material. More recently, grading of the joints in a continuous manner has become of interest.

However, the actual use of such joints has been limited because of a few concerns which this study has 

addressed. These concerns were addressed by comparing the stresses in single adhesive, bi-adhesive, linear, 

and exponential functionally graded adhesives (FGAs) using an analytical model. This unique model allows 

efficient analysis of adhesively bonded joints in the context of a finite element framework without meshing 

the 2-D or 3-D details of the joint geometry. The efficiency of the model facilitated vast parametric and 

optimization studies, while still having the ability to be placed within a larger, global structural-scale model.

The first concern was that gains to be made using FGAs were not worth the added manufacturing 

complications. It was shown that for the single strap joint configuration investigated, FGAs could reduce the 

maximum peel stress in the adhesive by up to 17% over the more compliant single adhesive joint without 

having an adverse effect on the adherend stress. Since the optimum grading for shear stress was not the 

same as that for peel stress, there was a slight increase in shear stress when peel stress was minimized.

However, all grading functions resulted in a significant reduction (> 45%) in peel stress over the stiffer 

single adhesive joint without adversely affecting the load carrying capability or the stress in the adherends.

Furthermore, decreasing the lower bound modulus of the grading caused even more dramatic stress 

reductions. Finally, it was shown that peel stress reductions comparable to those of the continuous functions 

studied here could be reached using a stepped function with around five or more different modulus 

adhesives, which is probably easier to manufacture.

The second concern about using FGAs was that the flow of adhesive during manufacturing would 

change the shape of the grading and cause the grading to be ineffective. Therefore, a sensitivity study was 

conducted on the three FGAs to see the effect of changing the shape of the grading (l) and the lower bound 

modulus (El) to reflect what might occur when adhesive is squeezed out of the joint. The maximum adhesive 

peel stress levels did increase for all FGAs when the lower bound modulus increased, but the exponential 

FGA proved to be quite tolerant to changes in grading shape. This could be the basis for a justification for 

using the more complicated exponential FGA over the bi-adhesive. Although the peel stress reductions were 

all very similar, the tolerability of the exponential grading to perturbation of grading shape which can occur 

due to adhesive flow might make it worth the extra complications. Also, it is recommended that care be 

taken to account for adhesive flow by adding more lower modulus adhesive at the end of the joint than is 

required, or using a grading system that does not change when the adhesive flows.

Third, functionally graded materials used in structures can become useless or detrimental when the 

loading conditions are changed. It was shown that changing the loading configuration does not change the 

stress magnitude gradient of the adhesive in the case of a single strap joint. Consequently, the optimum 

grading shape and sensitivity was shown not to be significantly impacted by changing the loading 

conditions. However, in the case of moment and shear loading, both ends of the joint are subjected to high 

peel stresses, thus suggesting for generality that both ends be always graded. It was also shown that the 

stress magnitude gradient remains the same for other joint types; high stress at the adherend discontinuities 

with low stress elsewhere; thus suggesting the possibility of using the same adhesive grading for a myriad of 

joint configurations. However, additional parameters such as adherend stiffness, overlap length, and 

adherend thickness must be taken into account to conclusively determine the universality of a specified 

grading.
NASA/CR—2012-217606 49

Additionally, proof-of-concept testing was conducted to show the potential of joints with FGAs.

Grading was done by strategically placing glass beads to change the stiffness of the adhesive along the joint.

Although not enough specimens were tested to result in significant conclusions, trends seem to indicate that 

these graded joints held over more load than specimens with no beads or uniformly distributed beads for 

failures dominated by cohesive and adhesive failures. Additionally, as expected, it was found that grading 

did not result in an increase in strength when the failure without grading occurred in the adherend.

Through theoretical and experimental means, this study has shown the practical benefits of and offered 

some guidelines on how to effectively use and design joints with FGAs. FGAs have the potential to spread 

the stress more evenly across the adhesive without unrecoverable damage, making them another means of 

improving adhesively bonded joints.

NASA/CR—2012-217606 51

Chapter 4

Progressive Failure of Adhesively Bonded Composite Joints

Enhanced finite elements are elements with an embedded analytical solution that can capture detailed 

local fields, enabling more efficient mesh independent finite element analysis. In Chapter 2, this method was 

applied to adhesively bonded joints. The adherends were modeled as composite Euler-Bernoulli beams, and 

the adhesive layer was modeled as a bed of linear shear and normal springs. The field equations were 

derived using the principle of minimum potential energy, and the resulting solutions for the displacement 

fields were used to generate shape functions and a stiffness matrix for a single joint finite element. In the 

current chapter, the capability to model non-linear adhesive constitutive behavior with large rotations is 

developed, and progressive failure of the adhesive is modeled by re-meshing the joint as the adhesive fails.

The joint element is compared with experimental results for various joint configurations, including double 

cantilever beam and single lap joints.

4.1 Introduction

The joint element has been shown to predict adhesive stresses quite well under linear elastic conditions, 

but these conditions are not sufficient to predict joint strength. Modern polymeric adhesives are usually 

highly nonlinear, causing linear elastic analysis to be insufficient. Furthermore, the eccentricity of many 

joint configurations results in large rotations early on in the loading [44]. Large rotations are also necessary 

to predict the onset of imperfection-induced buckling, an aspect of great concern in thin-walled aerospace 

structures.

Therefore, geometric nonlinear effects due to large rotations and material nonlinearity are pivotal in 

predicting the strength of a joint. This chapter will extend the previously created joint element to include 

these effects. Additionally, methods of growing an adhesive crack and adapting shape functions during the 

analysis will be presented in order to preserve the original intent of the joint element, which is to model a 

joint with very few elements.

Currently, the scientific community seems to model the progressive failure of joints using fracture 

mechanics methods like LEFM and cohesive zone models, or continuum mechanics, both with dense-mesh 

finite elements [45]. Since the joint element is merely a tool, it will accommodate using inputs derived from 

either of these philosophies to govern the stress-strain relation of the adhesive. A method will be shown of 

characterizing the adhesive layer using either bulk adhesive tensile data as would someone using continuum 

mechanics damage progression, or fracture mechanics inputs like mode I strength and fracture toughness.

The application of each will be demonstrated and results will be compared with experiments.

4.2 Formulation

The formulation of the joint element has been broken up into discrete parts, namely the co-rotational 

formulation, material nonlinearities, crack growth, adaptive shape functions, and adhesive constitutive 

modeling. Each section presents a formulation to address a certain aspect of the progressive failure of the

joints. The co-rotational formulation addresses large rotations in joint problems while material nonlinearities 

show how nonlinear constituents are modeled. The crack growth formulation deals with the failure of the 

adhesive layer and adaptive shape functions are used to reduce the number of joint elements required to 

model nonlinear adhesive behavior. Finally, the last section illustrates a few methods of defining the 

properties of the adhesive based on several different experimental techniques.

4.2.1 Co-rotational Formulation

Consider a structure consisting of N layers of thin plates under cylindrical bending joined together by

N-1 thin layers of a much more compliant adhesive material (see Figure 2-1). The plates are assumed to 

behave as “wide” Euler Bernoulli beams (hence the cylindrical bending assumption). The adhesive joining 
NASA/CR—2012-217606 52

the plates is modeled as a Winkler foundation. The plates can be isotropic, transversely isotropic, or a 

layered composite. The plates and adhesive are assumed to be under proportional loading, and are modelled 

as nonlinear elastic materials.

A co-rotational formulation will be used to capture large rotations, and has been primarily adapted from

prior work by Belutschko and Hsieh [46] and Crisfield and Moita [47]. This formulation tracks the rigid 

body rotation of an element through a local rotational coordinate system, and considers the rotations and 

deformations measured with respect to this rotated frame of reference to be small (Figure 4-1). The main 

benefit of this formulation is that the previously implemented code for the small rotation problem (see 

Section 2.2) can be utilized in subsequent calculations.

The element has 2N nodes located at the boundaries of the centerline of the plates (numbered as shown 

in Figure 2-2), and the nodal displacements are defined as:

1 T i N = é ù

ë û qq q q K K 4-1

where the superscript represents the adherend or plate, and 

1 231 2 3

i il il il ir ir ir = é ù qqqqqq ë û q 4-2

refers to the horizontal, vertical, and rotational displacements of the left and right nodes in plate i 

respectively.

4.2.1.1 Rigid Body Displacements

The element has a local rotated coordinate system, xˆ , which is rotated and translated relative to the 

fixed coordinate system, x , by an angle and a vector 1 qt respectively (Figure 4-1). The translation and 

angle will be properly defined later. The nodal displacements of the element in the fixed coordinate system

can be decomposed into rigid body displacements, qrig , and displacements which only cause deformation 

in the body, qdef , in the relation:

qq q = + rig def . 4-3

Figure 4-1. The nodal displacements can be broken up into two parts: a) rigid translation and rotations and

b) local deformations.

a) b)

z

x

zˆ xˆ

zˆ

xˆ

1 qt
NASA/CR—2012-217606 53

The rigid body nodal displacements, qrig , can be further decomposed into rigid body displacements 

resulting from rigid body rotation, qr , and displacements resulting from rigid body translation, qt :

q qq rig t r = + . 4-4

The translational displacements, qt , are defined as:

1 1 T

tt t = é ù

ë û qq q K 4-5

which is the horizontal and vertical displacements of the left node of the first plate and the rotation of the 

first adherend:

1 11

t 1 2 = é ù q q ë û q . 4-6

Although the rotation is not necessarily part of the rigid body translation, it is more convenient to insert it 

into the translational rigid body displacements because each adherend will be rigidly rotated by the angle .

To find the rigid body displacements due to the rotation of the element about the first node, consider the 

right node of the ith adherend, node ir (Figure 4-2). Initially, node ir can be located relative to the first node 

by a position vector ir x . When the element rotates about the first node by the angle , its new position 

relative to the first node can be expressed by an orthogonal transformation matrix as 1

T T xir , where 

1

0

0

0 01

c s

s c

é ù

ê ú = -ê ú

ê ú

ë û

T

4-7

and s and c denote the sine and cosine of the angle . Therefore, the displacement vector, ir qrot , of node ir

due to rigid body rotation can be expressed as

1 ( ) ir T q T Ix rot = - ir . 4-8

Figure 4-2. Displacements, ir qrot , of node ir due to a rigid body rotation of the joint element.

ir x

ir qrot

1

T T xir

node ir
NASA/CR—2012-217606 54

Translating this to all nodes and combining with Equation 4-4, the displacements due to rigid body rotation 

are

( ) T q q T IX rig t =+ - 4-9

where

1

1

é ù

ê ú = ê ú

ê ú

ë û

T

T

T

O

4-10

and the vector X is simply a collection of the initial x and z coordiantes of the nodes, and is defined 

explicitly as

[ 1 1 ]

T

XX X XX X X = l r K K il ir Nl Nr

4-11

where the first subscript identifies the plate number, and the following letter, either l or r, refers to the left or 

right node respectively. The nodal coordinate vector for the ith adherend and the left node is defined as

Xil il il = [ x z 0] 4-12

while the coordinate vector of the right node is defined in an identical fashion.

4.2.1.2 Determination of the Rotation Angle

If the rotation is not constant within the joint, the rotation angle is an approximation. Adhering to the 

conventional approach for co-rotational beam formulations, the rotation angle was chosen to be the rotation 

of the first adherend as shown in Figure 4-1a. To find the transformation matrix of Equation 4-7, the sine 

and cosine of the rotation angle can be expressed as

1

1

sin z l

s

l = = 4-13

and

1

1

cos x l

c

l = = 4-14

which are defined in terms of the nodal displacements of the first adherend by

1 1

1 11

r l

x l lq q =+ - , 

1 1

1 22

l r

z l qq = -

4-15

where l is the original length of the element and l1, 1x l , and 1z l refer to the current length of the 1st adherend

and the length decomposed into x and z components Figure 4-3.
NASA/CR—2012-217606 55

Figure 4-3. The initial and current lengths of the 1st adherend are used to determine the rotation angle.

4.2.1.3 Local Coordinate System

First, the internal force vector and stiffness matrix will be found in the local, rotating coordinate system.

The stress and strain of the adherends and adhesive are assembled together in one stress and one strain 

vector as shown:

1 1 1 ( 1) ˆ ˆˆ ˆˆ ˆ ˆ ˆ T

a i ai N a N N - - = é ù

ë û σ σσ σσ σ σ σ K K 4-16

and 

1 1 1 ( 1) ˆ ˆˆ ˆˆ ˆ ˆ ˆ T

a i ai N a N N - - = é ù

ë û ε εε εε ε ε ε K K 4-17

where the overbar caret denotes quantities in the local rotating coordinate system, and the adhesive/adherend 

strain and stress vectors are as defined in Equations 2-2, 2-3, 2-12, and 2-13. Using beam theory and 

assuming small strains from the rotated coordinate system, the strains are related to the adherend centerline 

displacements,uˆ , by the equation

ε Gu ˆ ˆ = 4-18

where the adherend centerline displacements are a collection of a vector of the centerline displacements of 

each adherend layer defined in Equation 2-33. Additionally, G is an assembly of the contributions of the 

adherend and adhesive layers assembled in the form

1

l

1

2

r q

1

1

r q

1

1

l

q

1

2

l

q

l
NASA/CR—2012-217606 56

[ ]

[ ]

[ ]

[ ]

[ ]

[ ]

1

1

1

( 1)

a

i

ai

N

a N

N

é ù

ê ú

=

é ù

ë û

ë û

G

G

G

G G

G

G

G

O

O

4-19

where the sub-matrices are defined in Sections 2.2.1 and 2.2.2. Furthermore, since the deflections in the 

local, rotated coordinate system are considered small, the shape functions derived for the linearly elastic 

case are used (Section 2.2.3). Using these shape functions, the local strain and displacements in the rotated 

coordinate system are related by the equation

ε Bq ˆ = ˆ 4-20

where B is defined as 

B = GN . 4-21

The principle of virtual work of the element can be written as:

( )0 Int Ext W W- =

4-22

and the internal work can be written as the internal nodal forces multiplied by the nodal virtual 

displacements, or the integral of the strain energy density over the volume of the element:

ˆ ˆ ˆ ( ) Int T Int

V dd d W = = qf q W dV 4-23

Where V is the volume of the element, and W ( ) qˆ is the strain energy density of the element resulting 

from a virtual displacement. Since the deformations are small relative to the rotated coordinate system, the 

internal virtual work can be rewritten as

ˆ ˆ Int T T

V d d W = q B σdV .

4-24

Assuming that external forces only occur as nodal forces and moments, the external virtual work of the 

element becomes

ˆ ˆ Ext T Ext d d W = q f . 4-25
NASA/CR—2012-217606 57

Finally, using Equation 4-22 and noting the fact that the virtual displacements are arbitrary, the resulting 

equilibrium equation is

ˆ ˆ T Ext

V

Bσ f dV= .

4-26

Now, the local internal nodal forces are

ˆ ˆ Int T

V

f B σ= dV 4-27

with the local stiffness matrix being given by

ˆ T

V

k B DB = dV 4-28

where D is an assembly of the adhesive and adherend matrices defined in Sections 2.2.1 and 2.2.2. Note 

that for linear elastic materials, the integration can be carried out analytically in the area, resulting in the 

relation given in Equation 2-47.

4.2.1.4 Global Coordinate System

Now we seek to find the residual and stiffness matrix in the global coordinate system. Since the internal 

work is not dependant on the frame of reference, one can write

ˆ ˆ T Int T Int d d qf qf = 4-29

where the nodal virtual displacements in the global frame are related to those in the local rotated coordinate 

frame through the equation:

ˆ rig dd d q q Tq = + 4-30

making Equation 4-29

ˆ ( ) T Int T T T Int

rig d dd qf q q Tf = - .

4-31

Since rigid body motion does not result in the generation of internal forces, 

ˆ 0 T T Int q Tf rig = 4-32

and Equation 4-31 becomes 

σ
NASA/CR—2012-217606 58

T Int T T Int ˆ d d qf qTf = . 4-33

With the virtual displacements being arbitrary, the internal nodal force vector in the global coordinate 

system becomes

Int T Int ˆ f Tf = . 4-34

To find the global tangent stiffness matrix, differentiation of Equation 4-34 gives

ˆ ˆ Int T Int T Int dd d f Tf T f = + . 4-35

The second term in the above equation becomes

ˆ ˆ ˆ ˆ T Int T T T

rig T f T k q T kT q T kT q d ddd == - .

4-36

The last term vanishes because, as before, displacements resulting in rigid body translation and rotation do 

not generate any internal force. The first term on the right side of Equation 4-35 is somewhat more difficult 

to obtain. The difficulty lies in the fact that T contains sines and cosines of , which in turn contain 1 q

and l . However, Crisfield [48] provides an approximation, which assumes that the extension l is small.

Based on this assumption, the first term in Equation 4-35 can be rewritten as

, , ˆ ˆ T Int T Int = f d Tf T f q q .

4-37

where

1,

,

1,

é ù

ê ú = ê ú

ê ú

ë û

T

T

T

O

4-38

and

1,

0

0

0 00

s c

c s

é ù

ê ú =- - ê ú

ê ú

ë û

T

.

4-39

Similarly, 

, 1, f f = é ù 0 0 q q ë û K 4-40

and
NASA/CR—2012-217606 59

1, [ ]

1

1

s c sc 0 0

l q =--

.

4-41

Combining all of these equations, the global tangent stiffness matrix can be written as a combination of the 

material stiffness, k mat , and the geometric stiffness, kgeo , in the relation

kk k = + mat geo

4-42

where

, , ˆ T Int k Tf

geo = q

4-43

and

ˆ T k T kT mat = .

4-44

Both the geometric and material stiffness matrices are functions of the nodal displacements, making the 

system of equilibrium equations nonlinear. The Newton-Raphson method can be utilized to find the 

solution. It was already noted that one of the benefits of this method is that the formulation of the linear 

element can be utilized. Another major advantage of this method lies in the fact that the local rotational 

frame stiffness and internal force vectors are not functions of the nodal displacements. Since numerical 

integration occurs while finding these vectors/matrices, the integration must only be carried out once during 

the analysis. This saves a considerable amount of computational time, especially for an element like the joint 

element, which requires more refined integration for the higher order shape functions.

4.2.2 Material Nonlinearities

Since modern polymeric adhesives often display highly nonlinear material behavior, it was necessary to 

include material nonlinearities in the joint element to estimate joint strengths more correctly. A particularly 

simple nonlinear elastic stress law was chosen:

σ σε ˆ ˆˆ = ( ) 4-45

where the stress is some general function of the strain. The only major change from the previous corotational formulation is that Equation 4-28 becomes

ˆ ( )ˆ T

V k B DqB = dV 4-46

where the local stiffness matrix in the rotated coordinate system is now a function of the local 

displacements.

Although it would be more correct to use an incremental flow type plasticity formulation that 

distinguishes loading and unloading stiffness, the simple nonlinear elastic relation, which assumes no 

permanent plastic strain, was chosen for several reasons. While this was chosen for simplicity sake, this 

decision can also be justified. The joint element is meant to be a design tool to give general approximations, 

so it is not expected that such a tool will be used in situations requiring unloading capabilities. Additionally, 

the nature of adhesively bonded joints is such that the high stresses occur in concentrated form at the joint 
NASA/CR—2012-217606 60

edges. Since the failing adhesive domain is eliminated in the iteration process (to be described later) the 

assumption of a nonlinear elastic type stress-strain law suffices for this modeling process since potential 

regions of “unloading” are minimal and contained in the regions which are eliminated. Thus, this 

assumption does lead to a meaningful rendition of the joint physics, yet facilitating an efficient (in the 

computational sense) solution strategy.

One other aspect worthy of discussion is the integration requirements for the nonlinear material 

formulation. When the adherends have a nonlinear stress-strain relationship, Equations 4-46 and 4-27 must 

be integrated over xˆ and zˆ at each Newton-Raphson iteration to allow a general stress-strain relationship.

This causes a considerable increase in computational time. However, there are some cases when this is not 

necessary. If only the adhesive layers have a nonlinear stress-strain relation and adhesive model 1 (Section 

2.2.2.1) is used, integration over zˆ can be avoided because the stress is constant through the thickness of 

the adhesive layer. Additionally, if the functions for the nonlinear stress-strain relations are known (and 

simple enough), integration over zˆ can be accomplished analytically. However, this would mean that the 

formulation is only good for that specific stress-strain relation, and cannot be extended to other general 

relations.

4.2.3 Crack Growth

When some user defined failure criterion is reached in some part of the adhesive layer, that portion of 

the adhesive is considered “failed” and can carry no load and has no stiffness. Setting the stress and stiffness 

of that portion of the adhesive to zero is an easy way to model the failure of the adhesive, but the shape 

functions for the joint element were not originally calculated based on a joint with failed adhesive, and 

cannot accurately model this new situation. Therefore, as with more traditional shape function prescribed 

finite elements, more elements are required to accurately find the solution. In the case of failed adhesive, a 

great number of elements may be needed, as will be illustrated in Section 4.3.3.

In order to increase the accuracy of the joint element after adhesive failure and crack growth, a method 

of removing the adhesive and adapting the mesh to the crack was devised. Since the joint element is meant 

to be used as a user defined element in a larger global assembly in commercially available finite element 

software, the mesh change would have to be strictly internal to the element so that the surrounding model 

does not have to change. Therefore, a sub-assembly method was devised to handle adhesive failure

(Figure 4-4) and is outlined Figure 4-5.

Figure 4-4. Diagram showing a) an uncracked joint element, b) a partially cracked element,

and c) a fully cracked joint element.

1

1

2

3

1

2

(a) Uncracked (b) Partially Cracked (c) Fully Cracked

lcrack

inner nodes

outer nodes
NASA/CR—2012-217606 61

Figure 4-5. Flow chart showing how cracked element Sub-Assembly is incorporated

into joint element solution procedure.

First, when failure in the adhesive is detected, the element is replaced by a sub-assembly with three 

elements as shown in Figure 4-4. The length of the crack determines the lengths of the sub-assembly 

elements. Within a Newton-Raphson type solver, the nodal displacements are prescribed (guessed) and the 

stiffness and internal force vector for the element are calculated. These vectors/matrices for all of the 

elements in the assembly are assembled, boundary conditions and loads are applied, and the residual (error 

of the initial nodal displacement guess) is calculated. If the residual isn’t within some tolerated state, a new 

nodal displacement “guess” is calculated based on the previous displacement, residual, and stiffness values 

and the whole cycle repeats.

In the case of a joint element with a crack, only the outer nodal displacements are prescribed since the 

global finite element assembly isn’t aware of the existence of the sub-assembly and the inner nodes.

Therefore, the sub-assembly becomes a nonlinear model within another nonlinear model and must be solved 

with its own Newton-Raphson type solution procedure. The prescribed nodal displacements of the outer 

nodes become the boundary conditions for the sub-assembly, and the whole system is solved using a 

nonlinear solver. When the desired error tolerance is reached, a stiffness matrix and internal force vector for 

the sub-assembly has been calculated. However, these quantities still have the inner degrees of freedom 

contained within. The force vector and stiffness matrix are then reduced using the Guyan Reduction Method

[49–51]. Once the internal degrees of freedom are removed, the stiffness matrix and force vector can be 

considered to be that of the equivalent joint element, and can be passed on to the global assembly.

After the global system is solved, there is a check to see if the crack has grown, or if new adhesive 

failure has been detected. If this is the case, the sub-assembly is adjusted by changing the lengths of the subassembly elements, and the global system is re-solved. This is done until no new adhesive failure occurs and 

Newton-Raphson

Increment load

yes

q j Outer node BCs

KSub , Int FSub

Adhesive Failure?

no

yes

create sub-assembly or adjust existing sub-assembly (lcrack)

Element

Sub-Assembly

Qj

q j

Sufficient accuracy?

Assemble and .

Int Fj K j

no

yes

k j , Int

j f

find 

Qj+1 from 

Kj and 

Rj

: j = j + 1

for each element

Apply Loads, BCs

Find . R j

Element cracked?

no

calculate

Sufficient accuracy?

Assemble and .

no

yes

,

calculate

find Qj+1 from Kj and Rj

Element

for each element

Apply Loads, BCs

Find . 

Newton-Raphson

Guyan

Reduction
NASA/CR—2012-217606 62

the crack is in equilibrium. A crack scaling constant, C5 , has been introduced to speed up or slow down 

crack growth as needed, and is used in the equation

5 ( ) cur cur cur prev

crack crack crack crack l l Cl l =- - 4-47

where prev

crack l is the previous crack length (prior to the global Newton-Raphson procedure) and curcrack l is the 

current crack length. Setting 5 C < 0 causes the crack to grow further than detected, and is useful when 

multiple iterations are needed to find crack equilibrium. Setting 5 C > 0 causes the crack to grow less than 

detected, and is necessary when crack overshoot is a concern. 

The advantage of this method is that fewer elements are needed in order to accurately capture crack 

growth. One can use the minimum elements needed to accurately capture the material and geometric 

nonlinear effects without crack growth being a factor. This can mean dramatically reducing the elements 

required, especially when there is little material nonlinearity and strains in the joint are small.

One of the major disadvantages of this method is the increased computational time. A local nonlinear 

problem must be solved within each iteration of the global nonlinear problem. Although the local nonlinear 

problem is always limited to three elements, it can significantly increase the runtime. Furthermore, the 

global load increment is repeated if the crack grows and the sub-assemblies need to be created or re-meshed.

Although the crack scaling parameter can significantly help in limiting the iterations needed to find crack

equilibrium this process can still be costly. However, the costs can be justified if joint strength prediction is 

of concern. Joint strength has been identified as a controlling factor in the ultimate load bearing capacity of 

many bonded structures.

4.2.4 Adaptive Shape Functions

When the adhesive in a joint element has a nonlinear stress-strain relationship as that shown in

Figure 4-6a, the shape functions used earlier may no longer be well-suited because they were derived for the 

linear-elastic case. For highly nonlinear materials, many elements may be required to accurately capture the 

joint behavior.

As the adhesive softens throughout the loading, the tangent modulus of the adhesive changes across the 

joint, as shown in Figure 4-6b. This resembles the functionally graded adhesive joints discussed in detail in 

Chapter 3. Using the non-constant coefficient ordinary differential solution procedure developed in Section 

3.2.1.1, the shape functions for the joint can be adapted to the adhesive softening at any point in the loading, 

which may decrease the number of elements needed.

Figure 4-6. For a single lap joint with (a) nonlinear material properties, the (b) tangent modulus of the

adhesive layer becomes similar to a functionally graded adhesive in advanced stages of loading.

0

2

4

6

0 3 6 9 12

Tangent modulus, GPa

Distance, x, along adhesive, mm

E

G

0

30

60

90

0 0.02 0.04 0.06 0.08

Stress, MPa

(a) Strain, mm/mm (b)
NASA/CR—2012-217606 63

To this end, adaptive shape functions were applied to the joint element. After each load increment, the 

shape functions for the next increment are calculated for a joint with a functionally graded adhesive. The 

function for the adhesive modulus, E(x), is based on the tangent stiffness of the adhesive in the prior load 

increment. In this way, the shape functions adapt to the softening of the adhesive and can represent the 

behavior of the joint with very few elements (see Section 4.3.4).

This has only been implemented for adhesive model 1 (Section 2.2.2.1) since the stress, and hence the 

tangent modulus, is a function of x only. However, a similar approach could be conceived for the other 

models by basing the adhesive modulus function on the adhesive centerline values (za=0). Furthermore, a 

similar method could be applied to the adherends to improve elemental convergence when the adherends are 

highly nonlinear.

4.2.5 Adhesive Model Characterization

One of the most important inputs for determining the strength of a joint is the characterization of the 

adhesive constitutive response. There have been many methods of characterizing the adhesive material, but 

two have emerged as the most common: bulk adhesive tensile test and fracture mechanics characterization 

tests (DCB, ENF, etc.). Therefore, the following sections outline methods of using both bulk adhesive 

tensile test data and fracture mechanics inputs to characterize the joint. Ultimately, the test data available 

and personal preferences of the user will decide which route to take.

4.2.5.1 Bulk Adhesive Tensile Characterization

One common way of characterizing adhesive materials is by performing tensile tests on bulk adhesive 

specimens, such as those depicted in Figure 4-7. The following section will outline an approximate method 

for modeling the adhesive based on such adhesive characteristic data, and will discuss the formulation and 

underlying assumptions involved.

Figure 4-7. Adhesive may be characterized by (a) experimental bulk adhesive tensile tests, then

(b) fitting a curve to the stress-strain plot.

Experimental Data

Curve Fit

mm

Axial Strain

( ) Bulk Bulk s e = f

(a) (b) Axial Strain Axial Stress
NASA/CR—2012-217606 64

Figure 4-8. Assuming that the adhesive is perfectly bonded to the adherends,

the adhesive can be considered a constrained body under triaxial stress.

If the adhesive is much deeper than thick (bai << ηai), it can be considered to be in a state of plane strain 

in the y direction (Figure 4-8), and the stress-strain relation shown in Equation 2-29 can be applied.

Furthermore, if we assume that the adhesive is perfectly bonded to the adherends and that the adherends are 

much stiffer than the adhesive (Eai << Ei), then it can be argued that the extensional strain in the adhesive is 

much smaller than the peel and shear components ( ˆ ˆˆ , xai ai ai e eg << ) which is the root of the assumption:

ˆ 0 xai . 4-48

This assumption gives rise to the common practice in adhesive joint analysis of ignoring the extensional 

stress and strain in the formulation. With the extensional strains being relatively small, the strain energy of 

the adhesive layer is virtually unaffected and does not necessarily need to be included. Although the 

extensional strain is negligible, the extensional stress (in both the x- and y-directions) is not insignificant, 

placing the adhesive in a state of triaxial stress [52]. Equation 2-29 can be used to write the extensional 

stress in terms of the peel strain:

4 ˆ ˆ s ne xai ai ai = C 4-49

which can be written in terms of the peel stress:

ˆ ˆ 1

ai

xai ai

ai

s s =

.

4-50

The same relation is true for the extensional stress in the y-direction, ˆ yai . This can be used to find the 

extensional stress without necessarily including it into the formulation. Furthermore, the peel stress becomes 

a function of the peel strain only:

xˆ

ˆai z

yˆ

ηai

bai

l
NASA/CR—2012-217606 65

1 ˆ ˆ . (1 2 )(1 )

ai

ai ai ai

ai ai

s e E n n = - +

4-51

This shows that the effective “resistance” to deformation in the z-direction is amplified by a factor that 

depends on Poisson’s ratio. Although this relation is intended for linear elasticity, the relation was assumed 

to hold for the nonlinear stress-strain relation as well. Therefore, the stress-strain relation was redefined as:

1 ˆ ( ) ˆ (1 2 )(1 )

ai

ai ai

ai ai

s e f n n = - +

4-52

which effectively increases the adhesive modulus. This relation can be used with adhesive models 1 and 2 

(Sections 2.2.2.1 and 2.2.2.2) to make them more accurate and to include the effects of having a triaxial 

stress state.

A Von Mises failure criterion was chosen for this particular formulation, although the same formulation 

could easily be altered for a different criterion [53]. Applying the notation for the adhesive layer, assuming 

the shear stresses in the xy and yz planes to be negligible, and using Equation 4-50, the Von Mises 

equivalent stress in terms of the shear and peel stress components are 

2

2 1 2 2 2 3

1

ai

vm ai ai

ai

s st æ ö = + ç ÷

è ø

.

4-53

Although the Von Mises equivalent stress is normally used to find the yield stress, in this case it will be 

assumed to hold for the entire nonlinear adhesive stress/strain response. Therefore, the Von Mises 

equivalent stress for a certain adhesive will be a nonlinear function of the adhesive strain found using bulk 

adhesive tensile tests (Figure 4-7):

( ) vm Bulk Bulk ss e = = f

.

4-54

To find the nonlinear curves approximating the peel and shear stress in the adherend, one more relation must 

be defined. A new variable will be introduced, i , which represents the ratio of peel to shear stress for 

adhesive layer i of a particular joint configuration:

ˆ

ˆ

ai

i

ai

= 4-55

allowing the shear stress to be defined as a function of the bulk stress:
NASA/CR—2012-217606 66

( )

2

2

2 1 2 2

1

ˆ

3 ai

ai

Bulk

ai

i

=

+

.

4-56

The method of finding the strain was a bit more arbitrary. Others have done this by utilizing a Von 

Mises strain criterion or similar methods [6,54,55]. For the current formulation, it was assumed that the bulk 

adhesive tensile specimen strain and the adhesive layer strains were linearly related to each other through 

the equations

6 ˆg e ai Bulk = C 4-57

and

7 ˆe e ai Bulk = C 4-58

where the constants C6 and C7 are found such that the initial slopes of the shear and peel stress-strain curves 

become the normal and shear modulus respectively.

For an actual joint, the ratio of peel to shear stress, i , not only varies across the joint, but changes 

during loading due to nonlinear geometric effects and nonlinear material effects. Therefore, this value will 

in actuality be a function of the joint geometry, loading, materials, and location within the adhesive in 

question. However, to simplify the determination of this value, it is proposed that one assume that the ratio 

of peel to shear doesn’t change significantly during the loading event and that only the stress at ends of the 

joint where the stress concentrations reside is important. The correctness of this first assumption will be 

tested later (Figure 4-33). Therefore, this value can be approximated by taking the ratio of the maximum 

peel to shear stress of the linearly elastic case as illustrated in Figure 4-9.

Figure 4-9. The peel to shear ratio for adhesive i can be approximated by dividing the

maximum peel stress by the maximum shear stress for the linear elastic adhesive case.

0 0.25 0.5 0.75 1

Stress

x/l

ˆaiMAX

ˆ

aiMAX

aiMAX

i

aiMAX
NASA/CR—2012-217606 67

For balanced joints with the same adherend materials and geometries, the maximum occurs on both 

ends of the adhesive and is identical on either end. However, for unbalanced joints, the stress concentrations 

at the ends of the adhesive can be of unequal magnitude. Finding the peel to shear ratio based on the higher 

and lower of the two stress concentrations can provide a good upper and lower bound to the nonlinear 

solution (see Section 4.3.4).

In order to approximate the Von Mises failure criterion for uncoupled shear and peel, an uncoupled 

strain-based criterion was chosen that simply considered the adhesive failed when

ˆ 1 ai

c

=

or

ˆ 1 ai

c

=

.

4-59

where c and c are critical peel and shear strain values. These values are found by applying Equations 

4-57 and 4-58 to the maximum strain of the bulk adhesive tensile test data.

Though it might seem unusual to use a strain-based criterion to approximate the Von Mises stress, it 

should be kept in mind that a Von Mises yield criterion was already applied to get from the bulk adhesive 

tensile test data to the peel and shear stress-strain relations. If the peel to shear ratio, i , was chosen 

correctly, both the shear and peel components should be close to their respective critical values at the same 

time.

4.2.5.2 Fracture Mechanics Characterization

The joint element model with adhesive model 1 is very similar to the Cohesive Zone Models (CZM)

[44,56–58] and is inherently suited for fracture mechanics-type inputs. One of the main differences between

most mainstream cohesive zone models and the joint element adhesive model lies in the thickness of the 

cohesive zone. Most CZM’s have no thickness, and lie at the interface between continuum elements. Since 

it has no thickness, a traction-separation law rather than a stress-strain law is defined for the CZM. Thus,

cracks in the center of the adhesive layer can be differentiated from cracks at the interface by placing CZM 

elements at different locations within the adhesive, although this is computationally very costly. The joint 

element, on the other hand, resembles a cohesive zone with an explicit thickness. The entire adhesive layer 

is a single cohesive zone, and cracks in the middle of the adhesive are not differentiated from those at the 

interface. The traction-separation law can be transferred approximately to a stress-strain law by dividing the 

separation by the thickness as shown in Figure 4-10.

Figure 4-10. Fracture mechanics properties such as critical stress

and fracture toughness can be used to form an adhesive

stress-strain law for the joint element.

Strain

Stress

Ea

IC

IC

ai

G
NASA/CR—2012-217606 68

For this type of adhesive characterization, the shear and peel responses are isolated and characterized in a 

series of experiments [59]. The peel and shear responses are considered to be uncoupled

(see Equation 2-21), and depend solely on the vertical and horizontal separations of the adherends 

respectively (Equations 2-22 and 2-23). Typically, a critical stress and fracture toughness are identified for 

Mode I and Mode II. Since the joint element model does not have continuum elements to represent the 

adhesive, it is recommended that the initial slopes of the stress-strain laws be set to the elastic modulus for 

peel and shear.

Finally, adhesive failure can be defined as occurring when

1.

i j

I II

IC IIC

G G

G G

æ öæ ö

ç ÷ç ÷ + =

è øè ø

4-60

The values of i and j can be chosen based on the preference of the element user.

4.3 Results and Validation

4.3.1 Geometric Nonlinearities

To validate the co-rotational formulation, several example joint and beam configurations were analyzed

using the joint element and compared with 2-D dense mesh finite element solutions with nonlinear 

geometric effects included or an analytical solution where applicable to demonstrate the joint element’s 

ability to capture large rotation situations and to demonstrate how many elements are typically required.

4.3.1.1 Beam Problems

The first example involves a class of beam elements, formulated like a joint element without an 

adhesive layer. The cantilevered beam is clamped at one end and displaced vertically at the other end in 

displacement-controlled loading. The beam is made of aluminum (E=70 GPa) and its shallow depth 

warranted a plane stress Euler-Bernoulli beam formulation for the joint element and 2-D plane stress 

elements for the dense mesh model. The geometric parameters are shown in Figure 4-11. A comparison of 

the load-displacement plots for the linear solution, 2-D dense mesh model, and joint element model is 

shown in Figure 4-12a. It can be seen that the joint element model agrees quite well with the 2-D dense 

mesh model. Additionally, a plot comparing the load-displacement curves for joint element models using 

different numbers of elements is shown in Figure 4-12b to show the number of elements required for a 

converged solution. Even with only two elements, the solution appears to be reasonably converged.

Figure 4-11. Geometric dimensions for aluminum beam problem for validation of the co-rotational formulation.

Δ, F

1 mm

1 mm

20 mm
NASA/CR—2012-217606 69

Figure 4-12. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the cantilever 

beam, along with (b) a convergence study for the joint element model.

Figure 4-13. Load vs displacement plots for a beam loaded in compression with a slight (a) moment and

(b) imperfection angle.

The next example features the same beam, except that it is compressed axially rather than loaded in 

bending to show how buckling loads can be found. Two different methods of applying imperfections are 

utilized. First, a slight moment is applied at the end of the beam (Figure 4-13a) and the load-displacement 

plots for different magnitudes of moments are compared with the well-known analytical solution [60]. 

Second, the beam is given a slight angle to induce buckling, and Figure 4-13b shows the results 

corresponding to different angles. In both examples, 50 joint elements were used and both agree well with 

the analytical solution, as expected.

4.3.1.2 Single Adhesive Layer Joints

The second type of example was single overlap joints, using the joint elements with a single adhesive 

layer and two adherends. These examples pointed more to the accuracy of the joint element formulation.

The unbalanced single lap joint illustrated in Figure 4-14 was pulled in a displacement-controlled 

manner. The adherends were titanium (E=110 GPa) and aluminum (E=70 GPa), with EA 9394 as the 

adhesive layer (E=4 GPa, G=1.79 GPa). As before, the shallow width of the joint required the use of a plane 

stress joint element formulation and the use of 2-D plane stress elements for the dense 2-D finite element 

mesh model. The joint elements used adhesive model 1, and the Euler-Bernoulli adherend model. The joint 

element model had 40 beam elements with one joint element, while 154,000 elements were used for the 2-D 

dense mesh model (Figure 4-15a). A comparison of the load-displacement plots of the different models is 

0

10

20

30

0 5 10

F (N)

Δ (mm)

Linear Solution

Abaqus 2-D

Joint Element

0

10

20

30

0 5 10

1

2

3

4

20

50

Δ (mm)

F (N)

# of Elements

(a) (b)

0

20

40

0 0.05 0.1

10

5

1

0.1

Δ (mm)

F (N)

M (N-mm)

0

20

40

0 0.05 0.1

1

0.5

0.1

0.05

Δ (mm)

F (N)

θ (°)

M F, Δ

θ

F, Δ

Analytical Solution Analytical Solution

(a) (b)
NASA/CR—2012-217606 70

shown in Figure 4-15b. The joint element model was able to replicate the response quite well, even with 

only one beam element rather than 40. The actual joint region requires fewer elements because all of the 

bending takes place outside of the overlap region. The increased flexural rigidity of the overlap region 

causes it to rotate rigidly rather than bend. Therefore, more elements are required outside the overlap regions 

to capture the nonlinear geometric effects of the joint.

The second single adhesive layer example is a layered beam subject to an axial compressive 

displacement. The adherends are aluminum, and the adhesive is EA 9394 (Figure 4-16). The width of the 

beam was 1 mm. Ten elements were found to produce a converged solution for the joint element model, and 

68,000 elements were used for the 2-D dense mesh model. The left end of the top adherend was clamped, 

while the left end of the bottom adherend was free to move vertically only. The right adherends were both 

constrained from rotating and a prescribed axial compressive displacement was applied. The beam was 

given an imperfection angle of 1° to induce buckling.

Figure 4-14. Single lap joint used to validate joint element co-rotational formulation.

Figure 4-15. Comparison of (a) joint element and 2-D dense mesh finite element representation of the joint in

Figure 4-14 and (b) the resulting load vs displacement plot.

Figure 4-16. (a) Geometric parameters and (b) buckled shape of layered beam under compressive loading.

40 mm 60 mm

5 mm F, Δ

1 mm

60 mm

0.5 mm

10 mm

Aluminum

Titanium EA 9394

81 Elements

154,000 Elements

Joint Element

Abaqus 2-D

0

2000

4000

6000

8000

0 1 2

Abaqus 2-D

Linear

1 Joint + 1 Beam Element

1 Joint + 40 Beam Elements

Displacement (mm)

Load (N)

Δ(mm)

F (N)

(a) (b)

F, Δ

0.1 mm

1 mm

20 mm

(a)

(b)

1°
NASA/CR—2012-217606 71

Figure 4-17. Load vs. displacement plot comparing a dense, 2-D finite element

mesh with the joint element using various adhesive and adherend models.

As can be seen in Figure 4-17, the standard model with Euler-Bernoulli adherends and adhesive model 

1 buckled too early. To see whether adherend shear deformations played a significant role, Timoshenko 

adherends with adhesive model 1 was compared. Adherend shear effects played an insignificant role, 

making little difference. A model with Euler-Bernoulli adherends and model 2 adhesive was subsequently 

compared and showed a much better match with the 2-D dense mesh model.

To explain why using adhesive model 2 over 1 made such a difference, one must remember the 

difference between the two models. Model 2 has the full shear strain to displacement relation, while the 

shear strain in model 1 is dependent on the difference in the difference in axial displacements of the 

adherends only. This is the difference between a Winkler Foundation with uncoupled shear and normal 

springs and a bed of shear and normal springs where the normal springs are attached to one another. This 

more accurate term for the shear shows how important shear effects are to buckling of sandwich-type 

structures. To further illustrate this point, the Young’s modulus and shear modulus of the adhesive were 

varied independently for the joint element model with Euler-Bernoulli adherends and adhesive model 1. The 

results of varying these parameters independently are shown in Figure 4-18. Increasing or decreasing the 

Young’s modulus by orders of magnitude has almost no effect on the response, while varying the shear 

modulus alone has a huge impact. Another notable observation is that there is a point where increasing the 

shear modulus doesn’t increase the buckling load. It can be assumed that this upper bound is the buckling 

load of the structure if it were a solid piece of aluminum. On the other side of the spectrum, when the shear 

modulus gets too low, the buckling load will approach two times the load of one adherend alone, as the 

adherends behave as independent beams.

0

200

400

600

800

1000

1200

1400

0 0.1 0.2 0.3 0.4 0.5

2-D Dense Mesh

Adhs. Model 1

Adhs. Model 2

Timoshenko, Adhs. Model 1

Load (N)

Displacement (mm) Δ (mm)

F (N)
NASA/CR—2012-217606 72

Figure 4-18. Load displacement plots for the layered beam under compressive loading showing the effects of 

varying the adhesive (a) Young's modulus and (b) shear modulus illustrating the importance of the shear

modulus on the buckling load.

Figure 4-19. Stress vs. strain approximation for aluminum beam depicted in Figure 4-11.

4.3.2 Material Nonlinearities

The material nonlinearity examples featured in this section only highlight adherend nonlinearity, rather 

than adhesive nonlinearity. Adhesive nonlinearity will be addressed in subsequent sections. This section 

only contains two brief examples, but the limitations and abilities of the joint element in modeling nonlinear 

adherends are shared by beam elements in general, and more in-depth discussion on these limitations and 

how to overcome them are dealt with extensively in literature [61–68].

The first example is the beam from Figure 4-11 where the adherend is modeled as being elasticperfectly plastic with a yield stress of 300 MPa, as featured in Figure 4-19. Large rotations were not 

considered in this analysis to isolate the effects of the material nonlinearity. As can be seen from

Figure 4-20a, the joint element compared quite well with the 2-D dense mesh finite element solution.

However, as shown in Figure 4-20b, it took a good number of joint elements to converge to a solution.

0

200

400

600

800

1000

1200

0 0.1 0.2 0.3 0.4 0.5

10E

E

.1E

0.01E

Load (N)

Displacement (mm)

0

200

400

600

800

1000

1200

0 0.1 0.2 0.3 0.4 0.5

100G

10G

G

.1G

0.01G

Load (N)

Displacement (mm) (a) (b) Δ(mm) F (N)

Δ(mm)

F (N)

0

100

200

300

0 0.01 0.02

σ (MPa)

ε (mm/mm)
NASA/CR—2012-217606 73

Figure 4-20. Comparison of (a) the linear solution, dense 2-D mesh, and joint element solution for the cantilever 

beam with nonlinear materials, along with (b) a convergence study for the joint element model.

Figure 4-21. Plots showing the effect of element size on the load-displacement response of the single lap joint 

featured in Figure 4-14with nonlinear adherends for the (a) joint element, and (b) dense 2-d finite element mesh.

The second example wasn’t quite so positive, and could serve as a good indication of why one would 

want to be careful when modeling joints with nonlinear adherends. The structure under question is the single 

lap joint shown in Figure 4-14, but with elastic-perfectly plastic adherends. The yield stress for the titanium 

was set at 1050 MPa, and the aluminum was at 300 MPa. Large rotations were considered in the analysis, 

and the adhesive was given linear material properties to isolate the effect of nonlinear adherends.

Figure 4-21a shows the load-displacement plot for the joint element model using different numbers of 

elements. As can be seen, none of the models are that far off of each other, but more elements are certainly 

necessary for a converged solution. However, the load-displacement plot did not resemble that of the 

Abaqus 2-D dense mesh model shown in Figure 4-21b. The load predicted by the 2-D dense mesh model 

drops after a peak, whereas the joint element model does not drop, but continues to hold more load. There 

are two explanations for this. First, since the stress-strain relation for the adherend is nonlinear elastic, 

unloading of the adherends is inaccurate. When the adherends first yield, the strain increases dramatically at 

one point (localization) while the rest of the adherend unloads. Since unloading is inaccurately captured in 

the joint element model, it continues to increase in load. The second discrepancy is that beam models still 

have the assumption that the displacement and strain vary linearly in the z-direction. Since this is not the 

case after yielding, the model is inaccurate after initial yielding.

0

2.5

5

0 2.5 5

Force (N)

Displacement (mm)

Linear

Abaqus 2-D

Joint Element

0

2.5

5

0 2.5 5

Force (N)

Displacement (mm)

1

2

3

4

10

20

50

# of Elements

(a) (b) Δ(mm) F (N)

Δ(mm)

F (N)

0

1000

2000

3000

0 0.5 1

1

5

10

40

80

Displacement (mm)

Load (N)

# of Joint 

Elements

0

1000

2000

3000

0 0.5 1

0.5

0.25

0.15

0.1

Displacement (mm)

Load (N)

Element Size 

(mm)

(a) (b) Δ(mm) F (N)Δ(mm)

F (N)
NASA/CR—2012-217606 74

Another observation about the 2-D dense mesh finite element model is that the solution continues to 

change when the element size is reduced. This is due to the stress singularity at the reentrant corners. As the 

element size is decreased, the stress concentration rises and the adherends yield sooner and more 

dramatically. Furthermore, the solution cuts off after the peak for an element size of 0.15 and 0.1 mm. This 

cutoff was due to the commercial FE analysis software, which ends the analysis after the step size has 

become too small. This is also probably due to the stress singularity at the reentrant corners, and illustrates 

some of the potential difficulties of modeling joints.

This last example illustrates why one should avoid using the joint element when failure of the joint is 

dominated by adherend yielding. It also brings out the need of applying some of the measures adopted for 

beam elements to the joint element to better capture the material softening of the adherends.

4.3.3 Crack Growth

To illustrate the benefits of growing a crack by re-meshing rather than just setting the failed adhesive 

stiffness and stress to zero, a bi-layered beam was pulled apart as shown in Figure 4-22a. The beam was

5 mm wide, and the adherends had a stiffness of 100 GPa. The adhesive had a Young’s modulus of 1 GPa, 

and was linear up to failure, which occurred at 500 MPa. The simplistic linear-until-failure adhesive was 

chosen because an analytical solution can be found and because it allows crack growth without material 

nonlinearity, isolating this aspect of the joint element.

Two different models were compared to show the benefits of re-meshing. First, rather than removing 

the adhesive and re-meshing, the stress and stiffness of the adhesive were simply set to zero when the stress 

reached 500 MPa. Second, the failed adhesive was removed and the element was replaced by a subassembly as illustrated in Figure 4-4b. The results of the two models with different ways of handling crack 

growth are shown in Figure 4-23. The benefits of re-meshing are clear. For the first model, the post-peak 

solution oscillates around the analytical solution with the oscillation amplitude reducing for more elements.

The second model with the re-meshing, on the other hand, is extremely close to the analytical solution with 

just a single element. There is some oscillation after the peak, but this is suspected to be caused by crack 

overshoot. This effect, however, disappears entirely with only four elements. This example dramatically 

shows that re-meshing the element to represent crack growth can result in huge elemental savings over 

zeroing the adhesive stiffness.

Figure 4-22 Example of the peeling of a (a) layered beam where the adhesive is modeled as (b) linear until failure to 

demonstrate the joint element crack growth ability.

0

250

500

0 0.25 0.5

Stress (MPa)

Strain (mm/mm)

F, Δ

0.1 mm

5 mm

(a) (b)

40 mm
NASA/CR—2012-217606 75

Figure 4-23. Load displacement plots for the peeling of a layered beam with different numbers of joint elements 

using (a) no re-meshing and (b) re-meshing.

Table 4-1. For an unbalanced joint, the peel to stress ratio is different on each side

and produces a different strength prediction.

Side of 

Adhesive

Stress 

Concentration 

ψ aY (MPa) aY (MPa) Predicted 

Strength (kN)

Left Higher 1.63 29.0 17.1 710

Right Lower 1.04 21.3 20.5 819

4.3.4 Adaptive Shape Functions

To show the effects of using adaptive shape functions, the single lap joint shown in Figure 4-14 with 

linear adherends and a nonlinear adhesive stress-strain relation is modeled with a dense 2-D mesh finite 

element model, the joint element model without adaptive shape functions, and the joint element model with 

adaptive shape functions. The adhesive had a bilinear bulk adhesive tensile test stress-strain relation similar 

to the plot in Figure 4-19 except that the linear properties are that of EA 9394 (E=4 GPa, G=1.5 GPa) and 

the bulk yield stress was 40 MPa. The procedure outlined in Section 4.2.5.1 was followed to find the peel 

and shear yield stress, aY and aY . The adhesive was allowed to yield indefinitely so that no crack would 

form or grow. This is an upper-bound prediction of joint strength according to the global yielding criterion 

proposed by Crocombe [69].

Since the joint was unbalanced, two peel to shear ratios were found; one on each side of the adhesive.

The left side was the side with the greatest magnitude of adhesive stress, while the right side was a bit lower.

Since the adhesive can yield indefinitely, the maximum load will not be reached until both sides of the 

adhesive begin to yield. Therefore, it was expected that the peel to shear ratio of the right side, the last side 

to yield, would result in the most realistic solution. The peel to shear ratios and peel and shear yield stresses,

along with predicted joint strengths, are shown in Table 4-1. A comparison of the load-displacement

response using the peel to shear ratio from the left (high ratio) and the right (low ratio) is shown in

Figure 4-24a. This is expected to provide bounds for the solution.

0

1000

2000

3000

4000

5000

6000

0 0.02 0.04 0.06 0.08 0.1

1

5

10

20

40

Load (N)

Displacement (mm)

# of Joint 

Elements

0

1000

2000

3000

4000

5000

6000

0 0.02 0.04 0.06 0.08 0.1

1

2

3

4

Load (N)

Displacement (mm)

# of Joint 

Elements

(a) (b)

Δ(mm)

F (N)

Δ(mm)

F (N)

analytical solution analytical solution
NASA/CR—2012-217606 76

Figure 4-24. Load displacement plots for the joint depicted in Figure 4-14 with an elastic perfectly plastic adhesive 

with yield stress of 40 MPa. Plot shows (a) results of basing the constitutive properties on the peel to shear ratio of 

the highest stressed side and the lower side, elemental convergence for the (b) 2-D dense mesh model, (c) joint 

element, and (d) joint element with adaptive shape functions.

The load-displacement plot for different sizes of elements using the 2-D dense mesh model is shown in 

Figure 4-24b. As with the nonlinear adherend solution (Figure 4-21), the reentrant corners caused stress 

singularities, which cause the solution to be mesh dependent. However, it appears that for the element sizes 

shown, the joint element predictions provide an upper and lower bound for the 2-D dense mesh solution.

The load-displacement results for the joint element model without adaptive shape functions are shown 

in Figure 4-24c. It takes a large number of elements to converge on a solution. In contrast, the joint element 

model with adaptive shape functions reaches a converged solution with just one element

(Figure 4-24d), revealing the benefit of using adaptive shape functions.

4.3.5 Experimental Comparison

Two different joint types were modeled with the joint element to show some of its the capabilities and 

compare it to experiments. First, double cantilever beam (DCB) specimens tested by Song and Waas [70]

were modeled with the joint element to validate the element with another finite element model and 

experimental data. The specimens were loaded and unloaded several times, which is the ideal situation to 

see the effects of assuming nonlinear elastic adhesive rather than elastic-plastic. Second, two additional 

DCB configurations were compared; one with a brittle adhesive and the other a much more ductile adhesive.

Third, a single lap joint was modeled to show how the joint element compares to experiments in a mixedmode test. The last two examples illustrate how bulk adhesive properties could be used to predict failure in a 

bonded joint.

0

200

400

600

800

1000

0 0.5 1

High Ratio

Low Ratio

Displacement (mm)

Load (N)

0

200

400

600

800

1000

0 0.5 1

1

10

20

40

80

Displacement (mm)

Load (N)

# of Joint 

Elements

0

200

400

600

800

1000

0 0.5 1

0.5

0.25

0.15

0.1

0.05

Displacement (mm)

Load (N)

Element Size 

(mm)

0

200

400

600

800

1000

0 0.5 1

1

2

3

Displacement (mm)

Load (N)

# of Joint 

Elements

Displacement (mm)

Load (N)

(a) (b)

Δ(mm)

F (N)

Δ(mm)

F (N)

(c) (d)

Δ(mm)

F (N)

Δ(mm)

F (N)
NASA/CR—2012-217606 77

4.3.5.1 Song and Waas DCB

The ability of the joint element to predict the force vs. displacement behavior was assessed through 

comparison with experimental results published by Song and Waas [70]. The stable progressive failure 

exhibited in DCB specimens is ideal to show the capability of the joint element to fail progressively. Most 

single lap joints exhibit no stable crack growth; the crack grows almost instantaneously after crack initiation.

Therefore, the DCB example illustrates the full capability of the joint element to grow a crack progressively.

This particular data set was chosen because the authors stated that the failure was fully interlaminar and 

cohesive, which is the type of failure currently modeled by the joint element.

Figure 4-25a shows a DCB specimen and the geometric parameters. The DCB specimens were 

constructed from two different 48 ply unidirectional composite laminates, E719/IM7 and E7T1/G40. The 

adhesive was one of the interlaminar matrix layers. The nonlinear stress-strain relation of the adhesive was 

based on GIc and Ea , shown in Figure 4-26, as was done by Song and Waas [70]. Normally, the mode one 

critical stress is preferred rather than the modulus, but this value was not provided by the authors. The 

geometric and material properties of the specimens are shown in Table 4-2. During the test, the loading was 

halted and the specimen was unloaded several times to measure the crack length within the specimen. The 

DCB specimen was modeled using one joint finite element accompanied by two beam elements on the top 

and bottom of the joint as shown in Figure 4-25b.

Table 4-2. Material properties and geometric parameters for Song/Waas [70] DCB specimens.

Specimen

Adhesive Adhr. Geometric Parameters (mm)

Ea

(GPa)

GIc

(N/m)

E

(GPa) l a b t η

E7T1/G40 4.1 335 116 200 52.6 15.5 4.65 35

E719/IM7 3.3 1130 135 200 35.5 15.1 3.23 6

Figure 4-25. A typical DCB specimen, (a) the geometric parameters and boundary conditions for the DCB 

specimens and (b) the mesh for the joint element DCB model.

l

b a

t

η

Beam Elements

Joint Element

x

F, Δ

(a)

(b)
NASA/CR—2012-217606 78

Figure 4-26. The stress-strain relation of the adhesive was defined based on the

Mode I energy release rate and modulus.

Figure 4-27. Load vs. displacement curves for a) E719/IM7 and b) E7T1/G40 DCB specimens tested by

Song and Waas [70] along with the present joint element model.

Figure 4-27a compares the experimental force-displacement responses for the E719/IM7 DCB 

specimen acquired by Song and Waas [70], along with the response predicted by the present joint element 

model. As can be seen, the present model was quite accurate at predicting the progressive failure of the joint 

based on the given material properties and parameters. The behavior of the E7T1/G40 DCB specimens, 

shown in Figure 4-27b, was not predicted as accurately by the joint element model. The joint element model 

predicted a stiffer elastic response and higher peak load, but the subsequent response is captured quite well. 

This comparison shows that the joint element can be used with fracture properties similar to discrete 

cohesive zone models in predicting the behavior of DCB joints.

Although these two DCB specimens were both loaded and unloaded several times, the non linear elastic 

material model used for the joint element adhesive was still sufficient to capture the overall behavior. This is 

because the stress concentration at the end of the joint causes the plastic zone in the adhesive to remain 

small. Therefore, the advantage of modeling the adhesive as elastic-plastic over non linear elastic is not 

great enough to justify the extra effort and complication. As long as the failure incurred due to crack growth 

is accounted for, the global response will be reasonably represented.

Strain

Stress

Ea

2

3

4

a a aa a

IC

E

E

G se e = - IC

GIC

0

100

200

300

0 2 4 6

Load, N

Vertical end displacement, mm

Joint Element

Song/Waas Experiment 0

60

120

180

0 1 2 3 4

Load, N

Vertical end displacement, mm 

Joint Element

Song/Waas Experiment

(a) (b) F (N)

Δ(mm)

F (N)

Δ(mm)
NASA/CR—2012-217606 79

4.3.5.2 Aluminum DCB Specimens

DCB specimens were manufactured and tested to compare two paste adhesives, EA 9394 and EA 

9309.3NA. The results of these tests were used to assess the ability of the joint finite element to predict the 

difference in performance of two adhesives with very different stress-strain relations based on material 

properties obtained from tensile tests. Solid cylindrical specimens with a 3.175 mm diameter and 3.175 mm 

long test section were machined out of cast adhesive cylinders. The specimens were tested at NASA Glenn, 

and digital image correlation (DIC) techniques (Figure 4-28a) were utilized to obtain the axial strain of the 

specimen at different loads. The strains for several points in the gauge section were averaged, and the stress 

was found by assuming constant stress in the cross section. Figure 4-28b shows characteristic stress-strain

data for the two adhesives, and the equations used to fit the data with a curve. The material parameters used 

in the equations are found in Table 4-3. The functions were chosen because they result in the same curve in 

compression and tension and seem to fit the data adequately.

Table 4-3. Material properties and geometric parameters of DCB specimens.

Specimen

Adhesive Adhr. Geometric Parameters (mm)

Ea

(GPa)

σu

(MPa) εfail υ

E

(GPa) l a b t η 

EA 9394 4.2 49.6 0.016 0.4 69 152.4 63.5 25.4 12.7 0.6

EA 9309.3

NA 2.7 41.3 0.068 0.42 69 152.4 63.5 25.4 12.7 0.55

Figure 4-28. Using a) cylindrical tensile specimens and digital image correlation, the b) stress-strain relations of the 

adhesives EA 9394 and EA 9309.3NA could be defined by curve-fitting experimental data.

(a) (b)

0

10

20

30

40

50

60

0 0.02 0.04 0.06

Stress, MPa

Strain, mm/mm

EA 9394

EA 9309.3NA

Curve Fits

b)

mm

Axial Strain

3

3

a fail u

a aa a

fail

E

E e s se e = -

tanh a au a

u

E ss e =

Stress (MPa)

Strain (mm/mm)
NASA/CR—2012-217606 80

Figure 4-29. Two adherends of an EA 9394 DCB specimen after complete failure. Adhesive found on both 

adherends indicates that failure occurred within the adhesive layer as desired.

Since the joint element considers only failure in the adhesive (cohesive failure), care was taken to 

ensure that the interface between the adherends and adhesive of the DCB specimens would not fail. The 

adherends were 7071 T6 Aluminum, and the surfaces to be bonded were sanded, etched in lye, and anodized 

in a sulfuric acid solution prior to bonding [71]. This treatment was sufficient to produce failures in the 

adhesive layer, as can be seen on the failure surfaces of a post-mortem specimen in

Figure 4-29. The failed specimen has adhesive covering both adherends, which means that the interface was 

not the plane of failure. Glass beads were used to maintain a consistent bond line thickness throughout the 

specimen, and pressure was applied to the specimen during curing. The specimens were allowed to cure for 

seven days at room temperature.

Three DCB specimens for each adhesive were tested on an Instron machine at 0.5 mm/min. All 

specimens failed cohesively like the specimen in Figure 4-29. The load-displacement curves for all six 

specimens are shown in Figure 4-30. The high strain-to-failure of EA 9309.3NA caused these specimens to 

hold over two times load of EA 9394 specimens. The EA 9394 specimens exhibited a load plateau rather 

than dropping in load after adhesive failure was initiated. It is possible that air bubbles in the adhesive 

caused the adhesive to fail prematurely, allowing the joint to not drop in load carrying capacity after failure 

initiated.

The joint element model was able to capture the behavior of the joints rather well. It was found that 

compliance in the experimental load train caused the models to over-predict even the initial linear portion of 

the loading. To compensate for this system compliance, the length of the adherends was increased by 7.5% 

for both DCB specimen types. This number was determined by fitting a linearly elastic model to the initial 

portion of the experimental force/displacement plot. These experiments were very effective in displaying the 

ability of the joint element to predict failure, along with showing how constitutive relations can be applied to 

get progressive failure of a thin adhesive layer.

4.3.5.3 Harris and Adams Single Lap Joint

Lastly, the joint element was compared with experimental data published by Harris and Adams [55] on 

single lap joints. The tests were carried out according to ASTM D1002-72 specifications. The geometric 

parameters are shown in Figure 4-31. The adhesive was MY750 and three different aluminum alloys served 

as the adherends. The only difference between the alloys was the 0.2% proof stress, as shown in 

Table 4-4. The adherends were modeled with an elastic-perfectly plastic stress-strain relation similar to 

the aluminum response shown in Figure 4-19. The adhesive, MY750, was characterized using bulk adhesive 

tensile tests, and the bulk adhesive stress-strain relation is shown in Figure 4-32a.

Adherends Adhesive

5 mm
NASA/CR—2012-217606 81

Figure 4-30. Load vs. displacement curves for DCB specimens with aluminum adherends and EA 9394 and EA 

9309.3NA adhesive, along with the joint finite element model prediction.

Figure 4-31. Geometric parameters for single lap joint tested by Harris and Adams [55].

Table 4-4. Material properties of the single lap joint adherends and adhesive [55].

E (GPa) υ 0.2% Proof Stress (MPa)

MY750 3.44 0.4 -

Aluminum 2L73 70 0.34 430

Aluminum BB2hh 70 0.34 220

Aluminum BB2s 70 0.34 110

0

300

600

900

1200

1500

1800

0 2 4 6 8

Load, N

Vertical end displacement, mm 

EA 9309.3NA

EA 9394

Δ (mm)

F (N)

12.7 mm 63.5 mm

1.6 mm F, Δ

25.4 mm

0.125 mm

Aluminum MY750
NASA/CR—2012-217606 82

Figure 4-32. (a) Stress-strain relation for bulk adhesive, along with peel and shear components for a single

lap joint with ψ=1.4, and (b) corresponding load-displacement plots.

Table 4-5. Experimental and predicted strengths of the single lap joint.

Adherend Experimental Strength (kN) Predicted Strength (kN)

2L73 4.8 ± 0.57 4.46

BB2hh 5.0 ± 0.38 4.52

BB2s 3.5 ± 0.32 5.00

The method outlined in Section 4.2.5.1 was followed to find the adhesive peel and shear stress-strain 

relation. First, the joint was analyzed with linear material properties and small rotations, and the peel to 

shear ratio, was found to be 1.4. Using this value, the Young’s modulus, and the Poisson’s ratio, the bulk 

adhesive tensile data was converted to the peel and shear stress-strain relations shown in

Figure 4-32a. Using this, the joint was modeled with 20 beam elements and one joint element and was 

loaded in a displacement controlled manner until the peak load had been reached. The load-displacement 

plots for the single lap joints with different aluminum alloys are shown in Figure 4-32b, and the results are 

compared with the experimental values found by Harris and Adams [55] in Table 4-5.

The joint with the 2L73 adherends failed without the adherends reaching the yield stress, while the 

BB2hh adherend joint had small amounts of adherend yielding and the BB2s joint was dominated by the 

effects of adherend yielding. Looking back to the single lap joint example of Section 4.3.2, the point was 

made that adherend yielding is not accurately captured by the current formulation of the joint element. As 

expected, the specimen with no signs of adherend yielding, 2L73, had a predicted strength well within the 

experimental error. The specimen with slight yielding, BB2hh, had a predicted strength slightly outside of 

the error range of the experiment. Finally, the BB2s adherend joint, being totally dominated by adherend 

yielding, had a predicted strength much higher than the experimental value. However, if one again uses the 

single lap joint of Section 4.3.2 as an example, one could easily imagine that if adherend plasticity were 

accounted for in a more accurate manner, the predicted peak load would be somewhere around the elbow 

where the slope first drops, around 3 kN. This would bring the prediction much closer to the experimental 

value. Unfortunately, as predicted in Section 4.3.2, the joints with more adherend yielding predict strengths 

increasingly deviating from the experimental value. 

If the elbow is taken to be the point of failure for the BB2s specimens, all three predictions would be 

lower than the experimental strength. There are several possibilities for this discrepancy. The first is that the 

actual joints had quite sizeable fillets at the ends of the adhesive. Although it has been shown that springtype joint models, like the joint element, predict stresses within the bondline similar to those in joints with 

fillets [6], the fillet might reduce the stress enough to increase the strength slightly. Furthermore, the peel to 

0

30

60

90

0 0.01 0.02 0.03 0.04 0.05

Stress (MPa)

Strain

Bulk

ˆa

ˆ

a

0

1000

2000

3000

4000

5000

6000

0 0.2 0.4 0.6

BB2s

BB2HH

2L73

Load (N)

Displacement (mm)

Adherend Type

(a) (b) F (N)

Δ(mm)
NASA/CR—2012-217606 83

Figure 4-33. Peel to shear stress ratio in adhesive layer of the single lap joint as a function of end displacement.

shear ratio, ψ, was only approximated base on the linear elastic joint. However, large rotations and the 

accompanying nonlinearities change the peel to shear ratio, making it a function of the loading. Figure 4-33

shows the value of ψ as a function of the end displacement, Δ. It can be seen that the peel to shear ratio 

drops early on in the loading. Therefore, ψ could be adjusted to yield a more accurate answer.

This comparison showed that, as expected, the joint element is less than accurate with regards to 

adherend material nonlinearity. However the method devised to use bulk adhesive tensile data appears to 

have been successful in approximating the strength of this single lap joint. For most advanced composite 

joints, the adherends display brittle failure, so capturing adherend yielding is of secondary importance.

However, a more precise model could be implemented to consider adherend damage.

0.7

0.9

1.1

1.3

1.5

0 2 4 6 8 10

ψ

Δ (mm)

NASA/CR—2012-217606 85

Chapter 5

Conclusion

In conclusion, a brief summary will be provided for the dissertation and the main points and conclusions 

will be highlighted. Following the summary, recommendations for future work to improve and expand the 

current models will be suggested.

5.1 Brief Summary 

In Chapter 2, a formulation was laid out to find the stiffness of a joint element under linearly elastic 

conditions. The joint element was created by finding the stiffness of a beam-on-elastic foundation type 

model, with any number of beams stacked and attached together by adhesive layers. The adherends were 

considered beams under cylindrical bending, and their displacement fields were reduced down to a function 

of x only. The adhesive was assumed to be linear in the transverse direction and written in terms of the 

adherend displacements. Using such a structural model allowed the governing equations to take the form of 

a system of ODEs, which were subsequently solved in a semi-numerical fashion to yield the stiffness 

matrix. The formulation was done in a general fashion, which allowed several different adherend and 

adhesive models to be inserted. Additionally, since modern composite joints are often very complex with 

tapers and other features, a basic “building block” framework was created to allow the modeling of realistic 

joints.

Several different joint configurations were used to validate the joint element by comparing them with 2-

D dense mesh finite element models. First, a parametric study compared several of the adhesive models for 

different joint geometries. It was found that the simplest model, which was primarily used for subsequent 

analysis, provided good results for thin, long adhesive layers. Additionally, more complex joints with steps 

and tapers were used for validation, illustrating the ability to model such joints. Finally, a joint with 

laminated composite adherends was modeled and found to agree well with a model generated by 

commercial sizing software.

Chapter 3 veered away from the formulation of the joint element, and demonstrated its usefulness in 

design and iteration type analysis. The model was used to show the practicality of using functionally graded 

adhesives in advanced composite joints to spread out the stress concentration and increase the overall joint 

strength. This was done by addressing concerns expressed in the adhesive bonding community, namely 

whether the gains are worth the increased complication, whether the flow of adhesive during manufacturing 

can neutralize or even negate the benefits of grading, and whether the same grading can be used for different 

loading scenarios.

These concerns were addressed by modeling a sample joint configuration with different types of grading 

functions. The analysis showed that adhesive stress reductions up to 17% were possible for the 

configuration studied. Furthermore, a sensitivity study was conducted to show the effect of changing the 

grade shape, which simulated the flow of adhesive during curing. It was shown that for some types of 

grading, the change in shape had only a minimal effect on the beneficial aspects of the grade. However, this 

was not true for all grading functions. Last of all, research revealed that for the configuration studied, the 

optimal grading did not change by changing the loading, nor did the optimum change for different joint 

types. This finding suggests that the same gradation could be used universally for any number of joint types, 

making a case for the mass-production of functionally graded adhesive tapes. All-in-all, the joint element 

was very useful for exploration of functionally graded adhesive joints in an extremely efficient manner.

As an additional demonstration, joints with functionally graded adhesives were manufactured and tested 

to failure. Although the methods of grading were not precise, this demonstrated the potential for such joints.

With test results alone, it was found that the graded adhesive specimens had an average strength increase of 

10% over uniform adhesive specimens.
NASA/CR—2012-217606 86

Finally, in Chapter 4, the linear elastic joint element concept was extended to include large rotations, 

material nonlinearity, and adhesive failure. Large rotations, which occur commonly in adhesively bonded 

joints, were handled through a co-rotational formulation. This formulation separated the displacements into 

rigid body displacement and local deformations about some rotated local coordinate system. The local 

deformations are assumed to be small, so the linear formulation can still be used. Material nonlinearities 

were included into the formulation. However, a nonlinear-elastic model was adopted for simplicity. It was 

shown through examples that, while this model was sufficient for the adhesive layers with high stress 

concentrations and often small plastic zones, it was not accurate for a description of the adherend materials, 

especially in the post-yielded stated. Problems arise with excessive adherend yielding and it is suggested the 

modeling of such joints with the joint element be avoided. On a positive note, the joint load associated with 

adherend yielding can be viewed as an upper limit load for the structural joint, predicted using the joint 

element.

Adhesive failure and crack formation and growth were accounted for through an internal re-meshing 

process. The element with an internal crack was replaced by a sub-assembly with the failed adhesive 

removed. This method added to the computational steps that needed to be taken during the analysis, but 

decreased the number of elements needed to capture progressive failure. To further decrease the number of 

elements required, a method of adaptive shape functions was employed. This used the formulation for the 

functionally graded adhesive element to solve for new shape functions at each load step based on the tangent 

modulus of the adhesive in the prior load step. For highly nonlinear adhesive materials, this method results 

in a dramatic reduction of elements needed.

Finally, methods of finding the nonlinear peel and shear stress-strain curves for the adhesive based on 

experimental procedures were outlined. First, using bulk adhesive tensile data, the response was broken up 

into shear and peel components for a certain joint configuration. This allowed the adhesive to be 

characterized with one test, but limited the shear and peel characterization to be specific to a certain joint 

type, geometry, and materials. Next, the resemblance of the adhesive model to cohesive zone models made 

it a natural candidate for fracture properties such as strength and fracture toughness. Tests are conducted to 

isolate the shear and peel “modes” and characterize them separately. This has the disadvantage of requiring 

more tests, but seems to have fewer assumptions involved. The joint element was constructed such that the 

user can use whichever way he/she feels is the most correct.

To cap off the study, the joint element was compared with several experimental data sets, including 

double cantilever beam and single lap joints. It was found to compare well when little or no adherend 

yielding was involved. The joint element was shown to have practical use in estimating the strength of an 

adhesively bonded joint.

5.2 Future Work

As with any worthy research topic, there is a vast amount of possibilities for further work to either 

extend the concepts or models outlined here or improve upon their current state. Here is a list of some of the 

main areas for future work.

5.2.1 3-D Plate Joint Element

The current joint element formulation used a beam type formulation, assuming plane stress or plane 

strain through the thickness. To make the joint element really useful for vehicle and structural designers, a 

3-D plate or shell type element needs to be constructed. This would require solving a system of PDE’s 

instead of ODE’s, and some way would need to be devised to convert the nodal displacements into edge 

boundary conditions. If this were done, the joint element would be truly ready for distribution and 

widespread use in the engineering community, and its utility would be many-fold.
NASA/CR—2012-217606 87

5.2.2 Joint Element Extensions

The general formulation of the joint element allows it to be extended into other applications. One of the 

easiest and most natural extensions is to create a Sandwich Element that can be used to analyze sandwich 

structures. Similarly, with the ability of having any number of “adherend” layers stacked on top of each 

other, the element could be used to find interlaminar stresses in a composite laminate, and even model 

delamination. If the 3-D plate joint element were created, a composite shell element would be extremely 

useful for designers in predicting the delamination of thin-walled composite structures.

5.2.3 Functionally Graded Adhesives

The study presented here on the usefulness of functionally graded adhesives was only an introduction 

and viability study to show that this concept warrants attention. Further work in this area would involve 

identifying a practical grading method which has the ability of providing very exact higher order gradations 

and demonstrating even greater benefits than the preliminary specimens shown herein. Current advances in 

nanotechnology may provide insights into manufacturing, in an automated manner, functionally graded 

polymer adhesives [72]. Furthermore, the adhesive would need to be characterized fully for the different 

grades. Then, the joint element model could be in conjunction with optimization methods to design the most 

beneficial grading functions. If the concept proves successful, mass-production would be the next logical 

step.

5.2.4 Adherend Yielding

Finally, the current joint element model showed a deficiency in its ability to represent adherend 

yielding. If joint performance beyond adherend yielding is required, then the adherend constitutive model 

would need to be represented through an appropriate incremental flow theory of plasticity approach with 

softening hinges [63], however, at the expense of computational efficiency. 

NASA/CR—2012-217606 89

References

[1] L. Hart-Smith, Adhesive Bonding of Composite Structures—Progress to Date and Some Remaining 

Challenges, Journal of Composites Technology and Research. 24 (2002) 133.

[2] O. Volkersen, Die Nietkraftverteilung in Zugbeanspruchten mit Konstanten Laschenquerschritten, 

Luftfahrtforschung. 15 (1938) 41–47.

[3] Goland, E. Reissner, The stresses in cemented joints, Journal of Applied Mechanics. 11 (1944) 

A17–A27.

[4] L.J. Hart-Smith, Adhesive-bonded single-lap joints, NASACR112236. (1973) 116.

[5] F. Delale, F. Erdogan, M. n. Aydinoglu, Stresses in Adhesively Bonded Joints: A Closed-Form 

Solution, Journal of Composite Materials. 15 (1981) 249 –271.

[6] F. Mortensen, O.T. Thomsen, Analysis of adhesive bonded joints: a unified approach, Composites 

Science and Technology. 62 (2002) 1011–1031.

[7] I.U. Ojalvo, H.L. Eidinoff, Bond Thickness Effects Upon Stresses in Single Lap Adhesive Joints, 

SAE International, Warrendale, PA, 1977.

[8] B. Bednarcyk, Y. Bansal, C. Collier, M.-J. Pindera, J. Zhang, Analysis Tools for Adhesively 

Bonded Composite Joints, Part 1: Higher-Order Theory, AIAA Journal. 44 (2006) 171–180.

[9] J. Zhang, B.A. Bednarcyk, C. Collier, P. Yarrington, Y. Bansal, M.-J. Pindera, Analysis Tools for 

Adhesively Bonded Composite Joints, Part 2: Unified Analytical Theory, AIAA Journal. 44 (2006) 

1709–1719.

[10] M.G. Ostergaard, A.R. Ibbotson, O.L. Roux, A.M. Prior, Virtual testing of aircraft structures, CEAS 

Aeronautical Journal. 1 (2011) 83–103.

[11] M. Eisenberger, D.Z. Yankelevsky, Exact stiffness matrix for beams on elastic foundation, 

Computers & Structures. 21 (1985) 1355–1359.

[12] M. Aydoğan, Stiffness-Matrix Formulation of Beams with Shear Effect on Elastic Foundation, 

Journal of Structural Engineering. 121 (1995) 1265–1270.

[13] P.A. Gustafson, A.M. Waas, A bonded joint finite element for a symmetric double lap joint 

subjected to mechanical and thermal loads, International Journal for Numerical Methods in 

Engineering. 79 (2009) 94–126.

[14] S.E. Stapleton, A. Waas, Macroscopic Finite Element for a Single Lap Joint, in: 

AIAA/ASME/ASCE/AHS/ASC 50th SDM Conference, Palm Springs, California, 2009.

[15] C.L. Dym, I.H. Shames, Solid mechanics: a variational approach, McGraw-Hill, 1973.

[16] M.Y. Tsai, J. Morton, An evaluation of analytical and numerical solutions to the single-lap joint, 

International Journal of Solids and Structures. 31 (1994) 2537–2563.

[17] C.-T. Chen, Linear system theory and design, Oxford University Press, 1999.

[18] C. Moler, C. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, TwentyFive Years Later, SIAM Review. 45 (2003) 3.

[19] ABAQUS User Manual v6.7-1, (2006).

[20] T.P. Lang, P.K. Mallick, Effect of spew geometry on stresses in single lap adhesive joints, 

International Journal of Adhesion and Adhesives. 18 (1998) 167–177.

[21] M.Y. Tsai, D.W. Oplinger, J. Morton, Improved theoretical solutions for adhesive lap joints, 

International Journal of Solids and Structures. 35 (1998) 1163–1185.

[22] J. Zhang, B.A. Bednarcyk, C. Collier, P. Yarrington, Y. Bansal, M.-J. Pindera, 3D Stress Analysis 

of Adhesively Bonded Composite Joints, in: AIAA/ASME/ASCE/AHS/ASC 46th Structures, 

Structural Dynamics and Materials Conference, Austin, Texas, 2005.

[23] L.J. Hart-Smith, D.A. Company, L.R. Center, Analysis and design of advanced composite bonded 

joints, National Aeronautics and Space Administration, 1974.

[24] L.J. Hart-Smith, Adhesive-bonded double-lap joints, NASA CR-112235. (1973).

[25] Q.G. Zeng, C.T. Sun, Novel design of a bonded lap joint, AIAA Journal. 39 (2001) 1991–1996.

[26] U.V.R.S. Turaga, C.T. Sun, Improved Design for Metallic and Composite Single-Lap Joints, 

Journal of Aircraft. 45 (2008) 440–447.
NASA/CR—2012-217606 90

[27] V. Birman, L.W. Byrd, Modeling and Analysis of Functionally Graded Materials and Structures, 

Applied Mechanics Reviews. 60 (2007) 195.

[28] Y. Liu, V. Birman, C. Chen, S. Thomopoulos, G.M. Genin, Mechanisms of Bimaterial Attachment 

at the Interface of Tendon to Bone, J Eng Mater Technol. 133 (2011).

[29] J. Boss, V. Ganesh, C. Lim, Modulus grading versus geometrical grading of composite adherends in 

single-lap bonded joints, Composite Structures. 62 (2003) 113–121.

[30] National Research Council (U. S.). Committee on Structural Adhesives for Aerospace Use, 

Structural adhesives with emphasis on aerospace applications : a report of the ad hoc Committee on 

Structural Adhesives for Aerospace Use, National Materials Advisory Board, National Research 

Council, M. Dekker, New York, 1976.

[31] C. Raphael, Variable-Adhesive Bonded Joints, Applied Polymer Symposium. 3 (1966).

[32] E. Sancaktar, S. Kumar, Selective use of rubber toughening to optimize lap-joint strength, Journal of 

Adhesion Science and Technology. 14 (2000) 1265–1296.

[33] I. Pires, L. Quintino, J.F. Durodola, A. Beevers, Performance of bi-adhesive bonded aluminium lap 

joints, International Journal of Adhesion and Adhesives. 23 (2003) 215–223.

[34] M.D. Fitton, J.G. Broughton, Variable modulus adhesives: an approach to optimised joint 

performance, International Journal of Adhesion and Adhesives. 25 (2005) 329–336.

[35] L.F.M. da Silva, M.J.C.Q. Lopes, Joint strength optimization by the mixed-adhesive technique, 

International Journal of Adhesion and Adhesives. 29 (2009) 509–514.

[36] S. Kumar, P.C. Pandey, Behaviour of Bi-adhesive Joints, Journal of Adhesion Science and 

Technology. 24 1251–1281.

[37] T. Vallée, T. Tannert, J. Murcia-Delso, D.J. Quinn, Influence of stress-reduction methods on the 

strength of adhesively bonded joints composed of orthotropic brittle adherends, International 

Journal of Adhesion and Adhesives. 30 (2010) 583–594.

[38] S. Kumar, Analysis of tubular adhesive joints with a functionally modulus graded bondline 

subjected to axial loads, International Journal of Adhesion and Adhesives. 29 (2009) 785–795.

[39] Aboudi, M.J. Pindera, S.M. Arnold, Elastic response of metal matrix composites with tailored 

microstructures to thermal gradients, International Journal of Solids and Structures. 31 (1994) 

1393–1428.

[40] S. Chapra, R. Canale, Numerical Methods for Engineers, Sixth Edition, 6th ed., McGraw-Hill 

Science/Engineering/Math, 2009.

[41] N.A. Apetre, B.V. Sankar, D.R. Ambur, Low-velocity impact response of sandwich beams with 

functionally graded core, International Journal of Solids and Structures. 43 (2006) 2479–2496.

[42] A.G. Salvi, A.M. Waas, A. Caliskan, Energy absorption and damage propagation in 2D triaxially 

braided carbon fiber composites: effects of in situ matrix properties, Journal of Materials Science. 

43 (2008) 5168–5184.

[43] AF 163-2 Technical Datasheet, (2009).

[44] M.S. Kafkalidis, M.D. Thouless, The effects of geometry and material properties on the fracture of 

single lap-shear joints, International Journal of Solids and Structures. 39 (2002) 4367–4383.

[45] M.D. Banea, L.F.M. da Silva, Adhesively bonded joints in composite materials: An overview, 

Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials Design and 

Applications. 223 (2009) 1 –18.

[46] T. Belytschko, B.J. Hsieh, Non-linear transient finite element analysis with convected co-ordinates, 

International Journal for Numerical Methods in Engineering. 7 (1973) 255–271.

[47] M.A. Crisfield, G.F. Moita, A unified co-rotational framework for solids, shells and beams, 

International Journal of Solids and Structures. 33 (1996) 2969–2992.

[48] M.A. Crisfield, Non-Linear Finite Element Analysis of Solids and Structures, Wiley, 1996.

[49] R.J. Guyan, Reduction of stiffness and mass matrices, AIAA Journal. 3 (1965) 380–380.

[50] S. Chen, H.H. Pan, Guyan reduction, Communications in Applied Numerical Methods. 4 (1988) 

549–556.
NASA/CR—2012-217606 91

[51] T.R. Chandrupatla, A.D. Belegundu, Introduction to Finite Elements in Engineering, 2 Har/Dsk, 

Prentice Hall College Div, 1996.

[52] E.E. Sechler, Elasticity in engineering, Dover Publications, 1968.

[53] R. Raghava, R.M. Caddell, G.S.Y. Yeh, The macroscopic yield behaviour of polymers, Journal of 

Materials Science. 8 (1973) 225–232.

[54] S. Gali, G. Dolev, O. Ishai, An effective stress/strain concept in the mechanical characterization of 

structural adhesive bonding, International Journal of Adhesion and Adhesives. 1 (1981) 135–140.

[55] J.A. Harris, R.A. Adams, Strength prediction of bonded single lap joints by non-linear finite 

element methods, International Journal of Adhesion and Adhesives. 4 (1984) 65–78.

[56] S. Li, M.D. Thouless, A.M. Waas, J.A. Schroeder, P.D. Zavattieri, Use of a cohesive-zone model to 

analyze the fracture of a fiber-reinforced polymer-matrix composite, Composites Science and 

Technology. 65 (2005) 537–549.

[57] D. Xie, A.M. Waas, K.W. Shahwan, J.A. Schroeder, R.G. Boeman, Fracture criterion for kinking 

cracks in a tri-material adhesively bonded joint under mixed mode loading, Engineering Fracture 

Mechanics. 72 (2005) 2487–2504.

[58] P.A. Gustafson, Analytical and Experimental Methods for Adhesively Bonded Joints Subjected to 

High Temperatures., (2008).

[59] P.A. Gustafson, A.M. Waas, The influence of adhesive constitutive parameters in cohesive zone 

finite element models of adhesively bonded joints, International Journal of Solids and Structures. 46 

(2009) 2201–2215.

[60] J.M. Gere, B.J. Goodno, Mechanics of materials, Cengage Learning, 2008.

[61] F. Armero, D. Ehrlich, Numerical modeling of softening hinges in thin Euler–Bernoulli beams, 

Computers & Structures. 84 (2006) 641–656.

[62] F. Armero, D. Ehrlich, An analysis of strain localization and wave propagation in plastic models of 

beams at failure, Computer Methods in Applied Mechanics and Engineering. 193 (2004) 3129–

3171.

[63] D. Ehrlich, F. Armero, Finite element methods for the analysis of softening plastic hinges in beams 

and frames, Computational Mechanics. 35 (2005) 237–264.

[64] B. Vratanar, M. Saje, A consistent equilibrium in a cross-section of an elastic–plastic beam, 

International Journal of Solids and Structures. 36 (1999) 311–337.

[65] M. Jirásek, Analytical and Numerical Solutions for Frames with Softening Hinges, Journal of 

Engineering Mechanics. 123 (1997) 8–14.

[66] J. Wackerfuß, Efficient finite element formulation for the analysis of localized failure in beam 

structures, International Journal for Numerical Methods in Engineering. 73 (2008) 1217–1250.

[67] P.L. Darvall, P.A. Mendis, Elastic-Plastic-Softening Analysis of Plane Frames, Journal of Structural 

Engineering. 111 (1985) 871–888.

[68] Z.P. Bažant, G. Pijaudier-Cabot, J. Pan, Ductility, Snapback, Size Effect, and Redistribution in 

Softening Beams or Frames, Journal of Structural Engineering. 113 (1987) 2348–2364.

[69] A.D. Crocombe, Global yielding as a failure criterion for bonded joints, International Journal of 

Adhesion and Adhesives. 9 (1989) 145–153.

[70] S. Song, A. Waas, Mode I failure of laminated polymeric composites, Engineering Fracture 

Mechanics. 49 (1994) 17–27.

[71] N. Burst, D. Adams, H. Gascoigne, Investigating the Thin-Film Versus Bulk Material Properties of 

Structural Adhesives, Journal of Adhesion. 87 (2011) 72–92.

[72] P. Podsiadlo, M. Qin, M. Cuddihy, J. Zhu, K. Critchley, E. Kheng, et al., Highly Ductile 

Multilayered Films by Layer-by-Layer Assembly of Oppositely Charged Polyurethanes for 

Biomedical Applications†, Langmuir. 25 (2009) 14093–14099.
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-04-2012 

2. REPORT TYPE

Final Contractor Report 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

The Analysis of Adhesively Bonded Advanced Composite Joints Using Joint Finite Elements 

5a. CONTRACT NUMBER

NNX09AM29H 

5b. GRANT NUMBER

NNC3-989 

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Stapleton, Scott, E.; Waas, Anthony, M. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 573754.01.01.02.03.02 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

University of Michigan 

Main Campus 

Ann Arbor, Michigan 48109 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18187 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/CR-2012-217606 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 39 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

This report was submitted as a dissertation in partial fulfillment of the requirements for the degree of Doctor of Philosophy to the University 

of Michigan, Ann Arbor, Michigan, February 2012. Brett Bednarcyk, Structures and Materials Division, NASA Glenn Research Center, 

organization code RXL, 216-433-2012, served as a member of the Ph.D. Committee. Copyright © 2012 by Scott E. Stapleton. 

14. ABSTRACT

The design and sizing of adhesively bonded joints has always been a major bottleneck in the design of composite vehicles. Dense finite element (FE) meshes are required to 

capture the full behavior of a joint numerically, but these dense meshes are impractical in vehicle-scale models where a course mesh is more desirable to make quick 

assessments and comparisons of different joint geometries. Analytical models are often helpful in sizing, but difficulties arise in coupling these models with full-vehicle FE 

models. Therefore, a joint FE was created which can be used within structural FE models to make quick assessments of bonded composite joints. The shape functions of the 

joint FE were found by solving the governing equations for a structural model for a joint. By analytically determining the shape functions of the joint FE, the complex joint 

behavior can be captured with very few elements. This joint FE was modified and used to consider adhesives with functionally graded material properties to reduce the peel 

stress concentrations located near adherend discontinuities. Several practical concerns impede the actual use of such adhesives. These include increased manufacturing 

complications, alterations to the grading due to adhesive flow during manufacturing, and whether changing the loading conditions significantly impact the effectiveness of the 

grading. An analytical study is conducted to address these three concerns. Furthermore, proof-of-concept testing is conducted to show the potential advantages of functionally 

graded adhesives. In this study, grading is achieved by strategically placing glass beads within the adhesive layer at different densities along the joint. Furthermore, the 

capability to model non-linear adhesive constitutive behavior with large rotations was developed, and progressive failure of the adhesive was modeled by re-meshing the joint 

as the adhesive fails. Results predicted using the joint FE was compared with experimental results for various joint configurations, including double cantilever beam and single 

lap joints.

15. SUBJECT TERMS

Bonded joints; Adhesives; Finite element method; Mathematical models; Nonlinearity; Shape functions; Structural design; Bending

models; Crack propagation; Weight reduction; Composite materials 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

104 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18