Scott E. Stapleton and Anthony M. Waas

University of Michigan, Ann Arbor, Michigan 

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Modeling Progressive Failure of Bonded Joints 

Using a Single Joint Finite Element 

NASA/TM—2010-216824

November 2010

AIAA–2010–3100
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 

Brett A. Bednarcyk

Glenn Research Center, Cleveland, Ohio

Modeling Progressive Failure of Bonded Joints 

Using a Single Joint Finite Element 

NASA/TM—2010-216824

November 2010

AIAA–2010–3100

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

51st Structures, Structural Dynamics, and Materials Conference

cosponsored by the AIAA, ASME, ASCE, AHS, and ASC

Orlando, Florida, April 12–15, 2010
Acknowledgments

The authors would like to thank Justin Bail of the Ohio Aerospace Institute for testing the adhesive tensile specimens at the NASA 

Glenn Research Center. 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. Additional fi nancial support was provided by the NASA 

Glenn Research Center through the GSRP Fellowship.

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://gltrs.grc.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 technical management. 
NASA/TM—2010-216824 1 

Modeling Progressive Failure of Bonded Joints 

Using a Single Joint Finite Element 

Scott E. Stapleton and Anthony M. Waas 

University of Michigan 

Ann Arbor, Michigan 48109 

Brett A. Bednarcyk 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

Enhanced finite elements are elements with an embedded analytical solution which can capture 

detailed local fields, enabling more efficient, mesh-independent finite element analysis. In the present 

study, an enhanced finite element is applied to generate a general framework capable of modeling an 

array of joint types. The joint field equations are derived using the principle of minimum potential energy, 

and the resulting solutions for the displacement fields are used to generate shape functions and a stiffness 

matrix for a single joint finite element. This single finite element thus captures the detailed stress and 

strain fields within the bonded joint, but it can function within a broader structural finite element model. 

The costs associated with a fine mesh of the joint can thus be avoided while still obtaining a detailed 

solution for the joint. Additionally, the capability to model non-linear adhesive constitutive behavior has 

been included within the method, and progressive failure of the adhesive can be modeled by using a 

strain-based failure criteria and re-sizing the joint as the adhesive fails. Results of the model compare 

favorably with experimental and finite element results. 

1.0 Introduction 

In the aerospace industry, fiber-reinforced polymer matrix composites (FRPCs) are gaining increasing 

use and attention because of their high strength-to-weight ratios, among other factors. FRPC joints 

perform much better with adhesive bonding rather than bolting or riveting because of their quasi-brittle 

nature (Ref. 1) and the ability of the bond to spread the load over a larger area leading to a lessening of 

the stress concentration (Ref. 2). Therefore, accurate analysis of adhesively bonded joints is becoming 

more critical than ever. 

Adhesive joints have traditionally been analyzed using two methods: analytical models and finite 

element analysis (Ref. 3). Analytical methods have been utilized to extract efficient closed-form solutions 

for adhesive single lap joint stresses. Classical formulas have been introduced by Volkerson (Ref. 4), 

Goland and Reissner (Ref. 5) and Hart-Smith (Ref. 6). More recently, refined analytical studies, carried 

out by Mortensen and Thomsen (Ref. 7) and Delale et al. (Ref. 8), have proven to be quite accurate in 

predicting stresses within adhesive joints. However, analytical methods are often limited by geometric 

assumptions used to obtain a closed form solution and are not as useful to designers for compiling 

vehicle-scale models that may contain multiple joints. Finite element analyses are widely utilized in 

industry, and can be used to assess joints with a wide variety of geometries and loading conditions. 

However, these methods can suffer from mesh-dependence and a lack of efficiency, which is especially 

crippling for initial sizing analysis and full vehicle-scale models (Ref. 9). 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 costs. Such tools can be used to make quick meshindependent assessments of bonded composite joints. Currently, such a capability is lacking, and joint 

assessment is typically performed late in the design cycle when structural changes that can lower the 

weight are much more difficult and expensive. 
NASA/TM—2010-216824 2 

Gustafson and Waas (Ref. 3) have merged analytical and finite element methods to perform efficient, 

mesh-independent finite element analysis of double lap joints to use for initial design and macroscopic 

vehicle modeling. Analytical models were embedded into a single finite element with minimal analyst 

input. Stapleton and Waas (Ref. 10) extended this method to single lap joints (Fig. 1), which are more 

complex due to the eccentricity of the axial load path. The adherends were treated as Euler-Bernoulli 

beams, and the adhesive was modeled as a bed of normal and shear springs. The current study extends the 

method to include a nonlinear constitutive model for the adhesive in conjunction with inputs to cohesive 

zone finite element modeling (Refs. 11 and 12) or a curve-fit to experimental test data. Additionally, a 

strain-based failure criterion is utilized to track damage in the adhesive and the joint element and adjacent 

