Scott E. Stapleton and Thomas Gries

RWTH Aachen University, Aachen, Deutschland

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Evan J. Pineda

Glenn Research Center, Cleveland, Ohio

Adaptive Shape Functions and Internal Mesh

Adaptation for Modelling Progressive Failure

in Adhesively Bonded Joints

NASA/TM—2014-218333

September 2014
NASA STI Program . . . in Prole

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientic 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 Ofcer. 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 signicant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signicant 

scientic 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. Scientic 

and technical 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. Scientic and 

technical ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientic and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientic, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientic 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 to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
Scott E. Stapleton and Thomas Gries

RWTH Aachen University, Aachen, Deutschland

Anthony M. Waas

University of Michigan, Ann Arbor, Michigan

Evan J. Pineda

Glenn Research Center, Cleveland, Ohio

Adaptive Shape Functions and Internal Mesh

Adaptation for Modelling Progressive Failure

in Adhesively Bonded Joints

NASA/TM—2014-218333

September 2014

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135
Acknowledgments

Portions of this work were 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 nancial support was provided by 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://www.sti.nasa.gov

Trade names and trademarks are used in this report for identication 

only. Their usage does not constitute an ofcial 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—2014-218333 1 

Adaptive Shape Functions and Internal Mesh Adaptation for 

Modelling Progressive Failure in Adhesively Bonded Joints 

Scott E. Stapleton and Thomas Gries 

RWTH Aachen University 

Aachen, Deutschland 52074 

Anthony M. Waas 

University of Michigan 

Ann Arbor, Michigan 48105 

Evan J. Pineda 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135 

Abstract 

Enhanced finite elements are elements with an embedded analytical solution that can capture detailed 

local fields, enabling more efficient, mesh independent finite element analysis. The shape functions are 

determined based on the analytical model rather than prescribed. This method was applied to adhesively 

bonded joints to model joint behavior with one element through the thickness. This study demonstrates two 

methods of maintaining the fidelity of such elements during adhesive non-linearity and cracking without 

increasing the mesh needed for an accurate solution. The first method uses adaptive shape functions, where 

the shape functions are recalculated at each load step based on the softening of the adhesive. The second 

method is internal mesh adaption, where cracking of the adhesive within an element is captured by further 

discretizing the element internally to represent the partially cracked geometry. By keeping mesh adaptations 

within an element, a finer mesh can be used during the analysis without affecting the global finite element 

model mesh. Examples are shown which highlight when each method is most effective in reducing the 

number of elements needed to capture adhesive nonlinearity and cracking. These methods are validated 

against analogous finite element models utilizing cohesive zone elements. 

Introduction 

With the increasing demand for fiber reinforced composites in lightweight aerospace structures, 

adhesively bonded joints are becoming more critical than ever. Bolted and riveted joints have proven to be 

poorly suited for composite materials (Ref. 1). Unlike traditional metals, brittle fibers often do not yield 

significantly to spread concentrated loads introduced by mechanical fasteners. Furthermore, bolts and rivets 

require holes in the material to be joined, which interrupts continuous fibers and introduces additional stress 

concentrations. Adhesive bonding is suitable for composite materials because it is less invasive, introduces 

load more gradually, and can often be much more cost effective. The adhesive market has indeed grown 

along with the advanced composite market, and the structural adhesive market in Europe has been 

forecasted to reach 67,000 tons by 2015; a growth of over 13 percent since 2008 (Ref. 2). 

However, adhesively bonded joints can be problematic to model. The models often do not scale because 

of fixed thickness requirements of the adhesive layers, making individual design for each joint necessary. 

Geometric discontinuities in adherends cause stress singularities in many models, thus non-traditional 

failure criteria or evaluation methods are often required. A lack of confidence in material models, failure 

criteria, and engineering experience results in gross over-design of joints along with safety requirements 

which sometimes require secondary mechanical fasteners, jokingly referred to as “chicken bolts” by many 

engineers. 
NASA/TM—2014-218333 2 

Adhesively bonded joints are typically analyzed using analytical models (closed-form) or numerical 

models (finite elements). Historically, analytical models (Refs. 3 to 12) were relied upon exclusively while 

computer capability was relatively small. Analytical models are fast, and have been used to conduct 

numerous parametric studies to further the understanding of the design of joints. However, assumptions are 

often made which allow closed-form solutions but limit joint geometry and materials which can be 

accurately analyzed. Furthermore, these analytical models do not couple well with larger component and 

vehicle models, limiting their usefulness. 

Finite element (FE) models, on the other hand, are general enough to allow a wide range of geometries 

and configurations, and can even couple joint analysis with larger models. The improving speed of 

computers makes this method more viable for joint analysis. However, there are some downsides to 

modeling joints with FE models. The reentrant corners often cause a geometric singularity, and numerous 

work has been conducted to create failure theories which account for this (Refs. 13 to 19). Failure methods 

investigated include stress-based, strain based, plastic energy density, and stress measured at a characteristic 

distance from the singularity. Furthermore, the extremely thin adhesive layer limits the size of elements 

which can be used to explicitly model the adhesive. This means that there are truly no coarse models, and 

coupling with vehicle-scale models can be problematic (Refs. 15 and 16). 

One relatively new technique for modeling progressive failure of adhesively bonded joints is 

progressive damage modeling incorporating fracture mechanics concepts. Interface elements using different 

methods such as discrete cohesive zone method (DCZM) or continuous cohesive zone method (CCZMare) 

used to resolve the stress singularity at material interfaces and reentrant geometrical corners, and allow the 

faces of the adherends to separate by treating the adhesive as a network of non-linear springs obeying a 

traction versus separation law (Refs. 20 to 26). These methods come in many different varieties, but most 

often involve stress-based initiation criterion and energy-based failure definitions. These methods have been 

shown to be extremely powerful for joints, but still have some drawbacks. Characterizing these interface 

elements requires a large amount of characterization tests, and appropriate handling of mode-mixity is also a 

subject within the cohesive zone modeling that has yet to be satisfactorily concluded (Ref. 25). Additionally, 

the cohesive laws utilized require an initial, numerical fictitious stiffness to prevent separation of the plies 

before delamination initiation. Furthermore, there is a maximum element length, based on the process zone 

