Adhesive bonding of aircraft primary structures has been in use for many years. For example, joining the stringers to skins of fuselage and wing structures, metallic honeycomb to the skins of elevators, ailerons, tabs, and spoilers constitute the main uses of adhesives in aircraft structures. Due to this increasing use of bonded structures in recent years, for weight saving, considerable work has been done in the fracture testing of different types of adhesive joints. However, most previous work on adhesive bonded joints deals mostly with Mode I fracture, and very little work appears to have been done on compressive delamination of adhesive joints. The Chow and Ngan test, consisting of two slender beams with a blister at the centre of each, which were bonded together and loaded in compression, is an exception. In the present work, a cohesive finite element model was developed for this blister test-piece, and the geometric non-linearity was incorporated in the strain/displacement relations. The crack propagation in the adhesive joint under compression for the proposed test-piece was found to agree with the available experimental observation. Most importantly, the delamination in adhesive joint under compression for different constrained cases was studied by cohesive finite element with different interface mesh densities. Stronger adhesive joint was achieved for less constrained Double Cantilever Beam (DCB) specimen.
Adhesive bonding has been used in the fabrication of primary aircraft fuselage and wing structures for many years . For example, joining the stringers to skins of fuselage and wing structures, metallic honeycomb to skins of elevators, ailerons, tabs, and spoilers constitute the main uses of adhesives in aircraft structures. Adhesive bonded aircraft structures are stable and durable, and, hence, this construction method has a good potential for future design programs. However, application of a new technology needs corresponding development of design and assessment methods. In industry, the earlier analytical design procedures have been replaced by the Finite Element Method (FEM) so that the complete structure including adhesive bonds can be simulated and assessed. The detailed FEM computations are based on detailed understanding of all relevant material behaviors including adhesives in joined or layered materials.
Delamination is one of the most common failure modes in layered materials, which may result from joint imperfections, edge effects, and various loadings. The presence of delamination can cause significant reduction in stiffness and strength of a joined structure and leads to a failure; hence a clear understanding of the failure behavior of the joined structure under shear/tension/compression is extremely important. Considerable work has been done in the compression delamination of composite structures [2-9], Mode-/ fracture of adhesive joints [10-14], and mixed mode fracture [15-21]; however, very little work appears to have been done on compressive delamination of adhesive joints.
Besides, different constrained ends of adhesive joints might change the fracture properties of the adhesive layer. In experimental study, capturing end constrained effects on fracture properties of the adhesive under compression requires proper complicated models of adherents and adhesive. On the other hand, the cohesive zone model can be introduced to describe the response of the adhesive layer and to simulate fracture process under compression for different constrained end cases.
The objective of this paper is, therefore, to study the delamination and growth of adhesively bonded joint under compression for different constrained end cases. Although most of the previous simulations of delamination growth were either for composite laminates or for the peel fracture of an adhesive joint, simulation for adhesively bonded joint under compression having different constrained ends is equally important. Chow and Ngan  proposed a test consisting of two slender beams with a blister at the centre of each and which were bonded together and loaded in compression. Simulation of this test along with those for the various constrained end effects on an adhesive joint under compression can be performed by using interface cohesive finite elements.
In the present work, a finite thickness interface cohesive element model is used to simulate the progressive delamination in adhesive joints under compressive in-plane loads for different constrained end conditions. For validation of the delamination model, a cohesive finite element model is developed for the blister test-piece, and the geometric non-linearity is incorporated in the strain/displacement relations. In addition, cohesive finite element model is validated with available literature . Finally, the delamination in adhesive joint under compression for different constraint cases is studied by cohesive finite element with different mesh densities of the interface.
Crack Propagation on Bonded Joint Under Compression
Chow and Ngan  proposed a test consisting of two slender beams with a blister at the centre of each and which were bonded together and loaded in compression. The geometry of the test-piece is depicted in Fig. 1, which consists of two slender beams bonded together with adhesive and loaded in compression. In Fig. 1, u is axial displacement, v is lateral displacement, a is the half crack length, q is quarter of the maximum distance between the curved beams, and P is the axial load. A blister arrangement is made at the centre of the beam simulating the bonding flaw, and the load displacement relationship is formulated.
Fig. 1 Slender beams with a blister at the center
The axial displacement, u, is given as
For small initial curvature and small deflection, both the (dy/dx)2 and (dv/dx)2 are much smaller than unity. From the first approximation of binomial expansion ,
where E is the modulus of elasticity, and / is the second moment of the area of the beam. Although the load displacement relationship is nonlinear for a constant a, for small P, Eq. (3) becomes linear.
The Cohesive Model
The idea for the cohesive model goes back to the strip yield models of Dugdale  and Barenblatt . Instead of letting stresses become singular, finite stresses are introduced in a cohesive zone ahead of at the crack tip, which are assumed to be same as the yield stress for a plane stress state  or as some function of the distance to the crack tip . This model has become very useful for practical applications, when numerical methods are used to solve nonlinear problems. The cohesive stresses or tractions have been introduced as functions of the local separation, of the material. This local separation is a vector having three components in mutually perpendicular directions for the selected orthogonal coordinate system. The cohesive model for crack propagation analysis of ductile materials is introduced by Needleman .
The interface element for cohesive finite element is an isoparametric element with a very small thickness as shown in Fig. 2. The thickness of the element is about 1/50 of the adherent thickness, and it is inserted as a numerical layer between the adherent layers. The strain vector, s, and the traction vector, t, are calculated in a local coordinate system (s,t,n), which is located on the element midplane. The local separation, 8, is calculated as the relative movement of the two surfaces from this midplane. The strain is then transformed to the global matrix, (x,y,z), to calculate the nodal reaction forces. The displacement vector, u, and the strain vectors are defined by
Fig. 2 Cohesive elements: a) 3D cohesive element, b) local coordinates direction, c) 2D cohesive element 
In an elastic region, if the deflection is large, nonlinearity is assumed to be for geometric nonlinearity and not due to a nonlinear stress-strain relationship. The strain displacement relationship for geometrically nonlinear case with respect to the reference coordinate system can be written as
Here i,j, and k run from 1 to 3, and Einstein’s summation role is used for repetitive indices. The indices 1, 2, and 3 represent the local s, t, and n directions respectively. The superscript ± represents the top/bottom surfaces with respect to the midplane. The coma between subscripts represents the partial derivative. From Eq. (6) replacing notation for local coordinate system, the strain displacement relationship for non-linear case is
where k runs from 1 to 3.
The material separation and thus damage to the structure is described by interface elements at the boundaries of the undamaged continuum elements. Thus, the mechanical behavior of the material is split into two parts, the damage-free continuum with an arbitrary constitutive law, and the cohesive zone specifying material damage and separation. Interface elements open according to some decohesion law and finally lose their stiffness so that the adjacent continuum elements get disconnected.
A bilinear constitutive law is chosen which is based on the 2-D cohesive model. This constitutive law relates the stresses to strains for any pure fracture or loading mode. In the case of mixed mode loading, the constitutive law relates the effective stresses to the effective strains, as shown in Fig. 3.
Fig. 3 Bilinear constitutive relationship
These quantities are defined as follows: The effective strain is a positive continuous quantity and when there is not any compressive normal strain, it is equal to the norm of strain vector and is defined by
where <> means if the inside parameter value of that is negative, it will be considered as zero. This constitutive law consists of three different parts :
i) If the effective strain is less than s0m (Fig. 3), the interface material behaves as linear elastic, and no damage is presented in the element.
ii) If the effective strain reaches s0m, the interlaminar damage initiates. After this point, the interface stresses decrease linearly.
iii) Strain /m refers to complete decohesion.
It is also necessary to specify the strains corresponding to initiation and completion of damage. The delamination initiation is predicted by the quadratic failure criterion. Damage is assumed to initiate when a quadratic interaction function involving the nominal stress ratios reaches a value of one. This criterion can be represented as
The constitutive equation for the elastic behavior can then be written as 
where n, s, and t are the local coordinate directions, and K is the stiffness matrix. This elasticity matrix (stiffness matrix, K) provides fully coupled behavior among all components of the traction and separation vectors. It is required to set the off-diagonal terms in the elasticity matrix to zero if uncoupled behavior between the normal and shear components is desired.
The Interface Properties
From Eq. (10), uncoupled constitutive behavior is defined as
The elasticity components for particular cases:means complete debonding between adherents,means perfect bonding between adherents. If the (interface thickness/adherent thickness) <<1, an approximation of stiffness values can be given by 
where E3, G13, and G23 are homogenized adherent moduli, and p is the interface thickness.
Cohesive Zone Modeling Validation
Double Cantilever Beam (DCB) Tests for Validation of Cohesive Zone Modeling
The DCB test specimen is generally used for the characterization of Mode-I fracture. The geometry and the loading conditions of the DCB configuration is shown in Fig. 4, and the finite element mesh is shown in Fig. 5. One layer of linear quadrilateral, type CPS4R (continuum plane stress 4 nodes reduced integration), elements are used for adherents, and one layer of linear quadrilateral, #508 COH2D4 (cohesive 2-dimensional 4 nodes), elements is used for the adhesive. The experiment consists of a load that is applied at the end blocks attached to the DCB specimen. The geometrical properties are the length L=203mm, the arm thickness h = 6.35mm, and width B = 25.4mm. The initial crack length a0 is about 55mm. The mechanical properties of the DCB specimen are E = 69 GPa, v = 0.3, GC = 1.6 N/mm .
Fig. 4 DCB configuration 
Fig. 5 Finite element model of the DCB configuration
The load history from the finite element simulation is compared to the experimental result of this DCB configuration tested at a cross head speed of 0.1mm/min. The corresponding load-displacement curves of the DCB tests are shown in Fig. 6. For the quasi-static loading, the model result was in excellent agreement with the experimental result .
Fig. 6 Load-displacement curves using the nominal interface strength 40 MPa for a DCB test
Compression Delamination Model
Since we have experimental results for compression delamination test for Aluminum adherents with initial crack at the center point along the longitudinal distance of the adhesive , we developed compression delamination model following the exact shape of the experimental set-up for proper validation of the model response. The finite element model of the system and the boundary and loading conditions is shown in the Fig. 7. The left end is fixed, and the x-directional displacement is applied at the right end. Three layers of linear quadrilateral, type CPS4R (continuum plane stress 4 nodes reduced integration), elements are used for Aluminum, and one layer of linear quadrilateral, type COH2D4 (cohesive 2-dimensional 4 nodes), elements is used for the adhesive. The geometrical and mechanical properties of the system: adherent thickness=6 mm, joint length on each side = 80 mm, E = 150 GPa, v = 0.25, and Gc = 0.352 N/mm.
Fig. 7 Finite element model of the compression delamination test
The linear relationship between the reaction force and free end displacement was observed (Fig. 8) before the crack started to propagate. Reaction force had linearly increased to a certain peak, where the crack started to propagate. The simulation result showed similar response to that given by Chow .
Fig. 8 Reaction force profile
Compression Delamination Under Different Constrained End Conditions
In this work, the DCB-specimen is studied under compression. The adhesive layer is modeled using cohesive elements. The effect of the compression and different constrained ends on the fracture is studied for different mesh densities of the cohesive interface. The geometry and the loading conditions of the DCB configurations are shown in Fig. 9 for fixed end and Fig. 10 for the top adherent with free end. The geometrical properties are the length Z=100mm; the arm thickness t1 = 2 mm, t2 = 4 mm, width B = 4 mm, and the initial crack length a0 = 50mm. The material properties for adherents are E1= E2=200 GPa, v = 0.3, and for adhesive are ultimate strength = 35 MPa, and fracture energy = 0.7 N/mm.
The result of the delamination under compression for constrained end for different mesh densities of the cohesive interface is shown in Fig. 11. The reaction force continued to increase to a certain peak before the delamination started to propagate. The reaction force gradually decreased after the peak, since the delamination occurred like the peel delamination. After the reaction force rises to a certain peak, crack propagates until the entire interface fails.The reaction force had increased to a certain peak before the delamination suddenly propagated. The reaction force rapidly decreased after the peak, since the delamination occurred in the entire interface, at the same time, after the maximum tolerable deformation. Since for the free end case, the entire layer failed at a certain peak (with no gradual decrease), the reaction force was larger in the free end case.
The delamination under compression in the adhesive joint was numerically investigated in this study. The delamination growth caused by compressive loads was modeled via softening behavior of cohesive interface elements. The effect of compression in crack propagation was studied and found to be in satisfactory agreement with experimental observations. A cohesive finite element model was developed for this blister test-piece, and the geometric non-linearity was incorporated in the strain/displacement relationship. The crack propagation through the adhesive joint under compression for the proposed test-piece was found to agree with the available experimental observation.
The effect of end constrained on the fracture resistance of the DCB specimen under compression was investigated. The numerical observations showed that for the fixed end DCB configuration, the reaction force had increased to a certain peak before the delamination started to propagate gradually. For the free end condition, on the other hand, the reaction force had increased to a certain peak before the delamination propagated rapidly. Since for the free end case the entire layer failed at a certain peak (with no gradual decrease), the reaction force was larger in the free end case.