beam elements are re-sized to account for the failed adhesive. The entire joint can then be replaced by a 

single joint finite element, while the remaining structure (outside the joint) is modeled using standard 

structural elements, for instance beam elements (Fig. 1). As a result, failure in the adhesive is built into 

the joint element. 

2.0 Analytical Formulation 

The joint element with adhesive failure capabilities is an extension of the linear elastic joint element 

derived elsewhere (Stapleton and Waas (Ref. 10)). The main difference is that this model now has the 

capability of handling adhesives with non-linear constitutive behavior and has a failure criteria 

implemented within the model. A flow chart of the joint element implementation is shown in Figure 2. 

The individual steps will be expanded upon in the following sections, but the basic layout consists of 

solving the linear problem first to obtain shape functions, defining a non-linear stress/strain relation for 

the adhesive and linearizing this relation at each load increment, calculating the joint stiffness matrix and 

force vector using the linear shape functions and strain energy for the non-linear case, and solving the 

resulting equations. Furthermore, the adhesive is checked for failure, and in the case of adhesive failure, 

the length of the joint element decreases to account for failed adhesive. 

2.1 Obtain Linear Shape Functions 

To model nonlinear adhesive in a joint element, the shape functions are obtained from the case of a 

joint with a linear adhesive. To do this, the adhesive and adherends were assumed to be linearly elastic. 

Figure 3 illustrates the required geometric parameters and material properties of a single lap joint. The 

width of the joint in the y-direction is b. The subscripts 1, 2, and a denote a variable associated with 

adherend 1, 2, or the adhesive. 

P

P

P P

Beam Elements Beam Elements

Joint 

Element

Figure 1.—The present study replaces a complex single lap joint mesh 

with a single joint element. 
NASA/TM—2010-216824 3 

Figure 2.—Flow chart of the joint element with adhesive progressive failure. 

Figure 3.—Required geometric and material parameters for overlap 

region of a single lap joint. 

Assuming the adherends behave like Euler-Bernoulli beams, the strain energy of the joint, U, is 

written as: 

ܷ௃௢௜௡௧ ൌ ଵ

ଶ ׬ ߪଵ௫߳ଵ௫ ܸ݀ ௏భ ൅ ଵଶ ׬ ߪଶ௫߳ଶ௫ ܸ݀ ௏మ ൅ ଵଶ ׬ ሺߪ௔௭߳௔௭ ൅ ߬௔௫௭ߛ௔௫௭ሻ ܸ݀ ௏ೌ (1) 

where ߪ௜௫ and ߳௜௫ represent the normal stress and strain in material i (1 or 2 for the adherends) in the xdirection, ߪ௔௭ and ߳௔௭ are the normal stresses/strains in the adhesive in the z-direction, ߬௔௫௭ and ߛ௔௫௭

represent the shear stress/strain in the adhesive on the xz-plane, and all integrals are taken over the 

volume, Vi of material i. 

Ea,Ga

l

1 E1

a 2 x

z2

η

t1

t2

E2

x

z1

Adhesive

Adherend 1

Adherend 2
NASA/TM—2010-216824 4 

Figure 4.—Geometric and material parameters for 

overlap region of a single lap joint. 

The two adherends are assumed to behave as beams, and the strain energy of the adherends can be 

found using standard Euler-Bernoulli beam theory. The derivation is not shown here, but the details are 

contained in Stapleton and Waas (Ref. 10). 

It should be noted that many adhesive joints involving bonded FRPCs are made up of thin adherends 

which are relatively long in the y-direction, causing them to behave more like plates in cylindrical 

bending rather than beams. Though not done herein, to model the adherends as wide plates in cylindrical 

bending, it is only a simple matter of replacing the modulus of elasticity, Ei with Ei/(1 – νi

2

), I = 1,2 and 

modeling the adhesive using plane strain, rather than plane stress assumptions. 

A diagram of the notation scheme for the adhesive is shown in Figure 4. It is assumed that the 

displacement varies linearly in the za-direction and that the adhesive and adherend are perfectly bonded at 

the interface. The displacements at the interface will be denoted by the subscript i. 

The adhesive is modeled as a bed of linear normal and shear springs. In terms of the interface 

displacements, the strains in the adhesive are: 

߳௔௭ ൌ ଵ

ఎ ൫ݓ௜ଵሺݔሻ െ ݓ௜ଶሺݔሻ൯ (2a) 

and 

ߛ௔௫௭ ൌ ଵ

ఎ ൫ݑ୧ଵሺݔሻ െ ݑ୧ଶሺݔሻ൯ (2b) 

and the stresses are: 

ߪ௔௭ ൌ ܧ௔߳௔௭ (3a) 

and 

߬௔௫௭ ൌ ܩ௔ߛ௔௫௭. (3b) 

The interface displacements are defined in terms of adherend centerline displacements using EulerBernoulli beam theory: 