size, required to obtain accurate results. Cohesive zone element utilizing shape functions that are enriched 

by an analytical solution (similar to the methodology presented in this work) have been used to alleviate 

mesh dependence and size requirements (Refs. 25 and 26). Most cohesive elements are formulated 

assuming a zero-thickness interface, and thus may not be adequate to model adhesive joints, especially if the 

adhesive layer is thick. 

Another fracture based technique for modeling crack propagation that can be applied to joint analysis is 

the virtual crack closure technique (VCCT) [KRUEGER, 2004]. With VCCT fracture toughness based 

criteria are used to determine if it is energetically favorable for a crack to propagate. Propagation is 

restricted to element boundaries and typically must be known a priori. Thus, modeling joints with a finite 

thickness adhesive may prove challenging with VCCT. 

All of the aforementioned methods are highly developed and have been shown to give a reasonable 

strength prediction for joints, but they are detailed models which require extremely fine meshes. Thin 

adhesive bonds, most often thinner than 1 mm, restrict the size of elements needed for the adhesive. The 

transition from the fine adhesive mesh to the coarser adherend mesh causes additional preprocessing work 

for the analyst. Therefore, joint design and analysis is typically completed after the global vehicle sizing on 

dense meshed sub-models, when design changes are expensive or impractical. 

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 (Ref. 27), an area that is getting increased 

attention in the aerospace industry with the aim of lowering design cycle and certification costs. 
NASA/TM—2014-218333 3 

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 (Refs. 28 and 29). More recently, Waas and Gustafson (Ref. 30) have created an element to 

capture the behavior of a double overlap joint subjected to mechanical and thermal loads. 

A general bonded joint finite element has been created (Refs. 30 to 33) wherein an entire bonded joint 

can be modeled with a single element. This joint element considers the adherends to behave like beams (or 

wide panels), and the adhesive to be made up of a bed of shear and normal nonlinear springs. The governing 

equations of this structural model are found and solved to produce enhanced shape functions for the joint 

element. Furthermore, the element has been generalized to allow multiple adherend/adhesive layers and ply 

drops/thickness tapers, providing the capability to model a variety of joint types with very few elements. 

This model was implemented in the software Joint Element Designer, which was written in C# and first 

conceived in a joint effort between the University of Michigan and NASA. However, this method loses its 

advantage when modeling highly nonlinear adhesives and trying to capture progressive failure. An increase 

in elements is required for an accurate solution, which goes against the philosophy of enhanced elements. 

This paper presents two methods which allow the bonded joint finite element to capture adhesive 

non-linearities and cracking without increasing the mesh needed for an accurate solution. The first method is 

the use of adaptive shape functions, where the shape functions are recalculated at each load step based on 

the softening of the adhesive. The second method is internal mesh adaption, where cracking of the adhesive 

within an element is represented by discretizing a cracked element into multiple elements in order to 

accurately represent the local, cracked geometry. Both of these methods were implemented in the Joint 

Element Designer software. Examples are shown which highlight the savings in elements, computational 

time, and integration points needed when using these methods and when the methods are particularly 

beneficial. The performance of the various joint element methodologies are compared to analogous models 

using CZM elements. 

Finally, while these methods are shown for an adhesively bonded joint but both methods have broader 

application. The adaptive shape functions could be used for any element where the material and geometric 

properties used to obtain the shape functions are changing. Updating the shape functions within an analysis 

would improve the ability of the shape functions to represent the deformation of the changing material/ 

geometry. Finally, adaptive mesh, which is not a new technique (Refs. 25, 26, 34, and 35), can be very 

effective in capturing discontinuities within an element without changing the elements relationship with the 

global model.

Formulation 

The formulation presented here is merely a summary of the formulation of the joint element, and the 

details presented are considered by the authors to be required for the understanding of the current study. For 

a more detailed formulation, the reader is encouraged to consult prior work (Ref. 36). 

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 (Figure 1(a)). The plates are assumed to behave 

as “wide” layered composite Euler Bernoulli beams (hence the cylindrical bending assumption). The axial 

displacement of adherend i, ui, varies linearly through the thickness and the transverse and rotational 

displacements, wi and wi,x, are assumed to be constant through the cross section. Additionally, it is assumed 

that the only significant stresses are the axial stresses induced by bending and axial deformation, so other 

stress components in the adherend can be neglected. The adhesive joining the plates is modeled as a Winkler 

foundation, where only peel and shear stress are included and the stresses are assumed to be constant though 

the thickness of the adhesive layer. With the displacement behavior known in the transverse, zi, direction, 

the problem can be formulated in terms of one dimension, x. 
NASA/TM—2014-218333 4 

Figure 1.—Adhesively bonded joint element (a) geometric parameters and (b) finite element discretization. 

Shape Functions for Linear Elastic Case 

The shape functions used for the joint element were derived by analytically solving the governing 

equations of the linear elastic joint model (Refs. 33 and 36). Using the principle of stationarity of 

potential energy and the aforementioned approximations, 2N fully coupled governing equilibrium 

differential equations are obtained. 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 (1) 

where u is a vector containing the adherend centerline vectors of all of the N adherends: 

1

T TT T

i N     uu u u (2) 

and the centerline vector of adherend i is defined as 

, ,, , () () () () () () T

i i i x i i x i xx i xxx   ux ux wx wx wx wx   u (3) 

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. 

One downside to this method worth mentioning is its dependence on solving the system of ODE’s to 

get the shape functions. This makes the extension of this method to a plate or shell-type element difficult 

because such an element would require the solution of a PDE, which can be significantly more difficult 

and time-consuming. 

The system of ordinary differential equations in Equation (1) is solved using the matrix exponential 

and rewritten in terms of the nodal degrees of freedom, q (Figure 1(b)), and shape functions, N, in the 

form: 

u Nq (4) 
NASA/TM—2014-218333 5 

The main thing to consider is that the shape functions were not prescribed, as is normally the case in finite 

elements, but derived by solving the governing equations for the stacked beam model. Therefore, the shape 

functions are a function of the geometric and material parameters of the adherend and adhesive layers: 

N N( , , ) ,( , , , ) Etb E G b i i i i N ai ai ai ai ai N  1.. 1.. 1 	  (5) 

Material and Geometric 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: 

 ( ) (6) 

