Alcune informazioni rapide sugli elementi finiti, da elaborare come meglio credi per non essere più colto in castagna … Cercherò di essere molto discorsivo e poco cattedratico: sono convinto che si può parlare chiaramente senza perdere di vista il rigore scientifico, quando serve, e qui più di tanto non serve. L’idea di ‘elemento finito’ non è venuta nuova con i computer, ma prima dei computer nessuno sapeva bene cosa farsene. Era una bella curiosità – credo americana degli anni ’40 – che qualche cervellone s’era studiato a tavolino. E guai a pensare che serva solo nel campo strutturale: l’idea è nata come concetto matematico, non come applicazione specifica ai materiali. Ha numerosissime applicazioni Rapidissima premessa. Un solido è un solido, ma matematicamente assomiglia per certi aspetti anche ad un liquido, ad un gas o perfino al vuoto. Basta pensare che si parla di ‘spazi’ che hanno determinate proprietà fisiche. Così, un solido può essere uno spazio di cui ogni punto, ogni elementino elementare, deve osservare regole note per comportarsi bene con i vicini, per comportarsi “da solido”. Se smette di comportarsi bene, si rompe. Allo stesso modo, uno spazio pieno di gas non è diverso da un solido, solo che le regole degli elementini di gas sono diverse (sennò non sarebbe un gas, sarebbe un solido). Ancora, uno spazio vuoto in cui è presente un campo elettrico, o elettromagnetico, non è esattamente “vuoto”: gli elementini hanno alcune caratteristiche, alcune proprietà ed alcni obblighi nei confronti di quelli accanto. Ad esempio, l’obbligo di trasmettere una data grandezza fisica (che potrebbe essere un vettore elettrico) in un dato modo al proprio vicino. Torniamo a palla al nostro solido. Parliamo di strutture, cioè di entità che hanno capacità di reggere un carico meccanico. I solidi, meccanicamente parlando, per essere tali devono obbedire a ‘equazioni di congruità’ fra lo sforzo applicato e le deformazioni che ne derivano. Messo in altre parole, se ad un blocco di gomma appoggio venti chili, quello si schiaccia, si deforma e protesta, ma alla fin fine regge i venti chili. E li regge assumendo una forma (cioè, con una deformazione) che non è campata a caso, ma sempre quella, date quelle condizioni. Questo vale per tutti i solidi, ovviamente. Non vale per i liquidi e per i gas – che reggono un carico meccanico con altri meccanismi, ma non hanno capacità di assorbire deformazioni non avendo forma propria. Una tipica, tipicissima equazione di congruità è la legge di Hooke, o della deformazione proporzionale. Se ti attacco cinque chili, tu pezzo di ferro di allunghi d’un centimetro; se te ne attacco dieci, i centimetri saranno due. Questo è un caso un po’ particolare di deformazione, in buona misura ed a veder bene è un’ invenzione che semplifica i calcoli. Però alcuni solidi, come i metalli, la seguono molto da vicino. Una deformazione invece è elastica se, tolto il carico, la deformazione rispetto alla situazione precedente si annulla: le materie plastiche e quelle biologiche, tanto per dire, di solito non sono né completamente elastiche né, tantomeno, proporzionali. Esistono poi altre equazioni di congruità che legano, in buona sostanza, la posizione reciproca dei vertici di un cubetto di materiale, una volta noti i carichi applicati sui vertici stessi . Questo cubetto, che può anche essere un tetraedro (il cartoccio di latte di una volta) oppure un altro tipo di geometria, lo chiamiamo elemento finito . Così è scritta in forma dotta e va da sé che vale anche negli altri campi d’applicazione: il problema è semplissimo a descrivere, adesso. Mentre nella analisi strutturale tradizionale si inventavano simulazioni, approssimazioni e scorciatoie per far sembrare casi semplici le strutture complesse, e si introduceva inevitabilmente un errore, con gli EF o FE sembra d’avere l’uovo di Pasqua. Prendo l’oggetto, lo scompongo idealmente in tanti cubettini con i vertici collegati fra loro e comincia a
dire: se sul cubettino numero 1 ci sono due chili in quella direzione e si sposta di un millimetro in base alle equazioni di congruità, , sul numero 2 a fianco ci saranno 8 etti e si sposterà di due millimetri. Quindi sul numero 3 …. E così via. Questa è l’analisi ad elementi finiti: una splendida idea assolutamente inutile fino a che non sono stati inventati i computer per metterla in pratica. A che serve calcolare laboriosamente una dozzina di variabili per ogni cubetto, se hanno le dimensioni di un pacchetto di sigarette e non descrivono per nulla l’oggetto complicato che hai in mano ? Se non fai la rete dei cubetti abbastanza fitta e verosimile (si chiama mesh) non ottieni nessuna informazione intelligente .Per altro, mi ricordo benissimo che all’università, per dare corpo e sostanza al fondamentamentale ed illeggibile Brebbia-Connor, Advances in Finite Elements Techniques mi toccò affrontare un calcolo a mano sennò non mi facevano dare l’esame. Come dire, volendo si può fare, ma l’utilità è discutibilissima. Prendiamo ad esempio, giusto per figura, il pezzo a fianco. Chiaro che adesso coi computer è tutto più semplice: lo disegni e la descrizione matematica te la fa la macchina. Chiaro anche che la quantità di cubetti – magari non sono cubetti, anche piramidini, cilindretti e così via, ci sono equazioni già pronte per molti mattoncini elementari – la quantità dei calcoli da fare, dicevo, cresce spaventosamente con il crescere della complessità dell’oggetto. Immaginiamo solo per caso di avere un quadrato su un foglio, e dividerlo in quattro sotto-quadretti. Passiamo da quattro vertici a otto. Se dividiamo ancora a metà, i quadratini diventano sedici, ed i vertici 64. E siamo su un piano ! se fosse un cubo, a tre dimensioni, il numero delle variabili cresce alla follia. Ma tant’è, adesso si può fare anche questo. Un oggetto stupido come quello a fianco come ridere ha 5000 elementi di vario tipo, a losanga, a piramide, a tetraedro, a cubetto, a cilindro e così via, e venti o trentamila vertici da processare uno dopo l’altro per farli andare d’accordo in base alle equazioni di congruità, dei carichi e dei vincoli di quel particolare caso. Il computer prende la geometria, la spezzetta, numera vertici ed elementi, fa i calcoli … insomma, è una pacchia. Il problema è spiegargli cosa deve fare, e questo non è proprio semplicissimo. Non a caso ho preso questo oggetto, però. C’ un basamento che potrebbe essere metallico, con un ‘cuscino’ giallino ed un supporto che immagineremo di nuovo metallico. Ora, se imponiamo (è facile a scrivere, ma prova a convincerne il computer !) che i vertici dei cubetti che stanno fra il giallino e le altre parti, quelli di frontiera, se cadono nella parte giallina seguono un set di equazioni di congruità, ed un altro dall’altro lato, abbiamo descritto una giunzione incollata. Non è una
descrizione perfetta, sia chiaro, ma per i fini strutturali è proprio quello che serve: sapremo come di ripartiscono gli sforzi e le deformazioni nella frontiera fra il metallo e la colla. Se invec che metallico fosse materiale composito, il set di equazioni sarebbe diverso ancora, ma il concetto è sempre quello. Cioè, è possibile valutare con questo sistema l’efficacia del giunto nel suddividere i carichi, la capacità del giunto adesivo di trasmettere gli sforzi e misurare quanto bene funziona, che è la domanda che ti ha colto in castagna. Se sai quanti kili regge l’adesivo, puoi anche giudicare a spanne se il giunto regge oppure no, e con che affidabilità. Non è facile, sia chiaro. Innanzitutto perché esistono biblioteche intere sulle caratteristiche dei materiali più comuni, cioè sulle loro equazioni di congruità, ma sul comportamento degli adesivi i dati sono scarsi e poco pubblici. Spesso, quasi sempre, sono altamente non-lineari nelle deformazioni, magari pure un po’ viscosi e via dicendo. Poi, va da sé, si fa l’ipotesi che il giunto sia perfetto, cioè non abbia lui, per errata applicazione, delle discontinuità che non sono rappresentate nel modello. Sembra un dettaglio da nulla, ma è la ragione per cui ci si fida assai dei materiali metallici – altamente affidabili e controllabili – e piuttosto poco di quelli compositi o plastici come gli adesivi. E se ci fosse una bollicina d’aria ? Da ultimo, sempre solo per metter le cose in chiaro, in genere le analisi FE partono da ipotesi ben descritte: questi carichi, questi vincoli e così sia. Su certe strutture va bene, non c’è problema: non è verosimile che un ponte ferroviario sia usato come pista d’ormeggio per palloni aerostatici (peso negativo!) e se anche fosse, sposta poco rispetto. Ma chi è pronto a giurare che, per Dio! le forze di contatto in un giunto asimmetrico di una cerniera di portiera d’auto sono millimetricamente quelle e non altre, a prescindere dalla pioggia, dalla violenza dello sbattimento dell’amante che scende inviperita, e pure della lavatrice caricata nel baule perché tanto era in offerta ma cavolo! quanto pesa ? Come dire, dove casca l’asino nelle analisi FE non è tanto e soltanto nel descrivere la geometria (che è facile), i materiali (è difficile, ma se hai i dati sperimentali giusti ci stiamo ancora): è soprattutto nel descrivere i carichi ed i vincoli che effettivamente agiscono nelle varie condizioni. Veniamo ad alcune tecnicalità, utili per tamponare pressappoco eventuali domande. Il programma che uso io è sostanzialmente FEMAP, che poi è l’interfaccia grafica: il cervello è il NASTRAN che fu scritto dalla NASA negli anni Sessanta, e che è ancora sempre lo stesso insieme di regole, di sequenze di calcolo e algoritmi di allora . Difficile ormai trovare qualcosa di più collaudato. Credo che molti siano i programmi in giro - ad es c’è il CATIA della Dassault, tanto per buttare là qualche nome, oppure ABAQUS, ed ancora l’ANSYS – ma tutti all’incirca derivano da NASTRAN. Il che all’atto pratico fa poca differenza perché quel che li distingue sono più che altro la veste grafica, e quindi la minore o maggiore facilità di modifica senza usare comandi e tastiera. E poi la disponibilità di moduli speciali che fanno operazioni bizzarre: ad esempio
le analisi nel dominio del tempo, la simulazione di fenomeni termodinamici o elettromagnetici, i compositi, l’aeroelasticità e così via. Da ultimo, fra le differenze importanti dei vari programmi, c’è la capacità di generazione autonoma ed intelligente delle mesh: i programmi più seri ti permettono sempre di fare modifiche ed aggiustamenti alla mesh, quelli più stupidi ti propinano una soluzione che devi farti andare bene per forza. Quelli intelligentissimi hanno strategie di intelligenza artificiale per risolvere i problemi al posto tuo, e questi sono proprio programmi simpatici e molto costosi. Per chi li deve e vuole pagare, va da sé. La fondamentale differenza sta nella capacità e velocità di gestione dei dati: un personal di adesso ha capacità di calcolo mille volte superiore all’IBM360 di trent’anni fa, e forse anche i più. Vuol dire gestione della grafica, innanzitutto, ma anche modelli molto più fini e complessi. Un mainframe attuale, alimentato con buoni modelli, ti riproducono perfino le deformazioni plastiche di una carrozzeria sotto urto, o gli effetti di una esplosione termonuclerare. Tanto per dire, non si fanno più da decenni esperimenti nucleari e ben pochi test d’urto sulle macchine, più che altro per pubblicità visto quel che costano. E’ solo merito dei computer e del fatto che così si risparmia e si evitano polemiche. Apro una parentesi macabro-matematica a proposito dei test d’urto: per vedere se un’auto è veramente sicura, cioè, quale lesione creava un dato incidente sui passeggeri, fino a pochi anni fa si usavano veri cadaveri pieni zeppi di strumenti (accelerometri, sensori di torsione e di spostamento etc). Ti lascio immaginare i test dei seggiolini per i bambini, che problemi di organizzazione ponevano per la messa a punto. Di quelli morali non ne parla nessuno, si vede che non c’erano. Per non parlare delle donne incinte. Mettere a punto un airbag funzionante deve essere stato un incubo da film dell’orrore: mica per altro, reperimento dei cadaverini a parte: sostanzialmente perché i cadaveri non sono riproducibili in serie, uno si comporta così, l’altro cosà, e se fai una modifica al seggiolino che succederà? Film dell’orrore di un ingegnere. Secondo me inventavano alla grande, altro che. Poi hanno inventato i manichini androidi, i dummies, quelli con i cerchi a scacchi sulla testa: funzionano ma non è ancora proprio un sistema perfetto: dopo ogni urto van controllati e rincontrollati. Il passo successivo è un bel modello matematico di essere umano ad elementi finiti, che puoi distruggere duemila volte ed è sempre come nuovo. Metti la fibbia qui, e la cerniera là, fai tutte le prove che vuoi ! Un bel passo in avanti, non c’è che dire, ma oggi siamo solo alla affidabile riproduzione di ferite. A morire questi dummies matematici ancora non hanno provato, per ora per loro è troppo complicato. Tutti i programmi bene o male parlano fra loro, specialmente nell’importare le geometrie che hanno ormai protocolli standardizzati: IGES, Parasolid, e così via. Spesso poi si può prendere una stadio intermedio da un programma e ‘iniettarlo’ in un altro: ad es i dati di mesh calcolati da uno possono essere sfruttati in un altro software, i risultati passati ad un programma di presentazione e postprocessing ancora diverso e così via. Questo all’incirca credo sia il più che si può dire sulle analisi FE di giunti incollati.
Application of interface finite elements to three-dimensional progressive failure analysis of adhesive joints J . P . M . GO N CË A L V ES 1 , M . F. S . F. 1
DE
MOURA 1 , A . G. MA G A L H AÄ E S 2 a nd P. M . S . T .
DE
CASTRO1
Departamento de Engenharia MecaÃnica e GestaÄo Industrial, Faculdade de Engenharia da Universidade do Porto, Rua Dr Roberto Frias, 4200-465 Porto, Portugal, 2Departamento de Engenharia MecaÃnica, Instituto Superior de Engenharia do Porto, Rua Dr AntoÂnio Bernardino de Almeida, 431, 4200-072 Porto, Portugal Received in final form 17 October 2002
ABSTRACT
The paper presents a new model for three-dimensional progressive failure analysis of adhesive joints. The method uses interface elements and includes a damage model to simulate progressive debonding. The interface finite elements are placed between the adherents and the adhesive. The damage model is based on the indirect use of fracture mechanics and allows the simulation of the initiation and growth of damage at the interfaces without considering the presence of initial flaws. The application of the model to single lap joints is presented. Experimental tests were performed in aluminium/epoxy adhesive joints. Linear elastic and elastoplastic analyses were performed and the predicted failure load for the elastoplastic case agrees with experimental results. Keywords adhesive joints; damage models; interface finite elements; threedimensional analysis.
NOMENCLATURE
GIc, GIIc
INDICES
v displacement at a given point u vector of relative displacements D matrix of interface stiffnesses s interface stress tensor u0, i relative displacement corresponding to stress strength umax, i relative displacement corresponding to complete failure ui current relative displacement E damage parameters matrix I identity matrix and GIIIc critical fracture energies Sj current damage surface uj; i intersections of the damage surface with the coordinate axes Duj1, i relative displacement increment c unknown parameter in the damage surface definition s and t tangential directions n normal direction i direction s, t or n j current increment INTRODUCTION
Correspondence: M. F. S. F. de Moura. E-mail:
[email protected]
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
The well-known advantages of adhesive joints over mechanically fastened joints (less sources of stress concentrations, more uniform load distribution, better fatigue properties and fabrication without distortion) led to their widespread use in structural applications.
479
480
J . P . M . G O N CË A L V E S e t a l .
Consequently, many researchers have dedicated their attention to the mechanical and structural behaviour of those joints.1 However, most of the publications are related with stress calculations only.2±9 The failure criteria for adhesive joints can be based on strength of materials or on fracture mechanics. Those based on strength of materials generally use maximum stress or maximum strain failure criterion.6, 10, 11 The fracture mechanics criteria are based on the correlation of the energy release rate with joint fracture,12±14 or the use of a generalised stress intensity factor.15 Some authors proposed different techniques to study progressive failure of adhesive joints. Cheikh et al.16 present a model based on the finite element method formulated in displacements that enables stress conditions to be imposed. They use the method to impose the equilibrium condition and to model the stress vector continuity at the interfaces. The authors concluded that this method can be a highly valuable tool in design assistance and they intend to apply it to the study of fissuring and rupture in bonded joints. Pradhan et al.17, 18 analysed adhesive bonded joints using the finite element method by computing the strain energy release rate for debonding. They considered the interfaces of adherent and adhesive as potential sources for the initiation of the debonding that leads to failure. Paired nodes are generated on the interfaces and released one after another in a proper sequence, modelling gradual crack growth. Bonded joints were studied for various ratios of elastic moduli, thicknesses of adherent and adhesive, and overlap length for eight different crack locations along the interfaces. Failure loads were predicted and material properties and geometries that might result in extending the load-carrying capacity of these joints were suggested. The authors concluded that this approach provides an analytical tool for predicting crack initiation and growth that leads to failure. Bogdanovich and Yushanov19 present a new development of threedimensional variational stress/strain and progressive failure analysis of adhesive composite joints. The progressive failure analysis uses the strain energy release rate to predict different scenarios of cohesive, adhesive or interlaminar crack propagation. Through a numerical example of a double lap adhesive joint of unidirectional composites, the authors illustrated the application of their model. They first identify the site of initial failure and failure mode by analysing the computed strains. They then perform a progressive failure analysis starting with the initial crack(s), which can be located between the adhesive and the adherent, inside the adhesive, or between the layers of the composite adherent. They concluded that the model is applicable for predicting the dominating crack propagation path in the double lap composite bonded joint. The objective of this paper is to present a method that can predict the strength of single lap adhesive joints. The
method uses interface elements and includes a damage model to simulate progressive debonding. The model for predicting the adhesive joint strength consists of placing interface finite elements between the adherents and the adhesive. A decohesion failure criterion that combines aspects of strength-based analysis and fracture mechanics is used to simulate debonding by considering a softening process. The proposed constitutive equations for the interface are phenomenological mechanical relations between stresses and interfacial relative displacements. With increasing interfacial relative displacements, the stresses across the interface reach a maximum, decrease, and vanish when complete decohesion occurs. The work of normal and tangential separation can be related to the critical values of energy release rates, which is a fundamental task of the proposed damage model. One advantage of the method is that the onset and propagation of damage at an interface between different materials can be predicted without the previous consideration of an initial crack. Also, complex moving mesh techniques are not necessary to simulate damage growth. The application of the model to single lap joints is presented. Linear elastic and elastoplastic analyses were performed and the predicted failure load for the elastoplastic case agrees with experimental results. INTERFACE ELEMENT
In order to model the interfaces behaviour, an 18-node isoparametric interface element (see Fig. 1), which is compatible with 27-node three-dimensional solid elements implemented in the ABAQUSÕ software,20 was developed. The element consists of a zero-thickness volumetric element in which the interpolating shape functions for the top and bottom faces are compatible with the kinematics of the elements that are being connected to it. The interface finite element formulation, detailed in previous works21, 22 can be viewed as a contact problem. Basically, the solution of the equations arising 16 15 17
7
14
18 10
8
6 13
11 1
9 12
5
4
2 3 Fig. 1 The interface element.
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
APPLICATION OF INTERFACE FINITE ELEMENTS
in contact problems consists of minimizing the potential energy subjected to certain kinematic constraints. This is done considering the variational method with penalty function formulation. The vector of relative displacements between two homologous points can be obtained from the displacement fields of the element faces (top and bottom) 8 9 8 9 8 9 < us = < vs = < vs = v
1 u u t vt : ; : ; : t; un vn top vn bot
si
st, i
i = n, s, t k = I, II, III Gkc
where s and t represent the tangential directions and n the normal direction. The stresses resulting from the relative displacements, defined above, are given by s Du where 2
ds D40 0
2 0 dt 0
3 0 05 dn
3
being (ds , dt) the shear and (dn) the normal interface stiffnesses. They represent the penalty parameter (N m 3) introduced by the user and they have to be carefully chosen in order to obtain a good performance of the model. As small values induce large interpenetrations and large values induce numerical problems, the optimum interface stiffnesses are the largest values that do not produce numerical problems. DAMAGE MODEL
Pure modes (I, II or III) The need for an appropriate constitutive equation in the formulation of the interface element is fundamental for an accurate simulation of the damage growth. A constitutive equation is used to relate the stresses s to the relative displacements u at the interface. For pure mode I, II or III loading, a linear elastic-linear softening behaviour is used (see Fig. 2). A high initial stiffness is used to hold both faces of the interface element together in the elastic range. After the interfacial normal or shear stresses attained their respective tensile or shear strengths, the stiffness is gradually reduced to zero, which means that the failure does not occur instantaneously. It is assumed that failure energy is dissipated as the crack grows. The relative displacement umax, i , for which complete failure occurs, is obtained by equating the area under the softening curve to the critical fracture energy (GIc, GIIc or GIIIc), which is a material property. The bilinear interfacial constitutive response shown in Fig. 2 can be implemented as follows:
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
481
u 0, i
u max, i
ui
Fig. 2 Softening stress/relative displacements relationship for pure modes (I, II or III).
. For ui < u0, i the constitutive equation is given by (2) . For u0, i # ui # umax, i we have
s
I
EDu
4
where I is the identity matrix and E is a diagonal matrix containing on the i-position the damage parameter, ei
umax; i
ui ui
umax; i
u0; i u0; i
5
which represents the damage accumulated at the interface, which is zero initially and reaches one when the material is fully damaged. The other two elements of the diagonal matrix are equated to one simulating an immediate cancellation of the strength in those directions. . For ui > umax, i all the penalty stiffnesses are set equal to zero. If crack closure is detected, interpenetration is prevented by reapplying, only the normal stiffness (dn). Frictional effects are neglected.
The properties required to define the bilinear interfacial softening behaviour are the initial stiffness (penalty parameter) included in matrix D, the fracture energy release rate (GIc, GIIc or GIIIc) and the corresponding nominal interlaminar tensile or shear strengths. The accuracy of the results of the analysis depends on the penalty stiffness that is chosen. High values avoid interpenetration of the crack faces but can lead to numerical problems. After a substantial number of numerical experiments,23 it was concluded that a penalty stiffness of 107 N mm 3 produced converged results while avoiding potential convergence problems during the non-linear procedure.
482
J . P . M . G O N CË A L V E S e t a l .
To include the damage model in a conventional nonlinear numerical analysis, a tangential modulus matrix is required. This matrix can be obtained from the variation of Eq. (4) ds d
I
EDu
I
which leads to ds
I E
E
ED du
6
u0 D du u u0
7
being u0 the vector of critical relative displacements. The tangential modulus matrix is the operator that relates ds with du.
Mixed-mode model In adhesive joints, a more complex stress state is present and damage propagation is more likely to grow under mixed-mode loading (mode I and II). Therefore, a formulation for interface elements should deal with mixedmode damage growth problems. The mixed-mode damage considered22, 24 is based on the pure mode model described above. Firstly we consider two limiting surfaces (see Fig. 3): Initial Damage Surface (IDS) defined by the critical relative displacements (u0, s , u0, n) and the Final Damage Surface (FDS) defined by the maximum relative displacements (umax, s , umax, n). The position of the vector of resulting relative displacements (components s and n) leads to three different behaviours:
1 Up to IDS the stresses are linearly related with the relative displacements by Eq. (2). 2 Between the IDS and the FDS there is a softening process simulating a gradual failure. The position in the softening process is defined by the current damage surface, which is a function of the resulting displacement vector u j with components s and n. A definition of a new surface Sj1, at a point, is necessary whenever the resulting relative displacements pass the current surface Sj. In this case, the damage surface reaches a new position (Sj1) by translation and rotation of the previous surface (Sj). The translation and rotation of the surface are defined by the increments of the relative displacement in the two directions (see Fig. 3). The new surface Sj1 is then obtained from the current surface Sj and displacements increment doing uj1; s uj1; n 1
8 uj; s c uj1; s uj; n c uj1; n where uj; i (i s, n) are the intersections of the damage surface with the coordinate system and uj1; i are the relative displacement increments. The factor c is the smallest positive root of Eq. (9), which is obtained from Eq. (8) c2
uj1; s uj1; n c
uj1; n uj1; s uj1; s uj1; n
uj1; n uj; s
uj; n uj1; s uj; s uj1; n The new values uj1; i
uj; i
uj1; i
c
uj1; i
uj1; s uj; n
uj; s uj; n 0
9
are, uj; i :
10
un
u max, n
uj *+1, n
c ∆u j +1, n
FDS
uj,*n
u0, n
∆u j +1, n
(u j +1, s, u j +1, n) ∆u j +1, s
(uj, s,uj, n) IDS
Sj
u0, s
Sj +1
uj,*s
u j*+1, s
c ∆u j +1, s
u max, s
us
Fig. 3 Mixed-mode damage model.
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
APPLICATION OF INTERFACE FINITE ELEMENTS
The two damage parameters obtained from Eq. (11) are not equal because they depend on the current displacements in each direction. Consequently, they do not reach the value one at the same time. It was considered that, when one of these parameters reaches one, the stiffnesses are set equal to zero in the two directions. 3 After FDS the situation is similar to the pure modes model when ui > umax, i.
It must be noted that since mode III is not considered, the respective shear component is treated as in the pure mode model, i.e., the respective strength is immediately cancelled when IDS is attained. Some typical testing examples used in a previous work22 showed good performance of the model. The tangential modulus matrix can also be obtained from Eq. (7). However, in this case, the matrix E is a function of the intersections of the damage surface with the coordinate system (uj1; i ) and their variation can be obtained from Eq. (10) duj1; i c duj1; i dc
uj1; i
uj; i
12
where dc can be obtained by differentiating the valid root of Eq. (9).
Table 2 Mechanical properties of the epoxy adhesive (Terokal 5045) E (MPa)
s0.2 (MPa)
smax (MPa)
emax (%)
tmax (MPa)
GIc (N mm 1)
GIIc (N mm 1)
437.4
12.7
15.0
11.3
20
0.6
1.8
140 120 100 Stress, MPa
The diagonal elements of matrix E of Eq. (4) are now defined by umax; i
uj1; i u0; i ej1; i
11 uj1; i
umax; i u0; i
483
80 60 40 20 0 0
0.05
0.1
0.15
Strain Adhesive
Aluminium
Fig. 5 Stress±strain curves of the adhesive and the aluminium.
SINGLE LAP JOINT ANALYSIS
The geometry of the single lap joint considered is illustrated in Fig. 4. Three different overlap lengths were considered (see Fig. 4). The materials used are aluminium and epoxy adhesive, whose mechanical properties are presented in Tables 1 and 2. The stress±strain curves for both materials are shown in Fig. 5. Tensile tests were z y
x
Adhesive
A
Adherent
Adherent 25
87.5
B L
L 12.5 1.5 16.0 0.25 20.0 1.5
87.5
Fig. 4 Single lap joint with dimensions in mm.
Table 1 Mechanical properties of the adherent (Aluminium) E (MPa)
smax (MPa)
emax (%)
62088
124.1
14.0
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
Fig. 6 Single lap joint failure.
performed, by applying displacement at a rate of 1 mm min 1 until failure. The maximum recorded load appeared before the complete separation of the two adherents. It was observed that damage propagates near the interfaces (see Fig. 6). However, a thin adhesive layer was always observed on the adherent's surface, which means we are in presence of cohesive failure. A three-dimensional finite element model was considered to predict the failure of the single lap joint. This choice is supported by the fact that there are three-dimensional effects on the stress distributions at the adhesive, and at the interfaces between the adhesive and the adherents. This has been shown in previous studies,7, 25 and can also be seen in Fig. 7 where normalised normal and shear stresses at one of the interfaces are shown. Figure 7a shows the normal stress szz at the top interface for a joint with 12.5 mm of overlap length and under a load of 1.8 kN. In Fig. 7b, the shear stress txz for the same joint and loading condition is presented. Note
484
J . P . M . G O N CË A L V E S e t a l .
(a) 12
Normalised stress
10 8 6 4 2 0 −2 −4 15 10 y (mm)
5 0
4
2
0
6
8
10
12
14
x (mm)
(b) Normalised stress
8 7 6 5 4 3 2 1 0 15 10 y (mm)
5 0
0
2
4
6
8
10
12
14
adjacent to the interfaces have 27 nodes and the remaining are 20 node solid elements. Each adherent is modelled with two layers of solid elements across its thickness and one layer of solid elements represents the adhesive. The interface elements, which include the mixed-mode damage model, allow the simulation of damage growth. Since a cohesive failure was experimentally observed, the mechanical properties of the adhesive were considered in the damage model. Since the maximum load appears before total failure and considering that crack initiation occurs at critical regions (denoted by A and B in Fig. 4), corresponding to the maximum stress concentrations,16 no vertical interface elements were considered in the middle of the overlap length to simulate the total separation. Figure 9 shows a typical case of damage propagation numerically obtained. The mesh considered was constructed so that, at any moment of damage propagation, there are at least two interface elements undergoing the softening process. Previous studies22, 26 have shown that this was a sufficient condition to obtain accurate results. The single lap joints were fixed at one end and loaded (applied tensile displacements) in the longitudinal direction. Two different analyses were performed: linear elastic and elastoplastic behaviour of both materials. Typical load±displacement curves for experimental and numerical analyses are shown in Fig. 10 for 12.5 mm overlap length. As expected, the elastoplastic analysis leads to a
x (mm)
Fig. 7 (a) Normalised stress szz at the top interface between adhesive and adherent. (b) Normalised stress txz at the top interface between adhesive and adherent.
20 Width, mm
that the plane y 0 mm corresponds to the longitudinal symmetry plane and thus only half of the overlap is considered in Fig. 7a and b. In both figures, threedimensional effects are noticeable for the maximum stress at the end of the overlap length, which can be seen to decrease towards the free edge of the joint. This behaviour agrees with the numerically predicted crack front at start of failure as can be seen in Fig. 8. The mesh used in the finite element model consists of 1024 brick elements of 20 nodes (reference C3D2020), 384 brick elements of 27 nodes (reference C3D2720) and 256 interface elements. The connection between the two types of brick elements is possible because the C3D27 element has a variable number of nodes (between 21 and 27). Specifically, the nodes located at the centre of the faces can be considered or not. In order to connect this element to an element with 20 nodes, the face that is common to both elements does not contain the central node. The interface elements are placed at the interfaces between the adherents and adhesives. The solid elements
25
15 10 5 0 87
88
89
90
a, mm
Fig. 8 Numerically predicted crack front at start of failure in the elastic analysis of the 12.5 mm overlap length joint (87.5 to 100 mm).
Fig. 9 Damage propagation in a single lap joint.
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
APPLICATION OF INTERFACE FINITE ELEMENTS
(a) 4000 3500
Load, N
3000 2000 1500 500 0 0
2
4
6 a, mm
8
10
12
0
2
4
6 a, mm
8
10
Displacement, mm
(b) 2.5 2 1.5 1 0.5 0 12
(c) 4500
4000 3500 3000
4000 3500 3000 2500 2000 1500 1000 500 0 0
0.2
0.4
0.6
2500
0.8
1
1.2
1.4
1.6
a, mm
(d)
2000 1500 1000 500 0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Displacement, mm
Elastoplastic
Elastic
Experimental
Displacement, mm
Load, N
2500
1000
Load, N
smoother relation between applied displacement and reaction force and shows better agreement with experimental results. Figure 11a±d shows the evolution of the crack length a as a function of the load or the applied displacement for the joint with overlap length equal to 16 mm. Figure 11a and b corresponds to the elastic analysis and Fig. 11c and d is the result of the elastoplastic analysis. In the elastic case, it can be observed that there is some propagation before the failure load is attained followed by an abrupt propagation. For the elastoplastic analysis the propagation occurs before the maximum load is attained. After this, the applied displacement is absorbed by the plastification and the interface stress levels do not increase. In Table 3 a comparison between experimental and numerical results is presented. Each experimental result is the average of five tests. The numerical results for elastic materials agree for 12.5 mm overlap length. However, for longer overlaps the errors are higher. This can be explained by materials' plastification related to higher failure loads. This problem can be overcome considering an elastoplastic analysis, which globally presents good agreement with the experimental results.
485
1.2 1 0.8 0.6 0.4 0.2 0 0
Fig. 10 Load±displacement curves obtained from experimental and numerical analysis of 12.5 mm overlap joint.
Table 3 Experimental and numerical results for maximum adhesive joint load Overlap length (mm)
Experimental Numerical ± Elastic Materials Numerical ± Elastoplastic materials
Pmax (N) Pmax (N) Error % Pmax (N) Error %
12.5
16
20
3390 3385 0.25 3287 3.1
3544 3850 8.0 3550 0.17
4218 4641 10.0 4014 4.8
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
a, mm
Fig. 11 (a) Load±crack length a curve obtained in numerical elastic analysis of 16 mm overlap joint. (b) Applied displacement±crack length a curve obtained in numerical elastic analysis of 16 mm overlap joint. (c) Load±crack length a curve obtained in numerical elastoplastic analysis of 16 mm overlap joint. (d) Applied displacement±crack length a curve obtained in numerical elastoplastic analysis of 16 mm overlap joint.
CONCLUSIONS
A method for predicting the strength of single lap adhesive joints is presented. The numerical three-dimensional model includes a previously developed interface element and a mixed-mode damage model. These features allow
486
J . P . M . G O N CË A L V E S e t a l .
the simulation of crack initiation and growth, which are important issues in the failure process. Single lap joints were studied considering three different overlap lengths and two different analyses concerning the materials' behaviour. The comparison between the numerical and experimental failure loads showed that the linear elastic case only agrees for an overlap of 12.5 mm. Longer overlaps (16 and 20 mm) are related with higher failure loads which induce higher normal stresses contributing to the plastification of the components. These effects are not considered in the elastic analyses and, consequently, only elastoplastic analyses present good performance. The joints studied do not include spew fillets. For single lap joints with spew fillet, the crack path is more complicated. Studies will be done in order to adapt the model presented in this paper to such cases. Acknowledgement This project was finalised in the context of FCT project EME 35975/99 (POCTI/43525/EME/2000-FEDER). REFERENCES 1 Adams, R. D., Comyn, J. and Wake, W. C. (1997) Structural Adhesive Joints in Engineering. Chapman & Hall, London. 2 Goland, M. and Reissner, E. (1944) The stresses in cemented joints. J. Appl. Mech. 66, A17±A27. 3 Adams, R. D. and Peppiatt, N. A. (1974) Stress analysis of adhesive-bonded lap joints. J. Strain Anal. 9, 185±196. 4 Cooper, P. A. and Sawyer, J. W. (1979) A Critical Examination of Stresses in an Elastic Single Lap Joint. Report no. TP-1507, NASA. 5 Tsai, M. Y. and Morton, J. (1994a) An evaluation of analytical and numerical solutions to the single-lap joint. Int. J. Solids Struct. 31, 2537±2563. 6 Harris, J. A. and Adams, R. D. (1984) Strength prediction of bonded single lap joints by non-linear finite element methods. Int. J. Adhesion Adhesives 4, 65±78. 7 Tsai, M. Y. and Morton, J. (1994b) Three-dimensional deformations in a single-lap joint. J. Strain Anal. 29, 137±145. 8 Pandey, P. C. and Narasimhan, S. (2001) Three-dimensional nonlinear analysis of adhesively bonded lap joints considering viscoplasticity in adhesives. Computers Struct. 79, 769±783. 9 Andruet, R. H., Dillard, D. A. and Holzer, S. M. (2001) Two- and three-dimensional geometrical nonlinear finite elements for analysis of adhesive joints. Int. J. Adhesion Adhesives 21, 17±34. 10 Crocombe, A. D., Bigwood, D. A. and Richardson, G. (1990) Analyzing structural adhesive joints for failure. Int. J. Adhesion Adhesives 10, 167±178.
11 Czarnocki, P. and Piekarski, K. (1986) Non-linear numerical stress analysis of a symmetric adhesive bonded lap shear joint. Int. J. Adhesion Adhesives 3, 157±160. 12 Hamaush, S. A. and Ahmad, S. H. (1989) Fracture energy release rate of adhesive joints. Int. J. Adhesion Adhesives 9, 171±178. 13 Anderson, G. P., Brinton, S. H., Ninow, K. J. and DeVries, K. L. (1988) A fracture mechanics approach to predicting bond strength. In: Advances in Adhesively Bonded Joints. ASME, New York, pp. 93±101. 14 Fernlund, G., Papini, M., McCammond, D. and Spelt, J. K. (1994) Fracture load predictions for adhesive joints. Composites Sci. Technol. 51, 587±600. 15 Groth, H. L. (1985) A method to predict fracture in an adhesively bonded joint. Int. J. Adhesion Adhesives 5, 19±22. 16 Cheikh, M., Coorevits, P. and Loredo, A. (2001) Modelling the stress vector continuity at the interface of bonded joints. Int. J. Adhesion Adhesives 21, 249±257. 17 Pradhan, S. C., Iyengar, N. G. R. and Kishore, N. N. (1994) Parametric study of interfacial debonding in adhesively bonded composite joints. Composite Struct. 29, 119±125. 18 Pradhan, S. C., Iyengar, N. G. R. and Kishore, N. N. (1995) Finite element analysis of crack growth in adhesively bonded joints. Int. J. Adhesion Adhesives 15, 33±41. 19 Bogdanovich, A. E. and Yushanov, S. P. (1998) 3-D progressive failure analysis of bonded composite joints. In: AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 2, AIAA, 1616±1626. 20 ABAQUS User's Manual, Version 5.6. Hibbitt, Karlsson & Sorensen, Inc., Pawtucket (1996). 21 De Moura, M. F. S. F., GoncËalves, J. P. M., Marques, A. T. and De Castro, P. M. S. T. (1997) Modeling compression failure after low velocity impact on laminated composites using interface elements. J. Composite Mater. 31, 1462±1479. 22 GoncËalves, J. P. M., De Moura, M. F. S. F., De Castro, P. M. S. T. and Marques, A. T. (2000) Interface element including point-to-surface constraints for three-dimensional problems with damage propagation. Engineering Computations: Int. J. Computer-Aided Engng Software 17, 28±47. 23 De Moura, M. F. S. F., GoncËalves, J. P. M., Marques, A. T. and De Castro, P. M. S. T. (1996) Elemento finito isoparameÂtrico de interface para problemas tridimensionais. MeÂtodos NumeÂricos para CaÂlculo y DisenÄo en IngenierõÂa 12, 447±466. 24 Schipperen, J. H. A and De Borst, R. (1997) Mode-I and mixed mode delamination in laminated composites. In: Computational Plasticity (Edited by D. R. J. Owen et al.). CIMNE, Barcelona, 1237±1243. 25 GoncËalves, J. P. M., De Moura, M. F. S. F. and De Castro, P. M. S. T. (2002) A three-dimensional finite element model for stress analysis of adhesive joints. Int. J. Adhesion Adhesives 22, 357±365. 26 Mi, Y., Crisfield, M. A., Hellweg, H-B. and Davies, G. A. O. (1998) Progressive delamination using interface elements. J. Composite Mater. 32, 1246±1272.
ß 2003 Blackwell Publishing Ltd. Fatigue Fract Engng Mater Struct 26, 479±486