ݑ୧ଵሺݔሻ ൌ ݑଵሺݔሻ ൅ ௧భ

ଶ ቀ

ୢ௪భሺ௫ሻ

ୢ௫ ቁ, ݑ୧ଶሺݔሻ ൌ ݑଶሺݔሻ െ ௧మ

ଶ ቀ

ୢ௪మሺ௫ሻ

ୢ௫ ቁ, 

 (4) 

ݓ୧ଵሺݔሻ ൌ ݓଵሺݔሻ, and ݓ୧ଶሺݔሻ ൌ ݓଶሺݔሻ. 

Using the principle of stationarity of potential energy, four fully-coupled governing equilibrium 

differential equations are obtained from the energy expression in Equation (1). Of the four governing 

equations, two correspond to the axial equilibrium, while two equations correspond to the transverse 

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

η

za

x

a

x

u ( ); i1 x

u ( ) i2 x

w ( ) i1 x

w ( ) i2 x

u ( ) a x z, a w ( ) a x, z u a a(x,za) wa(x,za) a
NASA/TM—2010-216824 5 

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: 

ሼݑሽ,୶ ൌ ሾܣሿሼݑሽ (6) 

where 

ሼݑሽ ൌ ሾݑଵ ݑଵ,௫ ݓଵ ݓଵ,௫ ݓଵ,௫௫ ݓଵ,௫௫௫ ݑଶ ݑଶ,௫ ݓଶ ݓଶ,௫ ݓଶ,௫௫ ݓଶ,௫௫௫ሿT (7) 

and 

 (8) 

with 

௕ீೌ ൌ௜ ߣ

ఎா஺೔

, ݅ ൌ 1,2 (9) 

Λ௜ ൌ ீೌ௕

ఎாூ೔

, ݅ ൌ 1,2 (10) 

Ω௜ ൌ ாೌ௕

ఎாூ೔

, ݅ ൌ 1,2 . (11) 

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

solution method. There are 12 eigenvalues of [A]: two real eigenvalues, four complex eigenvalues, and six 

repeating eigenvalues. Therefore, the solution is made up of two exponential terms, four exponential 

terms multiplied by a sine or cosine, and the six 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 Equation (6) 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 in Equation (6) can be written in terms of the matrix exponential, 

ሾ݁஺௫ሿ, as 

[ ] A 





































































































































































0 1 0 0 0 00 0 0 0 0 0

1 0 0 1

2

t

1 1 0 0 1 0 0 1

2

t

2 1 0 0

0 0 0 1 0 00 0 0 0 0 0

0 0 0 0 1 00 0 0 0 0 0

0 0 0 0 0 10 0 0 0 0 0

0 1

2

t

11 1 0 1

4

t

1

2

1 0 0 

1

2

t

11 1 0 1

4

t

11 t2 0

0 0 0 0 0 00 1 0 0 0 0

2 0 0 

1

2

t

1 2 0 0 2 0 0 

1

2

t

2 2 0 0

0 0 0 0 0 00 0 0 1 0 0

0 0 0 0 0 00 0 0 0 1 0

0 0 0 0 0 00 0 0 0 0 1

0 1

2

t

22 2 0 1

4

t

22 t1 0 0 

1

2

t

22 2 0 1

4

t

2

2

2 0
NASA/TM—2010-216824 6 

ሼݑሽ ൌ ሾ݁஺௫ሿሼܥሽ. (12) 

The matrix exponential can be expressed as the infinite series (Ref. 13)

ሾ݁஺௫ሿ ൌ ෍ ௫ೖ

௞! ሾܣሿ௞ ஶ

௞ୀ଴

. (13) 

To obtain faster convergence, a method of scaling and squaring (Ref. 14) 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 ሾ݁஺௫ሿ for k – 1 and k. The 

value of the acceptable error was set at 0.0001. The next step is to solve for the vector of constants, {C}, 

using the boundary conditions. It is at this point that the methods to model a single lap joint, double 

cantilever beam, and bonded doubler diverge. For the current study, a double cantilever beam (DCB) will 

be used to validate the joint element; therefore the boundary conditions will be those of the DCB (Fig. 5). 

Applying prescribed displacements to the joined nodes on the left gives the following system of 

equations: 

ሾܥܤ଴ሿሼݑ଴ሽ ൌ ሼqሽ (14) 

where ሼݑ଴ሽ is ሼݑሽ evaluated at x=0, 

ሼݍሽ ൌ ሾݍଵ ݍଶ ݍଷ ݍସ ݍହ ݍ଺ሿ், (15) 

and 

 . 

(16)

At the right end of the joint, the nodes are not joined to any other elements and are referred to as “free 

nodes” (Fig. 5). At these nodes, the natural boundary conditions obtained by minimizing the strain energy 

(Eq. (6)) are applied. These boundary conditions enforce a traction free boundary at the ends, and there is 