where the stress is some general function of the strain in the adhesive and adherends. 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. 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 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. 

The elemental stiffness matrix, k , is found by integrating the: 

( ) T

V dV  k B DqB (7) 

where B is the matrix linking the strains and nodal degrees of freedom (which contains the shape 

functions for the initial linearly elastic state: N), and D q( ) is the nonlinear tangent stiffness matrix. 

Additionally, a co-rotational formulation was used to account for large rotations in the system, as 

described in References 37 to 39 and adapted to the joint element in Reference 36. 

Adaptive Shape Functions 

When the adhesive in a joint element has a nonlinear stress-strain relationship, the shape functions 

obtained for a linearly elastic adhesive may no longer be well-suited. The shape functions were derived 

based on a solution of the governing equations of the problem, so the stiffness of the materials is used to 

acquire the shape functions. As the adhesive softens, the tangent modulus changes and the shape 

functions found before loading no longer represent the solution of the governing equations based on the 

softened adhesive tangent modulus. This is illustrated in Figure 2, where a joint with a highly ductile 

adhesive is loaded and the shape functions (Figure 2(d)) are required to change as the adhesive softens. 

For highly nonlinear materials, many elements may be required to accurately capture the joint behavior. 
NASA/TM—2014-218333 6 

Figure 2.—Demonstration of shape function progression during a loading scenario: (a) single lap joint under 

tension, (b) load versus displacement highlighting three points, and (c) tangent modulus of the adhesive and 

(d) selected shape functions at each of the three points. 

To address this issue, 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 where the stiffness of the 

joint changes along the x-direction. The function for the adhesive modulus, Eai(x), is based on the tangent 

stiffness of the adhesive in the prior load increment. The peel stress in adhesive i can be written in terms 

of the peel strain as 

( ) ii i  
 aaa E x (8) 

and similarly for the shear stress: 

( ) ii i   aaa G x (9) 
NASA/TM—2014-218333 7 

Figure 3.—To solve the governing equations for (a) an element with a softened adhesive, the joint element was 

(b) discretized into regions (in this case 6) within which the tangent modulus was considered constant. 

This causes the system of governing equations to now have non-constant coefficients (compare with 

Eq. (1)) 

, ( ) x u Au x (10) 

where the coefficient matrix, A( ) x , is now a function of x. Since the coefficient matrix is non-constant, a 

semi-numerical method of solution was adopted. The domain was split into segments in which the 

coefficient matrix, A( ) x , was considered constant and solved using the matrix exponential for that 

section (Figure 3). The sections were than linked together with boundary conditions, and an approximate 

solution of A( ) x was found (Ref. 32). 

The number of segments was determined by the number of integration points, so that the shape 

functions were refined with an increasing number of integration points. 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 very few elements, simplifying analysis steps like mesh generation and post-processing. 

The shape functions were re-calculated at each loading step in the analysis. This adds to the 

computational time, but can be extremely advantageous in reducing the required mesh for highly 

non-linear adhesives. 

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 traditional finite elements, more 

elements are required to accurately find the solution where the adhesive has failed. In the case of failed 

adhesive, a great number of elements may be needed, as will be illustrated later. 

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 as shown in Figure 4 and outlined in Figure 5. 

0

1

0 1 x /l

a) b)

0

1

0 1 x /l

tan

lin

E

E

A( )x

A1 A2 A3 A4 A5 A6 A7

A8

tan

lin

E

E
NASA/TM—2014-218333 8 

Figure 4.—Diagram showing (a) an uncracked joint element, (b) a partially cracked element, and 

(c) a fully cracked joint element, and the addition of internal sub-elements. 

Figure 5.—Flow chart showing how cracked element sub-assembly is incorporated into joint element solution 

procedure. 

As with standard finite element analysis, a Newton-Raphson type solver is used to solve the global 

system of equations. At load increment j, the global degrees of freedom (dofs), Q j , are guessed, and this 

is broken down into elemental dofs, q j . After applying the prescribed displacements to the joint element, 

the adhesive is checked for failure or cracking. 

If failure in the adhesive is detected, the element is replaced by a sub-assembly with three elements as 

shown in Figure 4(b). The length of the crack, lcrack, determines the lengths of the sub-assembly elements. 

The displacement of the inner nodes is unknown, therefore the sub-assembly becomes a nonlinear model 

within another nonlinear model and must be solved with its own Newton-Raphson type solution 

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 Kj Fj

no

yes

k j , Int

j f

find 

Qj+1 from 

Kj and 

Rj: j = j + 1

for each element

Apply Loads, BCs

Find . Rj

Element cracked?

no

calculate

Newton-Raphson

Guyan

Reduction

Qj

qj

Sufficient accuracy?

Assemble and . Int K j Fj

no

yes

k j , Int

j fcalculate

find Qj+1 from Kj and Rj

Element

for each element

Apply Loads, BCs

Find . R j
NASA/TM—2014-218333 9 

procedure. A sub-assembly Newton-Raphson type solver is employed, the elemental dof are used as outer 

node boundary conditions, and the stiffness, KSub , and internal force vector, Int FSub , for the sub-assembly 

are 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 (Refs. 40 to 42). 

Once the internal degrees of freedom are removed, the stiffness matrix and force vector, k j and Int

j f , can 

be considered to be that of the equivalent joint element, and are passed on to the global assembly. 

These vectors/matrices for all of the elements in the assembly are assembled, boundary conditions 

and loads are applied, and the residual, R j , (error of the initial nodal displacement guess) is calculated. If 

the residual is not 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. 

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 

sub-assembly elements, and the global system is re-solved. This is done until no new adhesive failure 

occurs and 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  
 (11) 

where prev

crack l is the previous crack length (prior to the global Newton-Raphson procedure) and cur

crack 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. 

Method 

The aim of this study was to introduce the concept of adaptive shape functions and an internal 

adaptive mesh for the joint element and show how these methods can be used to reduce the amount of 

elements or integration points needed to reach a converged solution for the failure of adhesively bonded 

joints. To demonstrate the effect that these methods have, two sets of examples will be introduced. First, 

the bonded joint element model will be compared with and without adaptive mesh and adaptive shape 

functions to show when these methods are most effective at reducing computational time and how 