a term for the axial force, transverse force, and moment at the free nodes. 

Figure 5.—Boundary conditions for a DCB type joint element: prescribed nodal 

displacements at the joined nodes (x=0) and stress free at the free nodes (x=l). 

[ ] BC  0





























































100000000000

000010000000

000001000000

001000000000

000000001000

000000000100

q2

q3

q q5 6

q1

q4

free nodes joined nodes
NASA/TM—2010-216824 7 

The boundary conditions are applied in the following equation: 

ሾܥܤ௟ሿሼݑ௟ሽ ൌ ሼ0ሽ (17) 

where ሼݑ௟ሽ is ሼݑሽ evaluated at x=l, {0} is a vector of zeros, and 

(18)

The first and fourth rows of the matrix in Equation (18) correspond to the axial force, the second and fifth 

correspond to the moment, and the third and sixth are the shear. Equation (12) is substituted into 

Equation (18) and Equation (17) and stacked on top of each other to get 

ሾζሿሼܥሽ ൌ ቄ

q

0

ቅ. (19) 

where 

ሾζሿ ൌ ൤

ሾܥܤ଴ሿ

ሾܥܤ௟ሿ

ሾ݁஺଴ሿ

ሾ݁஺௟ሿ

൨. (19) 

Using this relation, one can obtain an expression for {C}: 

ሼܥሽ ൌ ሾߞሿିଵ ቄ

ݍ

0ቅ . (20) 

The expression can be simplified further by removing the zeros from the right hand side along with six 

columns on the right side of ሾߞሿିଵ to get {C} in terms of the nodal degrees of freedom only. This reduced 

matrix will be denoted by ሾߞԢሿ . The constants are now given in terms of the nodal degrees of freedom:

ሼܥሽ ൌ ሾߞԢሿሼݍሽ . (21) 

This relation can be inserted into Equation (12) to get the adherend centerline displacements in terms of 

the nodal degrees of freedom, q1–6: 

ሼݑሽ ൌ ሾܰሿሼݍሽ (21) 

where the shape functions, ሾܰሿ, are defined as: 

ሾܰሿ ൌ ሾ݁஺௫ሿሾߞԢሿ. (21) 

2.2 Define Adhesive Stress/Strain Relation 

With the shape functions determined for a joint with a linear adhesive, the non-linear constitutive 

stress/strain relations of the adhesive need to be defined. This function can be based on constitutive 

relations, or fracture mechanics. A stress/strain relationship based on constitutive relations would take a 

tensile stress/strain curve with a fitting function to the data. A relation based on fracture mechanics would 

use a function which has its maximum stress at the mode one critical stress (ߪூ௖) and the area under the 

[ ] BC  l

















































































0 1 0 00 0 0 00 0 0 0

0 0 0 10 0 0 00 0 0 0

1

2

t

11 0 

1

2

t

11 0 0 1

4

t

1

2

1 0 -1 0 1

4

t

11 t2 0 0

0 0 0 00 0 1 00 0 0 0

0 0 0 00 0 0 00 0 1 0

1

2

t

22 0 

1

2

t

22 0 0 1

4

t

22 t1 0 00 1

4

t

2

2

2 0 -1
NASA/TM—2010-216824 8 

curve would be the critical strain energy release rate of the adhesive (GIC). Regardless of the method used 

to define the relation, the adhesive stress, ߪ ௔can be written as a function of the adhesive strain, ߳௔: 

ߪ ௔ൌ ݃ሺ߳௔ሻ. (22) 

Although this paper refers to the normal stress and strain in the adhesive only, the same derivation holds 

for the shear strain/stress relation. For simplicity, the subscript a will be dropped from the stress and 

strain symbols. 

It should also be noted that this stress/strain relation assumes no permanent plasticity, but resembles 

non-linear elasticity. Although it is acknowledged that this relation is not representative of the actual 

behavior of most adhesives, this relation was found to be sufficient for the purpose of estimating failure 

loads to aid in initial vehicle design and scaling. 

2.3 Linearize Stress/Strain Relation 

To simplify calculations and avoid the need for a non-linear solver, the stress/strain relation of the 

adhesive is linearized about the previous strain increment. The Taylor series expansion of the stress about 

a strain at load step n, ߳௡ is: 

ߪሺ߳ሻ ൌ ݃ሺ߳௡ሻ ൅ ௗ௚ሺఢ೙ሻ

ௗఢ ሺ߳െ߳௡ሻ ൅ ܪܱܶ) 23 (

where HOT represents higher order terms. To linearize, the higher order terms are ignored. 

2.4 Calculate the Adhesive Strain Energy 

The adhesive normal strain energy ܷ௡ାଵ at the next load step, n+1, is found as the strain energy from 

the previous increment plus the integral of the stress as a function of strain from the previous increments 

to the current increment: 

ܷ௡ାଵ ൌ ׬ ׬ ߪሺ߳ሻܸ݀߳݀ ఢ೙శభ

ఢ ௏ ೙ ೌ ൅ ܷ௡ . (24) 

Carrying out the inner integral gives: 

ܷ௡ାଵ ൌ න ቀ

ଵ

ଶ

ௗ௚ሺఢ೙ሻ

ௗఢ ሺ߳௡ାଵሻଶ ൅ ቀ݃ሺ߳௡ሻ െ ௗ௚ሺఢ೙ሻௗఢ

߳௡ቁ ߳௡ାଵ െ ݃ሺ߳௡ሻ߳௡ െ ଵ

ଶ

ௗ௚ሺఢ೙ሻ

ௗఢ ሺ߳௡ሻଶቁ ௏ೌ

ܸ݀ ൅ ܷ௡. (25) 

2.5 Perform Rayleigh/Ritz Using Linear Adhesive Shape Functions 

To obtain the stiffness and force matrices for the joint, the shape functions derived for the linear 

adhesive case (Eq. (5)) are used. Using Equation (2) and (4) and the shape functions derived for the linear 

adhesive case, the strain in the adhesive is found in terms of the nodal displacements, q1–6. The strain in 

the adhesive at the current, n+1, increment is written as a function of x and q1–6: 

߳௡ାଵ ൌ ݂ሺݍଵି଺, ݔሻ (26) 

while the displacements from the previous increment are used to define the adhesive strain at the previous 

increment, ߳௡, as a function of x only: 

߳௡ ൌ ݂ሺݔሻ. (27) 

The energy is then minimized, which yields the i

th, jth component of the contribution to local joint 

stiffness matrix from the adhesive, k, to be: 
NASA/TM—2010-216824 9 

݇௜,௝ ൌ డమ

డ௤೔ డ௤ೕ

න ଵ

ଶ

ௗ௚ሺఢ೙ሻ

ௗఢ ሺ߳௡ାଵሻଶ

௏ೌ

ܸ݀ (28) 

and the i

th component contribution of the adhesive to the local joint force vector, f, to be 

݂௜ ൌ డమ

డ௤೔ න ቀെ݃ሺ߳௡ሻ ൅ ௗ௚ሺఢ೙ሻ

ௗఢ ሺ߳௡ሻቁ ߳௡

௏ೌ

ܸ݀ (29) 

The contribution of the adhesive to the local force vector and stiffness matrix can be added to the 

contributions of the two adherends. Since only the adhesive is modeled as nonlinear, the contributions of 

the adherends to the stiffness matrix are the same as in the linear case. The local stiffness matrix and local 

force vector can be used to find the local nodal displacements, ሼࢗሽ: 

ሾ݇ሿሼݍሽ ൌ ሼ݂ሽ. (30) 

2.6 Assemble Global Matrices, Apply Loading, and Solve Global Equations 

Once the joint element stiffness matrix and load vector are found, they are assembled with the rest of 

the elements in the model. The loading increment is applied, and the system of linear equations is solved. 

For this particular study, an in-house finite element code was used to assemble and solve the finite 

element global equations. Future work will implement the method into an Abaqus (Ref. 11) user element 

to be integrated in global, vehicle scale models. 

2.7 Check for Adhesive Failure 

Once the nodal displacements are found, Equation (2a) and (2b) are used to find the strain in the 

adhesive as a function of horizontal position, x. Then, a strain-based failure criterion based on the failure 

strain, ߳௙௔௜௟, is used to determine if and how much the adhesive has failed (Fig. 6). 

If failure is detected, the joint element is shortened by the length of the failed adhesive region, and the 

adjacent beam elements are lengthened to compensate (Fig. 7), with no continuity of displacements 

imposed between them. After this step is completed, the steps described in Sections 2.1 to 2.7 are 

repeated until the joint reaches equilibrium. After equilibrium is reached and further failure no longer 

occurs, the load is increased by one increment, and the process is repeated. Through this method, the 

stresses, strains, loads, and displacements for the joint can be found at each load increment as the joint 

deforms non-linearly and progressively fails. 

Figure 6.—Example of modeling failure in the adhesive: the normal 

strain surpasses the failure strain at the right end of the adhesive, 

causing this region to be categorized as “failed”. 

za

x

failed adhesive region

Distance along the adhesive

Strain

Ԗfail 
NASA/TM—2010-216824 10 

Figure 7.—When failure occurs in the adhesive, 

the joint element is shortened and the adjacent 

beam elements are lengthened (DCB geometry 

assumed). 

3.0 Comparison with Published Results 

The ability of the joint finite element to predict the force versus displacement behavior was assessed 

through comparison with experimental and finite element results published by Song and Waas (Ref. 15). 

Results for a double cantilever beam (DCB) were chosen rather than a single lap joint because of the 

simple loading in the adhesive (no shear stress). 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 8(a) 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 non-linear stress/strain relation of the adhesive 

was based on GIC and Ea (Fig. 9). 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 1. 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 8(b). 

Figure 10 compares the experimental force/displacement responses for the E719/IM7 DCB specimen 

acquired by Song and Waas (Ref. 15), along with the response predicted by these authors using an FEA 

model including discrete cohesive zone elements and 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 two models used different adhesive stress/strain relations, which 

probably accounts for the difference in predicted responses. However, Song and Waas (Ref. 15) reported 

as adhesive fails...

as adhesive fails...

Beam Elements

Joint Element
NASA/TM—2010-216824 11 

 Figure 8.—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. 

Figure 9.—The stress/strain relation of the 

adhesive was defined based on the mode 

one energy release rate and modulus. 

TABLE 1.—MATERIAL PROPERTIES AND GEOMETRIC PARAMETERS 

FOR SONG AND WAAS (REF. 15) DCB SPECIMENS 

Specimen Adhesive Adherend Geometric Parameters 

Ea 

(GPa)

GIc 

(N/m)

E 

(GPa)

l 

(mm)

a 

(mm)

b 

(mm)

t 

(mm)

η

(μm)

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 

that their FEA took 1 to 2 hr to run each load loop, while the present analysis took only 30 to 60 s per 

load loop, a speed increase of a factor of 120. This demonstrates the usefulness of the joint finite element 

in performing quick, computationally-efficient analysis of joints, even in progressive failure scenarios. 

The plots in Figure 11 are helpful in visualizing how the stress and strain evolve in the adhesive with 

increasing end displacement. Figure 11(a) shows the load/displacement plot for the E719/IM7 DCB 

specimen. The circles show points during loading at which stress (Fig. 11(b)) and strain (Fig. 11(c)) 

snapshots were taken. The solid circles correspond to solid lines on the stress and strain plots, while the 

hollow circles correspond to the dotted lines on the plots. The stress in the adhesive rises at the left end of 

the adhesive until the maximum stress according to the stress/strain relation (Fig. 9) is reached, at which 

point the stress at the end of the adhesive reduces to zero. At this point, the adhesive starts to fail. The 

strain reaches a critical level, about 1.6 percent strain, before failure begins. It can be seen that the shape 

of the stress/strain relation will have a significant effect on the stress levels at the end of the joint. 

l

b a

t

η

Beam Elements

Joint Element

a)