effective they can be. Second, the bonded joint element model is validated against a commercial CZM 

model and the computational times will be compared. 
NASA/TM—2014-218333 10 

Joint Element Model Comparison 

First, the bonded joint element model was used with and without adaptive shape functions (ASF) and 

adaptive mesh (AM) to see the effect of these methods. Four variations of the joint element model were 

compared (Table 1). The first one, standard joint elements, has no ASF or an AM. The shape functions 

are calculated based on the initial, undamaged state of the element. The adhesive is considered to be 

failed when the criterion is reached, and load carrying capability of the failed section of adhesive is set to 

zero in both shear and peel. The ASF model is the same as the standard model except with the shape 

functions being recalculated as the adhesive softens, and the AM model is the same as the standard model 

except that the failure of adhesive results in an introduction of internal nodes within the element and the 

position of these nodes changes as the crack in the adhesive grows. Finally, the last model, adaptive shape 

functions and mesh (ASF+M) combines both techniques. 

Two different joints are analyzed, a double cantilever beam (DCB) and a single lap joint 

(see Figure 6(a) and (b), respectively). The DCB joint was selected because the crack typically grows 

progressively through the joint and the adhesive is loaded in mode I only. The single lap joint analyzed in 

this study results in catastrophic failures (no progressive failure) and mixed mode loading. The boundary 

conditions, loading and geometry of the two joints are shown in Figure 6. 

For each joint configuration, two different adhesives are analyzed to show how effective adaptive 

shape functions and adaptive internal mesh methods are for highly plastic and brittle adhesives. The 

adhesives had an initial Young’s modulus of 4 GPa, a maximum stress of 50 MPa, an initial shear 

modulus of 1.4 GPa, and a maximum shear stress of 29 MPa. The first adhesive displayed material 

nonlinearity only, and the stress-strain relationship was defined by a tanh function with the initial slope 

and maximum stress dictated by the aforementioned properties (see Figure 7(a)). There was no failure of 

the adhesive, it was allowed to strain indefinitely. This case is interesting because it reflects the global 

yielding failure theory introduced by Crocombe (Ref. 13), which states that the upper bound of the 

strength of a joint can be found by determining the point at which the whole adhesive layer has a tangent 

modulus of zero. The second adhesive used in the study had a linear stress-strain relationship up until 

failure (Figure 7(b)). As before, the slope of the stress-strain relationship was defined previously, along 

with the stress at failure. 

TABLE 1.—FOUR JOINT ELEMENT MODEL 

VARIATIONS COMPARED 

Joint element 

model variation 

Adaptive shape 

functions 

Adaptive 

mesh

Standard - - 

ASF X - 

AM - X 

ASF + M X X 

Figure 6.—Joint configurations used in the joint model comparison: (a) double 

cantilever beam (DCB) joint and (b) single lap joint (SLJ), where both joints 

are 20 mm wide (not to scale). 

20 mm 50 mm

F, 

0.2 mm

5 mm

20 mm 50 mm F, 

5 mm

0.2 mm

a)

b)
NASA/TM—2014-218333 11 

Figure 7.—Adhesive stress/strain relations used to demonstrate the convergence improvements using adaptive 

shape functions and adaptive mesh for adhesives displaying (a) material nonlinearity and (b) brittle material 

behavior. 

Two values were compared between the variations of the model: integration points and computational 

time to convergence. Convergence is reached when the error of the analysis is less than 1 percent. The 

error is defined in one of two ways. For joint configurations which display progressive crack growth 

(i.e., DCB joint configuration with both adhesives and single lap joint with the adhesive with nonlinear 

material properties but no failure) the error, e1, is defined as 

1

1

1 () () 100

( )

N

i i

i i

F F

e

N F






    (12) 

where N is the number of load increments, ( ) F i is the end force at load increment i, and ( ) F i is the 

end force for the fully converged solution at load increment i. This is a measure of the average percent 

deviation of the load-displacement plot in the load-direction. When the joint fails catastrophically, or with 

no progressive crack growth, the error is taken as simply the percent deviation of the maximum load, 

or 

max max

2

max

100

F F

e

F






 (13) 

where e2 is the error, Fmax is the maximum load (i.e., strength) and Fmax

 is the strength of the fullyconverged solution. The analyses were all run on the same personal computer with 8 GB of RAM. While 

the actual computational time should be considered more of a qualitative property, the computational 

times of the different models relative to one another should be considered significant.

Finally, to identify circumstances under which using adaptive shape functions becomes highly 

advantageous, a single lap joint with an adhesive containing a trapezoidal stress-strain relation and 

different levels of strain energy before failure as shown in Figure 8 was analyzed. The initial slope and 

maximum stress of the stress-strain relation is the same as that defined for the previously studied 

adhesives and the final slope was a negative value of the initial slope. The total area under the curve, 

(WIc for mode I and WIIc for mode II) was set as an integer multiple of the area the curve for the equivalent 

triangular stress-strain relation (WaI for mode I and WaII for mode II). Mode I and mode II strain energies 

were both increased, so that Wc/Wa represents the value of WIc/WaI and WIIc/WaII simultaneously. AM and 

0

25

50

0 0.02 0.04 0.06

Stress (MPa)

Strain

0

25

50

0 0.02 0.04 0.06

Stress (MPa)

a) b) Strain

Peel

Shear

Peel

Shear
NASA/TM—2014-218333 12 

Figure 8.—Trapezoidal stress/strain relations for (a) normal and (b) shear used to show the benefits of using 

adaptive shape functions. The length of the horizontal region of the trapezoid was extended or shortened to 

get different Wc/Wa ratios. The current energy (c) defined for mode I and II. 

ASF+M models were used, and computational time (t) and integration points (n) needed to reach 

convergence was compared in order to see how adaptive shape functions effects the computational time 

for adhesives with different amounts of critical strain energy. 

The adhesive was defined as failed when the following relation was satisfied: 

1 I II

Ic IIc

W W

W W

  (14) 

where WI and WII are defined as the current strain energy, as shown in Figure 8. It should be noted that 

this is synonymous to mixed-mode laws used for cohesive zone element models, except with fracture 

energy rather than strain energy. Since the adhesive of the joint element was assigned a thickness, we 

speak of everything in terms of stress/strain rather than traction/separation. However, since the 