b)

x

Strain

Stress





Ea

ൌ െ



2





3



GIC

σ = EaԖ

Ic

a

G

E2 Ԗ3
NASA/TM—2010-216824 12 

Figure 10.—Load versus displacement curves for E719/IM7 DCB 

specimen tested by Song and Waas (Ref. 15), along with their 

Finite Element Model and the present joint element model. 

Figure 11.—The solid and hollow circles in the a) load/displacement plot for 

the E719/IM7 DCB specimens correspond to the solid and dotted curves 

in the b) stress and c) strain plots. 

-0.4

0

0.4

0.8

1.2

1.6

0 5 10 15 20 25 30 35 40

Strain, %

Distance, x, along adhesive, mm 

-10

0

10

20

0 5 10 15 20 25 30 35 40

Stress , MPa

Distance, x, along adhesive, mm 

0

50

100

150

200

250

0123456

Load , N

a) Vertical end displacement, mm 

b)

c)

1

2 3

4

1

2 3 4

1

2

3 4
NASA/TM—2010-216824 13 

Figure 12.—Load versus displacement curves for E7T1/G40 DCB 

specimen tested by Song and Waas (Ref. 15), along with their 

Finite Element Model and the present joint element model. 

The behavior of the E7T1/G40 DCB specimens, shown in Figure 12, was not predicted as accurately 

by the joint element model. The joint element model predicted nearly the same elastic slope as the 

Song/Waas FEA model, except with a more conservative prediction of the maximum load reached. 

Song/Waas used an interatomic force law to get the constitutive relations, while the present model used 

the adhesive modulus and critical strain energy release rate. It appears that for this case, the two methods 

resulted in very similar predicted responses. This comparison shows that the joint element can compete 

with discrete cohesive zone models in predicting the behavior of DCB joints. 

4.0 Comparison with Experimental Results 

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 constitutive 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 the NASA Glenn Research Center, and digital image correlation (DIC) techniques (Fig. 13(a)) 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 13(b) 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 2. The functions were 

chosen because they result in the same curve in compression and tension and seem to fit the data 

adequately. 

0

20

40

60

80

100

120

140

160

180

01234

Load, N

Vertical End Displacement, mm 

Joint Element

Song/Waas FEA

Song/Waas Experiment
NASA/TM—2010-216824 14 

Figure 13.—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. 

Figure 14.—A thin adhesive layer with an applied 