displacements in the adhesive are assumed to be linear in the thickness direction of the adhesive, the 

fracture energy and separation can be obtained by simply scaling the strain energy and strain by the 

adhesive thickness. If the model is extended to include displacements of a higher order through the 

thickness, the traction/separation and stress/strain relations can no longer be linked. Therefore, we use the 

stress/strain relation for generality. 

Comparison With Cohesive Zone Model 

The joint element models were also compared to a cohesive zone model in Abaqus/Standard v6.11-1 

finite element software (Ref. 7) using pre-existing modules within the software. The implicit, static finite 

element solver was used and a finite element mesh convergence study was performed to ensure numerical 

convergence was achieved. The adherends in the Abaqus models were discretized with 2D, linear, 

quadratic, plane strain, incompatible modes, CP4EI elements. These elements were chosen to better 

represent the linear bending action in the adherends. The adherends were assigned linear elastic material 

properties. 

Cohesive elements (COH2D4) were used to represent the adhesive layers in the Abaqus models. 

These 2D elements were given a finite thickness according to the thickness of the adhesive and the 

constitutive response was governed using a traction-separation description. Element tractions are related 

to the nodal separation through a traction-separation law. Typically, the area under the traction-separation 

curve is equal to the fracture toughness of the material. For these analyses, a bi-linear traction-separation 

law was chosen. A simple, maximum stress criterion was used to mark the initiation of softening within 

the element, and the power law, mixed mode implementation within Abaqus was used to couple the 

opening and shear fracture modes. Both the Abaqus and joint models used the same adhesive strengths 

used in the earlier examples. The mode I and mode II fracture toughnesses employed in the Abaqus model 

were chosen such that they were equal to the critical strain energy densities, used in the joint models, 

scaled by the thickness of the adhesive (GIC = 625 N/m, GIIC = 601 N/mm). 

0

25

50

0 0.02 0.04 0.06

WaI WaII

0

25

50

0 0.02 0.04 0.06 0.08

0

25

50

0 0.02 0.04 0.06

WIc WIIc

a) b)

Wi

Wic

c) a,  a

a,a

a

a (MPa)

a

a (MPa)
NASA/TM—2014-218333 13 

The Abaqus built-in cohesive elements were considerably less flexible than the joint element, so the 

joint element failure criterion and constitutive relations were altered to match the cohesive elements. For 

the case of single mode loading (DCB: mode I only), the constitutive relation shown in Figure 9(a) was 

used for both models (except the cohesive zone model was scaled by the adhesive thickness to get a 

traction-separation law). The critical strain energies were the same as for the trapezoidal example when 

Wia/Wic = 5. 

Since the joint element model considers the adhesive to be made of uncoupled shear and normal 

springs, the constitutive law had to be adjusted to give similar results to that of the cohesive element 

model for mixed-mode loading (SLJ). The built in cohesive zone model used a quadratic power law 

softening criterion, where softening was initiated when 

2 2

1 a a

Ic IIc

 

 

       

  

 (15) 

To adapt the joint element constitutive relation to this softening criterion, the shear to peel ratio,  , was 

assumed to be constant throughout the loading and defined as 

aMax

aMax





 (16) 

where  aMax and aMax  were the maximum peel and shear stresses in the adhesive found through a 

linearly elastic analysis. Finally, using Equations (15) and (16), the adjusted peel and shear stresses at 

softening,  c and c  , were found using 

222

1 1

Ic Ic IIc



 

          

      and     Ic IIc . (17) 

For the cohesive zone model, the downward slopes for mode I and mode II are adjusted so that both 

modes reach zero traction when Equation (15) (peel and shear stress being replaced with mode I and 

mode II fracture energies) is satisfied, assuming that the ratio of peel to shear occurring at the onset of 

softening remains constant. For the joint element, this translates into cracking of the adhesive when: 

2 2

1 I II

Ic IIc

W W

W W

       

  

 (18) 

Assuming the ratios of peel to shear stress and peel to shear strain remain constant, this yields the 

following adjusted critical strain energies: 

2 2 2 1 1 a

Ic Ic a IIc

E

W W GW

                  and 2

a IIc Ic

a

G W W

E 

  . (19) 

These relation yields the mixed-mode adjusted constitutive relation shown in Figure 9(a). 
NASA/TM—2014-218333 14 

Figure 9.—Constitutive law joint element for comparison with cohesive elements, (a) with single-mode constitutive 

laws and (b) mixed-mode adjusted constitutive law for  = 1.08. 

Results and Validation 

Adaptive shape functions (ASF) and an adaptive internal mesh (AM) were introduced to the joint 

element to reduce the number of elements needed to capture adhesive nonlinearity and crack growth when 

modeling an adhesively bonded joint. The following section illustrates the benefit of these two methods in 

reducing computational time, elements required, and integration points required. The first section 

compares the joint element model with and without these features for different joint configurations and 

adhesive types, while the final section compares the model with a commercial code cohesive zone model. 

Joint Element Model Comparison 

Double Cantilever Beam (DCB) 

The first joint configuration to be studied was the DCB configuration. For this configuration, the 

adhesive only experiences mode I loading and crack growth occurs in a progressive manner. This analysis 

is useful in showing how the different variations of the joint element model behave under progressive 

failure and crack growth. The number of elements, number of integration points, and computational time 

for both adhesives are shown in Table 3. Corresponding load versus displacement plots are shown in 

Figure 10. Since the failure is progressive and maximum load is only one integral part of the results, the 

error measurement e1 is used and convergence occurs when the error falls below 1 percent. 

The first adhesive considered was the ductile adhesive with no failure, but softening only. Since no 

crack forms for this adhesive, adaptive mesh is not needed and only standard elements and elements with 

adaptive shape functions are compared. The adaptive shape functions were extremely effective in 

reducing the number of integration points and computational time required to get a converged solution, 

with a computational time two orders of magnitude smaller. However, even using adaptive shape 

functions, the converged solution could not be captured with merely one element, but three were required. 

The second adhesive considered was a brittle adhesive, with the stress-strain relation for mode I and 

mode II remaining linear until failure. The standard joint elements and adaptive shape function models 

both performed similarly, with similar computational times. However, the adaptive mesh method was 