stress is in a state of triaxial stress. 

TABLE 2.—MATERIAL PROPERTIES AND GEOMETRIC PARAMETERS OF DCB SPECIMENS 

Specimen Adhesive Adherend Geometric Parameters 

Ea

(GPa)

σu

(MPa)

Ԗfail υ E 

(GPa)

l 

(mm)

a 

(mm)

b 

(mm)

t 

(mm)

η 

(mm)

EA 9394 4.2 49.6 0.016 0.4 69 152.4 63.5 25.4 12.7 0.6 

EA 9309.3NA 2.7 41.3 0.068 0.42 69 152.4 63.5 25.4 12.7 0.55 

The tensile loading of a thin adhesive layer with relatively large in-plane dimensions differs greatly to 

that of a solid cylindrical specimen because the adhesive layer is extremely thin in one direction, and 

constrained from lateral displacement by the top and bottom adherends. Because of these conditions, the 

adhesive is effectively a body in plane strain in the two directions perpendicular to the adhesive thickness 

(Fig. 14). The adhesive is constrained from contracting (Poisson’s effect) in the x and y-directions while 

being loaded in the z-direction, which induces a stress in all three directions, commonly called a state of 

triaxial stress (Ref. 16). To find the stress/strain relation for a material under triaxial stress, consider first 

an isotropic, linearly elastic material. The normal stress in the z-direction is: 

ߪ௭ ൌ ா

ଵାఔ ቂ߳௭ ൅ ఔଵିଶఔ ሺ߳௭ ൅ ߳௫ ൅ ߳௬ሻቃ. (31) 

z

y

x

σz

σz
NASA/TM—2010-216824 15 

The adhesive can be assumed to be in a state of plane stress in the xy-plane, and the strains ߳௫ and ߳௬ can 

be set to zero. Then, the normal stress in the z-direction reduces to: 

ߪ௭ ൌ ଵିఔ

ሺଵାఔሻሺଵିଶఔሻ ܧ߳௭. (32) 

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 non-linear stress/strain relation. Therefore, the stress/strain relation was redefined as: 

ൌ ߪ

ଵିఔ

ሺଵାఔሻሺଵିଶఔሻ ݃ሺ߳ሻ , (33) 

which effectively increases the adhesive modulus. 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 wouldn’t be the point of failure. 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. 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 15. 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 15. The force/displacement curves for all six 

specimens are shown in Figure 16. The high strain-to-failure of EA 9309.3NA caused these specimens to 

hold over two times the amount of load as the EA 9394 specimens. The EA 9309.3NA specimens held so 

much load that the aluminum adherends deformed plastically before the adhesive failed. This can be seen 

in the plot as a gradual rounding prior to a drop in load. 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 lengthened 

by 7.5 percent 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. While not exactly matching the 

Figure 15.—Two adherends of an EA 9394 DCB specimen after complete failure. The fact that 

adhesive is found on both adherends shows that failure occurred within the adhesive layer 

as desired. 

Adherends Adhesive

Air bubbles at crack front

5 mm
NASA/TM—2010-216824 16 

Figure 16.—Load versus displacement curves for DCB specimens 

with aluminum adherends and EA 9394 and EA 9309.3NA 

adhesive, along with the joint finite element model prediction. 

experimental data, the model was rather accurate considering the aforementioned adherend plasticity and 

air bubbles in the end of the adhesive. 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. 

5.0 Conclusion 

This study showed how a single finite element could be used to capture the behavior of an entire joint 

by using an analytical method to solve for the appropriate shape functions rather than prescribe the shape 

functions using a simple polynomial. Additionally, progressive failure was included by defining the stress 

as a non-linear function of the strain and enforcing a strain-based failure criterion. Failure of the adhesive 

was approximated by shortening the joint element by the length of the failed adhesive and lengthening the 

adjoining beam elements. The accuracy of the joint element model was assessed through comparison with 

finite element and experimental DCB test results published by Song and Waas (Ref. 15), and 

experimental DCB tests conducted by the authors. Results showed that the joint element model could 

capture the force/displacement behavior of the joint quite well, with a significant computational efficiency 

advantage compared to a traditional FEA. The assessment also demonstrated how the stress/strain relation 

of the adhesive could be obtained using fracture mechanics or constitutive relations determined by bulk 

tensile tests. When using tensile tests, a factor was used to amplify the stress/strain relation to account for 

additional deformation resistance due to a triaxial state of stress. With proper implementation of this 

method, a single finite element can be used in a broader, global-scale model to capture the behavior and 

failure of joints. 

References 

1. Hart-Smith, L., “Adhesive Bonding of Composite Structures- Progress to Date and Some Remaining 

Challenges,” Journal of Composites Technology and Research, Vol. 24, No. 3, 2002, pp. 133–151. 

2. Frostig, Y., Thomsen, O., Mortensen, F., “Analysis of Adhesive-Bonded Joints, Square-End, and 