pivotal in reducing the computational time significantly. Since there was no nonlinearity before failure in 

the adhesive, the adaptive shape functions did not affect the analysis except by adding to the 

computational time since the shape functions had to be recalculated at every load step. 

0

10

20

30

40

50

0 0.05 0.1 0.15 0.2

Stress (MPa)

Strain

0

10

20

30

40

50

0 0.05 0.1 0.15 0.2

Stress (MPa)

a) b) Strain

IIc 

 Ic

 Ic 

 1.08

IIc 

WIc

WIIc WIIc

 WIc


NASA/TM—2014-218333 15 

TABLE 2.—COMPARISON OF NUMBER OF ELEMENTS (nel), NUMBER OF INTEGRATIONS POINTS (n), 

AND COMPUTATIONAL TIME (t) FOR THE CONVERGED SOLUTION (e1  1 PERCENT) FOR THE 

DCB JOINT WITH DIFFERENT MODELS AND ADHESIVES 

Model 

Ductile adhesive Brittle adhesive 

el n n t (sec) el n n t (sec) 

Standard joint elements 40 680 4435 40 680 1988 

Adaptive shape functions (ASF) 3 25 71 1 1025 1498 

Adaptive mesh (AM) ---- ---- ---- 1 65 33 

ASF + M ---- ---- ---- 1 65 112 

Figure 10.—DCB joint configuration with load versus displacement plot for the models with (a) a ductile and 

(b) a brittle adhesive. 

Single Lap Joint (SLJ) 

The single lap joint was the next joint configuration to be considered, which is a more realistic joining 

configuration containing mode I and mode II loading in the adhesive, and catastrophic crack growth 

rather than progressive crack growth, except for the case of the nonlinear adhesive with no failure. 

When catastrophic failure occurred, it was more practical to use the error value e2, since it deals with the 

maximum load only. The number of elements, number of integration points, and computational time for 

both adhesives are shown in Table 3. Corresponding load versus displacement plots are shown in 

Figure 11. 

The first adhesive considered for the single lap joint configuration was the nonlinear adhesive with no 

failure. The maximum load is reached when the adhesive becomes fully softened in shear (the tangent 

shear modulus is nearly zero along the entire length of the adhesive). Similar to the DCB configuration, 

models with adaptive mesh were not considered because no crack growth modeling was required. As can 

be seen, using adaptive shape functions for this case made a huge difference, with two orders of 

magnitude less integration points required and three orders of magnitude faster analysis. Obviously, for 

this case, it would be highly advantageous to use adaptive shape functions. 

The second adhesive considered for the single lap joint was the brittle adhesive, which is linear until 

failure. As with the DCB configuration with the same adhesive, adaptive shape functions are not helpful 

because the adhesive has no nonlinearity, so recalculating the shape functions at every step is 

unnecessary. However, the adaptive mesh model was extremely fast for this configuration. Since the 

adhesive stress-strain relation was linear, the solution could be found in one Newton-Raphson iteration, 

with the internal mesh creation only adding a small amount of additional time. 

0

0.05

0.1

0.15

0.2

0 0.05 0.1 0.15 0.2

F (kN)

(mm)

Stress

Strain

Adhesive Material

F, 

0

0.5

1

1.5

2

0 0.5 1 1.5

F (kN)

a) b) (mm) Stress

Strain

Adhesive Material

F, NASA/TM—2014-218333 16 

TABLE 3.—COMPARISON OF NUMBER OF ELEMENTS (nel), NUMBER OF INTEGRATIONS POINTS (n), 

AND COMPUTATIONAL TIME (t) FOR THE CONVERGED SOLUTION (e1/e2  1 PERCENT) FOR THE 

SLJ JOINT WITH DIFFERENT MODELS AND ADHESIVES 

Model 

Ductile adhesive Brittle adhesive 

el n n t (sec) el n n t (sec) 

Standard joint elements 80 1360 35255 40 680 11521 

Adaptive shape functions (ASF) 2 10 71 40 680 12364 

Adaptive mesh (AM) -- ---- ------- 1 33 2 

ASF + M -- ---- ------- 1 33 25 

Figure 11.—SLJ joint configuration with load versus displacement plot for the models with (a) a ductile and 

(b) a brittle adhesive. 

Trapezoidal Adhesive With Variable W

Considering the foregoing results, one can hypothesize that the additional cost of calculating the 

adaptive shape functions at every load step is beneficial for adhesives with high fracture toughness. To 

further investigate this hypothesis, the adaptive mesh model was compared with and without adaptive 

shape functions for an adhesive with a trapezoidal stress-strain relation and variable amounts critical 

strain energy W, which in this study is synonymous with the fracture toughness normalized by the 

adhesive layer thickness. The strain energy for each mode (WIc and WIIc) is normalized by the strain 

energy of an adhesive with a triangular stress-strain relation with the same critical stress and slopes. The 

computational times (t) and integration points (n) needed for a converged solution (e2  1 percent) for the 

two models are compared in Figure 12. 

In Figure 12 the computational time and integration points for the converged solution of the adaptive 

mesh model is normalized by that of the adaptive shape function + adaptive mesh model. For values of 

t 1, the ASF+M model is faster than the AM model. The same can be said about integration points, n . 

This plot shows that for low Wc values (Wc/Wa  3), it is not advantageous to use the adaptive shape 

functions. The two models appear to have similar results for Wc/Wa = 3, and for Wc/Wa  3 using adaptive 

shape functions improves both the computational time and number of integration points needed. For 

Wc/Wa = 7, the converged solution was reached with the aid of adaptive shape functions in 1/30th the 

computational time. It is expected that the time savings would continue to increase with an increase in Wc. 

With increasing toughness of modern aerospace grade adhesives, it is expected that adaptive shape 

functions would be beneficial for many, if not most, adhesive bonding configurations. 

0

2

4

6

0 0.02 0.04 0.06 0.08 0.1

F (kN)

(mm) Stress

Strain

Adhesive Material

F, 

0

10

20

30

0 0.1 0.2 0.3 0.4 0.5

F (kN)

a) b) (mm) Stress

Strain

Adhesive Material

F, NASA/TM—2014-218333 17 

Figure 12.—Comparison of the computational time (t) and integration points 