Spew-Fillet- High-Order Theory Approach,” Journal of Engineering Mechanics, Vol. 125, No. 11, 

1999, pp. 1298–1307. 

0

300

600

900

1200

1500

1800

02468

Load, N

Vertical end displacement, mm 

EA 9309.3NA

EA 9394

Joint Element Model
NASA/TM—2010-216824 17 

3. Gustafson, P.A. and Waas, T.A. “A Macroscopic Finite Element for a Symmetric Double-lap Joint 

Subjected to Mechanical and Thermal Loading,” Proceedings of the AIAA/ASME/ASCE/AHS/ASC 

48th SDM Conference, April 23–26, 2007, Honolulu, HI, AIAA–2007–2308. 

4. Volkersen, O., “Die Nietkraftverteilung in Zugbeanspruchten mit Konstanten Laschenquerschritten,” 

Luftfahrtforschung, Vol. 15, 1938, pp. 41–47. 

5. Goland, M. and Reissner, E. “The Stresses in Cemented joints,” Journal of Applied Mechanics,

Vol. 11, 1944, pp. A17–A27. 

6. Hart-Smith, L. J. “Adhesive bonded single lap joints,” NASA CR–112236, Douglas Aircraft, 

McDonnell Douglas Corporation, Long Beach, CA, 1973. 

7. Mortensen, F. and Thomsen, O.T., “Analysis of Adhesive Bonded Joints: a Unified Approach,” 

Composites Science and Technology, Vol. 62, 2002, pp. 1011–1031. 

8. Delale, F., Erdogan, F., and Aydinoglu, M.N., “Stresses in Adhesively Bonded Joints: A ClosedForm Solution,” Journal of Composite Materials, Vol. 15, 1981, pp. 249–271. 

9. Zhang, J, Bednarcyk, B.A., Collier, C, Yarrington, P., Bansal, Y., and Pindera, M.J., “3D Stress 

Analysis of Adhesively Bonded Composite Joints,” Proceedings of the AIAA/ASME/ASCE/AHS/ 

ASC 46th SDM Conference, April 18–21, 2005, Austin, TX, AIAA–2005–2021. 

10. Stapleton, S.E. and Waas, A.M., “Macroscopic Finite Element for a Single Lap Joint,” Proceedings of 

the AIAA/ASME/ASCE/AHS/ASC 50th SDM Conference, May 4–7, 2009, Palm Springs, CA, 

AIAA–2009–2449. 

11. ABAQUS, Inc, ABAQUS user Manuel v 6.7-1, Electronic Version, 2006. 

12. Tsai, M.Y. and Morton, J., “An evaluation of Analytical and Numerical Solutions to the Single-Lap 

Joint,” Int. J. Solids Structures, Vol. 31, No. 18, 1994, pp. 2537–2563. 

13. Chen, C., Linear System Theory and Design, 3rd ed., Oxford University Press, New York, NY, 1999. 

14. Molder, C. and Van Loan, C., “Nineteen Dubious Ways to Compute the Exponential of a Matrix, 

Twenty-Five Years Later,” SIAM Review, Vol. 45, No. 1, 2003, pp. 3–49. 

15. Song, S.J., and Waas, A.M., “Mode I Failure of Laminated Polymeric Composites,” Engineering 

Fracture Mechanics, Vol. 49, No. 1, 1994, pp. 17–27. 

16. Sechler, E., Elasticity in Engineering, Wiley: Chapman & Hall, New York, NY, 1952. 
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-11-2010 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

Modeling Progressive Failure of Bonded Joints Using a Single Joint Finite Element 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Stapleton, Scott, E.; Waas, Anthony, M.; Bednarcyk, Brett, A. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 727950.04.03.22 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-17453 

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/TM-2010-216824 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Category: 39 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Enhanced finite elements are elements with an embedded analytical solution which can capture detailed local fields, enabling more efficient, 

mesh-independent finite element analysis. In the present study, an enhanced finite element is applied to generate a general framework 

capable of modeling an array of joint types. The joint field equations are derived using the principle of minimum potential energy, and the 

resulting solutions for the displacement fields are used to generate shape functions and a stiffness matrix for a single joint finite element. 

This single finite element thus captures the detailed stress and strain fields within the bonded joint, but it can function within a broader 

structural finite element model. The costs associated with a fine mesh of the joint can thus be avoided while still obtaining a detailed 

solution for the joint. Additionally, the capability to model non-linear adhesive constitutive behavior has been included within the method, 

and progressive failure of the adhesive can be modeled by using a strain-based failure criteria and re-sizing the joint as the adhesive fails. 

Results of the model compare favorably with experimental and finite element results. 

15. SUBJECT TERMS

Composites; Sandwich panels; Preliminary design; Local buckling; Structural sizing; Laminates 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

23 

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