(n) needed to reach convergence (e2  1 percent) for the joint element 

model with adaptive mesh (AM) and with both adaptive mesh and 

adaptive shape functions (ASF+M) for different stress/strain relations. 

Comparison With Cohesive Zone Model 

In order to compare the joint element model with more standard adhesive joint finite element 

techniques, a DCB and SLJ were modeled using the cohesive zone elements in Abaqus and the joint 

element models (ASF+M and AM). The coarsest finite element meshes for the Abaqus model and the 

ASF+M model are shown in Figure 13 for both joint configurations. The meshes shown represent the 

coarsest meshes that provided a converged solution (e1/e2  1 percent). Both Abaqus models contained 

2800 elastic, CPE4I elements, and 100 COH2D4 cohesive zone elements. The number of elements, 

number of integration points, and computational time for both adhesives are shown in Table 4, and the 

corresponding load versus displacement plots are shown in Figure 14. Computational time was not 

explicitly compared because the vast difference between an optimized commercial code and a flexible 

research code does not provide for a good comparison. However, the number of elements and total 

number of integration points should provide a good comparison. Since both joint configurations with the 

aforementioned adhesive had significant nonlinearity and adhesive cracking, both the adaptive shape 

functions (ASF) and adaptive mesh (AM) were beneficial in providing a converged solution with a 

minimum number of elements. 

When considering the adaptive shape functions, one might conclude that with enough integration 

points, the nonlinearity of any adhesive could be captured with one element. This may be true in theory, 

but numerical stability issues caused a requirement of more elements. With very long elements and highly 

non-linear tangent stiffness (especially negative stiffness), the chances increase that the matrix 

exponential used to solve Equation (10) becomes unstable and does not converge at large values of x. 

Therefore, more elements are needed to shorten the shape functions and keep them in the stable regime. 

Future work will include applying improved matrix exponential solution procedures to correct this 

problem. 

0

5

10

15

20

25

30

1234567

W/Wa

W



Wa



F, 

1

AM

ASF M

n

n

n 



AM

ASF M

t

t

t 

NASA/TM—2014-218333 18 

Figure 13.—Mesh used for the DCB (a, b) and SLJ (c, d) configurations for the 

Abaqus finite element simulations (a, c) and the ASF+M model (b, d). 

TABLE 4.—COMPARISON OF NUMBER OF ELEMENTS (nel), NUMBER OF INTEGRATIONS POINTS (n), 

AND COMPUTATIONAL TIME (t) FOR THE CONVERGED SOLUTION (e1/e2  1 PERCENT) 

FOR THE JOINTS WITH DIFFERENT MODELS 

Model Double cantilever beam Single lap joint 

el n n el n n

Adaptive mesh (AM) 320 5440 160 1360 

ASF + M 20 340 4 51 

Abaqus 2800 11400 2800 11400 

Figure 14.—Comparison of load versus displacement plots of different models for (a) DCB and (b) SLJ 

configurations. 

0

3

6

9

12

0 0.05 0.1 0.15 0.2

F (kN)

(mm)

0

0.2

0.4

0.6

0 0.5 1 1.5 2

AM

ASF+M

Cohesive Elements

F (kN)

(mm) a) b)

F, F, NASA/TM—2014-218333 19 

Summary and Conclusions 

This study deals with semi-analytical finite elements with shape functions which are determined

based on a model rather than prescribed, specifically bonded joint finite elements. Two methods were 

demonstrated to decrease the number of elements needed for a converged answer. 

The first is adaptive shape functions (ASF). Since the shape functions are determined from a model, 

the material properties are an input parameter affecting the shape functions directly. When the material 

softens, i.e., the material properties change, the shape functions can be updated based on the new softened 

properties. This causes the shape functions to adapt to the material, making them very representative of 

the actual structural behavior and allowing a single element to remain accurate even with highly nonlinear 

materials. 

The second method for decreasing the number of elements needed was using an adaptive internal 

mesh (AM). When a portion of adhesive was cracked within a joint element, the joint element was 

re-meshed with beam elements and a joint element. This re-meshing occurred internally within the 

element, requiring no alteration of the global model. This method can capture progressive failure of an 

adhesive with very few elements, but requires a nonlinear finite element model to be solved within 

cracked element. 

Joint element models with and without the two methods were compared with each other for a double 

cantilever beam and single lap joint configuration with brittle and ductile adhesives. It was shown that the 

adaptive shape function method became more beneficial with increasing ductility of the adhesive. With 

sufficient ductility, the time spent to recalculate the shape functions at every load increment became 

insignificant when compared to the reduction in mesh and integration points. The adaptive mesh method 

was beneficial for all cases with progressive and catastrophic material cracking. 

Finally, the joint element model with both methods was compared with a standard, built-in cohesive 

zone model in Abaqus. Using both adaptive shape functions and mesh, less than 3 percent of the elements 

and integration points were needed, which is a significant gain. 

Both methods have potential for broader applications within finite elements. The adaptive shape 

functions concept can be utilized whenever the shape functions have some sort of material or geometrical 

parameter with them to improve the accuracy of shape functions during changes in materials and 

geometry. Furthermore, the adaptive internal mesh can and has been used to capture the introduction and 

propagation of discontinuities within an element without requiring a global re-meshing of the model. 

References 

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

Challenges. Journal of Composites Technology and Research 2002;24(3):133. 

2. Samantha Bell. Structural adhesives deliver lighter composite structures with performance and costs 

benefits. Reinforced Plastics 2012. 

3. Volkersen O. Die Nietkraftverteilung in Zugbeanspruchten mit Konstanten Laschenquerschritten. 

Luftfahrtforschung 1938;15:41–47. 

4. Goland, Reissner E. The stresses in cemented joints. Journal of Applied Mechanics 1944;11:A17–

A27. 

5. Hart-Smith LJ. Adhesive-bonded double-lap joints. NASA CR-112235 1973. 

6. Hart-Smith LJ. Adhesive-bonded single-lap joints. NASA CR-112236 1973:116. 

7. Delale F, Erdogan F, Aydinoglu M n. Stresses in Adhesively Bonded Joints: A Closed-Form 

Solution. Journal of Composite Materials 1981;15(3):249–271. 

8. Mortensen F, Thomsen OT. Analysis of adhesive bonded joints: a unified approach. Composites 

Science and Technology 2002;62(7–8):1011–1031. 

9. Frostig Y, Thomsen OT, Mortensen F. Analysis of Adhesive-Bonded Joints, Square-End, and SpewFillet—High-Order Theory Approach. Journal of Engineering Mechanics 1999;125(11):1298–1307. 
NASA/TM—2014-218333 20 

10. Tsai MY, Morton J. An evaluation of analytical and numerical solutions to the single-lap joint. 

International Journal of Solids and Structures 1994;31(18):2537–2563. 

11. Adams RD, Peppiatt NA. Stress Analysis of Adhesive Bonded Tubular Lap Joints. The Journal of 

Adhesion 1977;9(1):1–18. 

12. Tsai MY, Oplinger DW, Morton J. Improved theoretical solutions for adhesive lap joints. 

International Journal of Solids and Structures 1998;35(12):1163–1185. 

13. Crocombe AD. Global yielding as a failure criterion for bonded joints. International Journal of 

Adhesion and Adhesives 1989;9(3):145–153. 

14. Harris JA, Adams RA. Strength prediction of bonded single lap joints by non-linear finite element 

methods. International Journal of Adhesion and Adhesives 1984;4(2):65–78. 

15. Bednarcyk B, Bansal Y, Collier C, Pindera M-J, Zhang J. Analysis Tools for Adhesively Bonded 

Composite Joints, Part 1: Higher-Order Theory. AIAA Journal 2006;44(1):171–180. 

16. Zhang J, Bednarcyk BA, Collier C, et al. Analysis Tools for Adhesively Bonded Composite Joints, 

Part 2: Unified Analytical Theory. AIAA Journal 2006;44(8):1709–1719. 

17. Zhang J, Bednarcyk B, Collier C, et al. 3D stress analysis of adhesively bonded composite joints. 

InVol 4.; 2005:2750–2773. 

18. Camanho PP, Tong L. Composite joints and connections: Principles, modelling and testing. 

Woodhead Publishing; 2011, 544 p. 

19. Towse A, Potter KD, Wisnom MR, Adams RD. The sensitivity of a Weibull failure criterion to 

singularity strength and local geometry variations. International Journal of Adhesion and Adhesives 

1999;19(1):71–82. 

20. Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crackpropagation analysis. International Journal for Numerical Methods in Engineering 1999;44(9):1267–

1282. 

21. Camanho PP, Dávila CG. Mixed-Mode Decohesion Finite Elements for the Simulation of 

Delamination in Composite Materials. Hampton, Virginia: Langley Research Center; 2002. 

22. Krueger R. Virtual crack closure technique: History, approach, and applications. Appl. Mech. Rev. 

2004;57(2):109–143. 

23. Xie D, Waas AM. Discrete cohesive zone model for mixed-mode fracture using finite element 

analysis. Engineering Fracture Mechanics 2006;73(13):1783–1796. 

24. Goutianos S, Sørensen BF. Path dependence of truss-like mixed mode cohesive laws. Engineering 

Fracture Mechanics 2012;91:117–132. 

25. Guiamatsia I, Ankersen JK, Davies GAO, Iannucci L. Decohesion finite element with enriched basis 

functions for delamination. Composites Science and Technology 2009;69(15–16):2616–2624. 

26. Guiamatsia I, Davies GAO, Ankersen JK, Iannucci L. A framework for cohesive element 

enrichment. Composite Structures 2010;92(2):454–459. 

27. Ostergaard MG, Ibbotson AR, Roux OL, Prior AM. Virtual testing of aircraft structures. CEAS 

Aeronautical Journal 2011;1(1-4):83–103. 

28. Eisenberger M, Yankelevsky DZ. Exact stiffness matrix for beams on elastic foundation. Computers 

& Structures 1985;21(6):1355–1359. 

29. Aydoan M. Stiffness-Matrix Formulation of Beams with Shear Effect on Elastic Foundation. 

Journal of Structural Engineering 1995;121(9):1265–1270. 

30. Gustafson PA, Waas AM. A bonded joint finite element for a symmetric double lap joint subjected 

to mechanical and thermal loads. International Journal for Numerical Methods in Engineering 

2009;79(1):94–126. 

31. tapleton SE, Waas AM. 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. 

32. Stapleton SE, Waas AM, Arnold SM. Functionally graded adhesives for composite joints. 

International Journal of Adhesion and Adhesives 2012;35:36–49. 
NASA/TM—2014-218333 21 

33. Stapleton SE, Waas A. Macroscopic Finite Element for a Single Lap Joint. In 

AIAA/ASME/ASCE/AHS/ASC 50th SDM Conference. Palm Springs, California; 2009. 

34. Rudraraju S, Salvi A, Garikipati K, Waas AM. Predictions of crack propagation using a variational 

multiscale approach and its application to fracture in laminated fiber reinforced composites. 

Composite Structures 2012;94(11):3336–3346. 

35. Rudraraju S, Salvi A, Garikipati K, Waas AM. Experimental observations and numerical simulations 

of curved crack propagation in laminated fiber composites. Composites Science and Technology 

2012;72(10):1064–1074. 

36. Scott E. Stapleton, Anthony M. Waas. The Analysis of Adhesively Bonded Advanced Composite 

Joints Using Joint Finite Elements. Cleveland, OH: NASA Glenn Research Center; 2012. 

37. Belytschko T, Hsieh BJ. Nonlinear transient finite element analysis with convected coordinates. 

International Journal for Numerical Methods in Engineering 1973;7(3):255–271. 

38. Crisfield MA. Non-Linear Finite Element Analysis of Solids and Structures. Wiley; 1996, 362 p. 

39. Crisfield MA, Moita GF. A unified co-rotational framework for solids, shells and beams. 

International Journal of Solids and Structures 1996;33(20–22):2969–2992. 

40. Guyan RJ. Reduction of stiffness and mass matrices. AIAA Journal 1965;3(2):380–380. 

41. Chen S, Pan HH. Guyan reduction. Communications in Applied Numerical Methods 1988;4(4):549–

556. 

42. Chandrupatla TR, Belegundu AD. Introduction to Finite Elements in Engineering. 2 Har/Dsk. 

Prentice Hall College Div; 1996,461 p.