Dinamica dell`impatto - Gruppo Italiano Frattura

UNIVERSITÀ DEGLI STUDI DI CASSINO
FACOLTÀ DI INGEGNERIA
TESI DI DOTTORATO IN INGEGNERIA CIVILE E MECCANICA
DINAMICA DELL’IMPATTO: INTERPRETAZIONE,
MODELLAZIONE E SIMULAZIONE NUMERICA DEL
COMPORTAMENTO MECCANICO DEI METALLI
Andrew Ruggiero
Università Degli Studi di Cassino
Facoltà di Ingegneria
Andrew Ruggiero
Dinamica dell’impatto: interpretazione,
modellazione e simulazione numerica
del comportamento meccanico dei metalli
Tesi di Dottorato in
Ingegneria Civile e Meccanica
XVIII ciclo
Coordinatore del corso:
Relatore:
Prof. Elio Sacco
Prof. Nicola Bonora
Indice Generale
SOMMARIO .............................................................................................................................................. 4
BIBLIOGRAFIA ......................................................................................................................................... 9
INDICE DELLE FIGURE...................................................................................................................... 10
INDICE DELLE TABELLE................................................................................................................... 14
1
INTRODUZIONE .......................................................................................................................... 15
2
ONDE DI SOLLECITAZIONE NEI SOLIDI............................................................................. 16
2.1
INTRODUZIONE ........................................................................................................................ 16
2.2
EQUAZIONE DELLE ONDE ......................................................................................................... 17
2.2.1
Tensione generata dall’impatto ......................................................................................... 18
2.2.2
Riflessione di onde elastiche alle interfacce ...................................................................... 19
2.2.3
Riflessione e trasmissione di onde elastiche in una discontinuità meccanica.................... 23
2.2.4
Tensione uniassiale............................................................................................................ 24
2.2.5
Deformazione uniassiale.................................................................................................... 27
BIBLIOGRAFIA ....................................................................................................................................... 33
3
MODELLAZIONE COSTITUTIVA............................................................................................ 34
3.1
INTRODUZIONE ........................................................................................................................ 34
3.2
MODELLI STRAIN RATE SENSITIVE ........................................................................................... 34
3.2.1
Modelli di resistenza formulati su basi fisiche................................................................... 36
3.2.2
Modelli di resistenza fenomenologici................................................................................. 38
3.2.3
Modello di resistenza di Johnson e Cook........................................................................... 38
3.3
3.3.1
MODELLI DI DANNEGGIAMENTO DUTTILE NEI METALLI .......................................................... 39
Modello di danno duttile non lineare................................................................................. 43
BIBLIOGRAFIA ....................................................................................................................................... 46
4
STRUMENTI DI SIMULAZIONE NUMERICA PER L’ANALISI DEI FENOMENI
DINAMICI ............................................................................................................................................... 48
4.1
INTRODUZIONE ........................................................................................................................ 48
4.2
ANALISI DINAMICA IN MSC.MARC ......................................................................................... 49
4.2.2
Houbolt Operator............................................................................................................... 52
4.2.3
Central Difference Operator.............................................................................................. 52
4.2.4
Damping............................................................................................................................. 53
4.3
ANALISI DINAMICA IN AUTODYN............................................................................................. 54
4.3.1
Metodo d’integrazione esplicito......................................................................................... 55
4.3.2
Viscosità artificiale ............................................................................................................ 56
2
4.4
IMPLEMENTAZIONE NUMERICA DEL MODELLO DI DANNO NON LINEARE .................................. 57
BIBLIOGRAFIA ....................................................................................................................................... 59
5
TAYLOR TEST ............................................................................................................................. 60
5.1
ANALISI TEORICA DEL TEST DI TAYLOR ................................................................................... 60
5.2
SIMULAZIONE NUMERICA DEL TAYLOR TEST ........................................................................... 62
5.3
ANALISI CRITICA DEI MECCANISMI DI PROPAGAZIONE DELLE ONDE DURANTE IL TEST DI
TAYLOR. ................................................................................................................................................ 69
5.4
5.4.1
ANALISI NUMERICA DEI MECCANISMI DI DANNEGGIAMENTO ................................................... 72
Effetto della dimensione del grano .................................................................................... 75
BIBLIOGRAFIA ....................................................................................................................................... 78
6
HOPKINSON BAR ........................................................................................................................ 80
6.1
PRINCIPIO DI FUNZIONAMENTO ................................................................................................ 80
6.2
SIMULAZIONE NUMERICA DELLA HOPKINSON BAR .................................................................. 83
6.2.1
Prova di compressione....................................................................................................... 83
6.2.2
Prova di trazione................................................................................................................ 90
6.3
CONCLUSIONI ........................................................................................................................ 104
BIBLIOGRAFIA ..................................................................................................................................... 105
7
FLYER PLATE IMPACT TEST................................................................................................ 106
7.1
SIMULAZIONE NUMERICA DEL FLYER PLATE IMPACT TEST ................................................... 109
7.2
ANALISI DELLO “SPALL SIGNAL” .......................................................................................... 113
7.2.1
Modello numerico ............................................................................................................ 115
7.3
EFFETTI GEOMETRICI SUL PROCESSO DI FRATTURA PER SPALLING ......................................... 117
7.4
RE-SHOCK EXPERIMENT......................................................................................................... 123
7.4.1
Fenomenologia del re-shock............................................................................................ 125
BIBLIOGRAFIA ..................................................................................................................................... 132
8
CONCLUSIONI ........................................................................................................................... 133
3
Sommario
Nel presente lavoro di tesi si è analizzata la risposta meccanica dei metalli in condizioni
d’impatto veloce. A tale scopo, si sono utilizzati gli strumenti della simulazione
numerica per analizzare tre configurazioni sperimentali largamente diffuse per la
caratterizzazione della risposta meccanica dei materiali in regime dinamico: il Taylor
Test, la Hopkinson Bar e il Flyer Plate Impact Test.
Ad oggi, il limite maggiore nell’utilizzo dei codici numerici come strumenti di
previsione del comportamento dei componenti meccanici in condizione d’impatto
veloce è dato dalla scarsa disponibilità di modelli costitutivi in grado di descrivere il
comportamento dei materiali in regime dinamico. Tali modelli devono essere in grado
di tenere in conto gli effetti, sulla resistenza del materiale, della deformazione, della
velocità di deformazione, della temperatura, del danneggiamento e, per impatti
iperveloci, della pressione idrostatica. Nel presente lavoro si sono utilizzati due modelli
indipendenti per trattare in modo disaccoppiato gli effetti dovuti ai parametri citati. Gli
effetti della deformazione, della velocità di deformazione e della temperatura sulla
resistenza del materiale sono stati descritti con il modello fenomenologico di Johnson e
Cook, [1], mentre gli effetti del danneggiamento sono stati descritti con un modello di
danno non lineare per rottura duttile nei metalli, [2]. Il modello costitutivo così
composto è stato implementato nei codici commerciali, MSC.Marc e Autodyn, utilizzati
per le simulazioni numeriche.
Il Taylor Test, [3], consiste nel far impattare, a velocità nota, un provino di forma
cilindrica contro una parete rigida. Il valore della tensione di snervamento in regime
dinamico è correlato, attraverso una semplice analisi monodimensionale, alla velocità
d’impatto e alla deformata del provino. Nel corso degli anni molti lavori hanno avuto
come obiettivo quello di superare alcune delle limitazioni date dalle ipotesi, di seguito
elencate, alla base della teoria di Taylor:
‰ la propagazione delle onde all’interno del cilindro sia monodimensionale;
‰ il materiale abbia un comportamento rigido perfettamente plastico ed
indipendente dalla velocità di deformazione, σ=σ(ε);
‰ il flusso plastico sia incompressibile;
4
‰ la deformazione elastica sia trascurabile.
Nel 1954, Lee e Tupper, [4], presentarono una modifica alla formulazione di Taylor per
tenere includere nell’analisi la deformazione elastica. Raftopoulos e Davis, [5],
inclusero la deformazione elastica e l’incrudimento. Jones et al. [6], proposero una
nuova equazione del moto per la parte indeformata del provino. Nel 1981, Erlich et al.
[7], proposero una tecnica alternativa, denominata “Rod on Rod” (ROR), in cui
l’impatto avviene tra due cilindri di medesimo materiale e uguale diametro, in modo da
eliminare le incertezze derivanti dalla mancata conoscenza delle condizioni di attrito
delle superfici a contatto.
Il punto più critico di tale approccio, comunque, sta nel fatto che la condizione di
unidimensionalità dello stato di sforzo non è verificata. È infatti notorio che le onde di
rilascio, che dal bordo della superficie d’impatto propagano radialmente, si
sovrappongono in prossimità dell’asse di simmetria, dando luogo ad un’onda di
tensione.
L’analisi del test, per mezzo degli strumenti della simulazione numerica e del modello
costitutivo implementato, ha portato all’individuazione di due diversi modi in cui i
meccanismi di danneggiamento duttile possono aver luogo. Per valori elevati della
deformazione di soglia il danno è causato da grandi deformazioni plastiche, in stato di
bassa triassialità dello stato di sforzo, che avvengono in prossimità della zona di
contatto, tardi nel processo di deformazione. Per valori della soglia di deformazione
relativamente basse il danno, al contrario, si sviluppa con un basso livello di
deformazione plastica, ad elevata triassialità dello stato di sforzo, nelle prime fasi del
processo di deformazione
Si è inoltre individuata una relazione tra il valore della deformazione di soglia, che è
uno dei parametri del modello di danno, e la dimensione media del grano.
La Hopkinson Bar è, ad oggi, la tecnica sperimentale più utilizzata per la
caratterizzazione della risposta meccanica del materiale in regimi di velocità di
deformazione che vanno da 102 a 104 s-1. Anche per questa tecnica l’ipotesi
fondamentale è che lo stato di sforzo possa essere assunto uniassiale.
Nel presente lavoro si è verificato che tale ipotesi può essere considerata vera sia per la
classica configurazione a compressione della Split Hopkinson Pressure bar, sia per una
5
configurazione alternativa, proposta da Staab e Gilat, che permette di effettuare la prova
direttamente in trazione. Si è verificato, inoltre, che con un’attenta progettazione della
prova, ovvero della determinazione della geometria del provino in relazione alla
geometria delle barre e alle impedenze meccaniche dei materiali utilizzati, la
deformazione e la velocità di deformazione possono essere ritenute, con buona
approssimazione, uniformi all’interno del provino. Questo è di importanza rilevante, in
quanto, permette di accettare i risultati ottenuti con questa tecnica sperimentale come
identificativi del comportamento meccanico del materiale in regime dinamico.
Si è verificato inoltre che i limiti maggiori in tale tecnica sperimentale sono dati dalla
difficoltà di mantenere costante, durante l’intera durata della prova, la velocità di
deformazione. In particolare si è osservato come nella prova di compressione la velocità
di deformazione decresce a causa dell’aumento di sezione per effetto Poisson. Nella
prova a trazione, al contrario, è proprio l’effetto di strizione a consentire il
mantenimento di un’elevata velocità di deformazione ad un valore pressoché costante.
A valle di tali verifiche si è provveduto alla progettazione e realizzazione di una barra di
Hopkinson a trazione in grado di caratterizzare il comportamento meccanico in regime
dinamico di materiali metallici quali il rame, gli acciai, alcune leghe di nichel, quali il
waspaloy etc.
La terza configurazione analizzata è quella del Flyer Plate Impact Test, che è l’unica
tecnica sperimentale che permette la caratterizzazione meccanica dei materiali per
velocità di deformazione superiori a 104 s-1.
La tecnica consiste nel realizzare un impatto planare, a velocità nota, tra due dischi
sottili. Un rapporto elevato tra il diametro dei dischi e il loro spessore (D/h>10)
garantisce uno stato di deformazione uniassiale in prossimità dell’asse di simmetria dei
dischi.
È
importante
sottolineare
che,
diversamente
dalle
configurazioni
precedentemente illustrate nelle quali si cercava, con maggiore o minore successo, di
realizzare uno stato sforzo che potesse essere assunto con buona approssimazione
unidimensionale, nel Flyer Plate Impact Test si verifica, a tutti gli effetti, uno stato di
deformazione unidimensionale.
Per tale tecnica, quindi, è disponibile, in forma esatta, una trattazione teorica che può
essere utilizzata per la verifica ed il confronto con i risultati numerici.
6
Tale test è largamente utilizzato anche per il particolare tipo di rottura che è in grado di
produrre nel disco bersaglio. La rottura, denominata spalling, avviene per una trazione
localizzata provocata dalla sovrapposizione dell’onda riflessa dalla superficie libera del
target e della sopraggiungente onda di rilascio.
Nell’esperimento, la misura avviene mediante la rilevazione, ad esempio attraverso
tecniche d’interferometria laser, del profilo di velocità di un punto situato sulla
superficie posteriore del disco bersaglio. La lettura del profilo permette di ricavare tutte
le informazioni necessarie ad identificare il comportamento meccanico del materiale.
Nel presente lavoro si è, dapprima, verificato che il modello numerico realizzato
riuscisse a riprodurre tutte le caratteristiche chiave dell’esperimento quali: i tempi di
arrivo del precursore elastico e dell’onda plastica sulla superficie libera del bersaglio;
l’intensità e la durata del plateau di massima velocità; i tempi di arrivo delle onde di
rilascio, elastica e plastica; l’arrivo dell’onda generata dalla creazione della superficie
libera dovuta alla rottura per spall. Si è poi passati all’analisi e all’interpretazione dei
fenomeni di deformazione e rottura del processo d’impatto che ha permesso di
raggiungere i seguenti risultati.
Si è identificato un processo di dissipazione nel meccanismo di separazione delle
superfici di rottura che porta ad una discordanza tra lo “spall signal” calcolato e quello
misurato sperimentalmente.
Si è effettuato uno studio parametrico degli effetti geometrici associati al processo di
rottura per spall, in impatti planari. Tale studio ha permesso di individuare un criterio
geometrico per la valutazione delle condizioni di spalling e per la determinazione della
posizione di primo innesco.
Si è infine analizzata una configurazione sperimentale che consiste nel posizionare sulla
parte posteriore del disco proiettile, un disco di maggiore impedenza meccanica. In tale
modo, al posto dell’onda di rilascio, si genera un’onda di compressione che
sovrapponendosi all’onda di compressione generata dall’impatto, porta al fenomeno del
re-shock.
Secondo la teoria delle onde di sforzo, la risposta di un materiale elasto-plastico, al reshock dovrebbe essere interamente plastica, poiché lo stato del materiale si dovrebbe
trovare sulla superficie di snervamento. Dagli esperimenti si osserva invece la presenza
7
di un gradino che precede l’arrivo del ricaricamento plastico, comunemente
riconosciuto come un inaspettato precursore plastico. Le interpretazioni del fenomeno
presentate in letteratura sono tutte basate sull’ipotesi che meccanismi fisici, che hanno
luogo alla micro o meso scala, portano lo stato del materiale all’interno della superficie
di snervamento.
In questo lavoro si è data un’interpretazione alternativa del fenomeno, basata su
considerazioni alla macroscala, capace di giustificare la presenza del gradino e di
dimostrare che non è un precursore elastico. Tale interpretazione, inoltre, chiarisce
anche la ragione per cui la ricompressione ed il corrispondente rilascio non debbano
essere, come dimostrato dagli esperimenti, sincroni.
8
Bibliografia
[1]
Johnson, G. R. and Cook, W. H., A constitutive model and data for metals
subjected to large strains, high strain rates and high temperatures, Proc. 7° Int.
Symp. On Ballistics, pp. 541-547, Netherlands, 1983.
[2]
Bonora, N., (1997), Int. J. of Fracture, 88, 359-371.
[3]
Taylor, G. I., “The use of flat ended projectiles for determining dynamic yield
stress: 1. Theoretical considerations” in Proc. R. Soc. Lond. Ser. A., vol 194, pp.
289-300, 1948.
[4]
Lee, E., and Tupper, J. Appl. Mech., 63-70 (1954).
[5]
Raftopoulos, D., and Davids, N., AIAA J. 5, 2254 (1967).
[6]
Jones, S. E., Gillis, P. P., and Foster, J. C., Jr., J. of Appl. Phys. 61, 499-502
(1987).
[7]
Erlich, D. C., Shockey, D. A., and Seaman, L., “Symmetric Rod Impact
Technique for Dynamic Yield Determination”, in Shock Waves in Condensed
Matter-1981, AIP Conference Proceedings 78, Menlo Park, CA, 1981, pp. 402406.
9
Indice delle Figure
FIGURA 2.1 - FORZE AGENTI SULL’ELEMENTO DI MASSA, [5]...................................................................... 17
FIGURA 2.2 – SCHEMA DI UN CORPO RIGIDO CHE IMPATTA UN CILINDRO A VELOCITÀ v0 , [5]. ................... 19
FIGURA 2.3 – RIFLESSIONE DI UN’ONDA LONGITUDINALE SU UNA SUPERFICIE LIBERA. .............................. 20
FIGURA 2.4 – RIFLESSIONE DI UN’ONDA LONGITUDINALE SU UNA SUPERFICIE LIBERA ORTOGONALE ALLA
DIREZIONE DI PROPAGAZIONE DELL’ONDA, [1].................................................................................. 20
FIGURA 2.5 – RIFLESSIONE DI UN’ONDA LONGITUDINALE SU UNA SUPERFICIE FISSA ORTOGONALE ALLA
DIREZIONE DI PROPAGAZIONE DELL’ONDA, [1].................................................................................. 22
FIGURA 2.6 – RIFLESSIONE E TRASMISSIONE DI UN’ONDA AD UNA DISCONTINUITÀ MECCANICA. ............... 23
FIGURA 2.7 – CURVA TENSIONE DEFORMAZIONE E PROFILI D’ONDA PER UN MATERIALE BILINEARE. ......... 25
FIGURA 2.8 – CURVA TENSIONE DEFORMAZIONE E PROFILO D’ONDA PER UN MATERIALE ELASTO-PLASTICO,
SECONDO LA “RATE INDEPENDENT THEORY”.
................................................................................... 26
FIGURA 2.9 – CURVA SFORZO-DEFORMAZIONE PER I MATERIALI ELASTICO PERFETTAMENTE PLASTICO ED
ELASTICO CON INCRUDIMENTO LINEARE, IN CASO DI STATO DI SFORZO UNIASSIALE OVVERO DI
DEFORMAZIONE UNIASSIALE, [1]....................................................................................................... 29
FIGURA 2.10 – CURVA SFORZO-DEFORMAZIONE IN STATO DI DEFORMAZIONE UNIASSIALE, PER VALORI
DELLA PRESSIONE ESTREMAMENTE ELEVATI, [5]. ............................................................................. 31
FIGURA 3.1 – EFFETTO DELLA VELOCITÀ DI DEFORMAZIONE SULL’ALLUMINIO COMMERCIALE CARICATO A
TAGLIO.............................................................................................................................................. 35
FIGURA 3.2 - EFFETTO DELLA TEMPERATURA SUL TITANIO α. .................................................................... 35
FIGURA 3.3 - VARIAZIONE DEL VALORE DELLO SNERVAMENTO DI UN ACCIAIO BASSO LEGATO CON LA
VELOCITÀ DI DEFORMAZIONE E LA TEMPERATURA............................................................................ 37
FIGURA 3.4 - RIDUZIONE DELLA DUTTILITÀ AL CRESCERE DELLA TRIASSIALITÀ DELLO STATO DI SFORZO
PER L'ACIAIO SA537, [14]. ................................................................................................................ 41
FIGURA 3.5 - EVOLUZIONE DEL DANNO, NORMALIZZATO RISPETTO AL DANNO CRITICO, IN FUNZIONE DELLA
DEFORMAZIONE PLASTICA, PER DIVERSI TIPI DI METALLI. ................................................................. 44
FIGURA 4.1 – DISTORSIONE DI UNA MESH LAGRANGIANA, [3]. ................................................................... 55
FIGURA 4.2 – TIPICA PROCEDURA DI “REZONING”, [3]................................................................................ 55
FIGURA 4.3 - SCHEMA LOGICO PER IL CALCOLO NUMERICO DELLA VARIABILE DI DANNO. LA PROCEDURA
INIZIA ALLA FINE DI OGNI INCREMENTO, QUANDO SONO GIÀ STATE CALCOLATE TUTTE LE VARIABILI
GLOBALI E ED È RIPETUTA PER OGNI PUNTO DI GAUSS DI OGNI ELEMENTO ATTIVO............................ 58
FIGURA 5.1 - SCHEMATIZZAZIONE DEL CILINDRO DI TAYLOR: (A) DURANTE LA DEFORMAZIONE; (B) AL
TERMINE DELLA DEFORMAZIONE.
..................................................................................................... 60
FIGURA 5.2 - DETTAGLIO DELLA MESH NELLA ZONA D'IMPATTO, PER IL RAME OFHC. .............................. 63
FIGURA 5.3 - ELEMENTO SEMINFINITO USATO NELLA MODELLAZIONE DELL'INCUDINE. ............................. 64
FIGURA 5.4 - PROFILO DELA DEFORMATA PER IL RAME OFHC (L0=25,4MM; V0=190M/S) OTTENUTA
ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* ) . .............................................................................. 66
10
FIGURA 5.5- PROFILO DELLA DEFORMATA PER IL FERRO ARMCO (L0=12,6MM; V0=279M/S) OTTENUTA
ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* ) . .............................................................................. 66
FIGURA 5.6 - PROFILO DELLA DEFORMATA PER L’ACCIAIO (L0=8,1MM; V0=343M/S) OTTENUTA ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* ) . ................................................................................................... 67
FIGURA 5.7 - PROFILO DELLA DEFORMATA PER IL RAME OFHC (L0=25,4MM; V0=190M/S) OTTENUTA
ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) . ............................................................. 67
FIGURA 5.8 - PROFILO DELLA DEFORMATA PER IL FERRO ARMCO (L0=12,6MM; V0=279M/S) OTTENUTA
ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) . ............................................................. 68
FIGURA 5.9 - PROFILO DELLA DEFORMATA PER L’ACCIAIO (L0=8,1MM; V0=343M/S) OTTENUTA ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) . .................................................................................. 68
FIGURE 5.10 A, B, C E D - GENERAZIONE, PROPAGAZIONE E SOVRAPPOSIZIONE DELLE ONDE DI PRESSIONE IN
UN ROR TEST A DIVERSI ISTANTI DI TEMPO DURANTE IL PROCESSO DI DEFORMAZIONE.
................... 71
FIGURA 5.11 - MAPPA DI DANNO PER IL RAME OFHC (L0=25,4MM; V0=190M/S) OTTENUTA ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* ) . ................................................................................................... 73
FIGURA 5.12 - DIAGRAMMI TENSIONE DEFORMAZIONE, AL VARIARE DELLA LEGGE COSTITUTIVA, PER UN
PUNTO APPARTENENTE ALLA SUPERFICIE DI CONTATTO PER L’IMPATTO DI UN CILINDRO DI FERRO
ARMCO. .......................................................................................................................................... 73
FIGURA 5.13 - MAPPA DI DANNO PER IL RAME OFHC (L0=25,4MM; V0=190M/S) OTTENUTA ASSUMENDO
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) . .................................................................................. 74
FIGURA 5.14 - DEFORMAZIONE DI SOGLIA IN FUNZIONE DELLA DIMENSIONE MEDIA DEL GRANO. .............. 76
FIGURA 5.15 - DEFORMATE E MAPPE DI DANNO CALCOLATE A CONFRONTO CON I RISULTATI SPERIMENTALI.
.......................................................................................................................................................... 77
FIGURA 6.1 - SCHEMATIZZAZIONE DELL'APPARATO E DELLA STRUMENTAZIONE DI UNA CONFIGURAZIONE
CLASSICA DELLA HOPKINSON IN COMPRESSIONE, [1]........................................................................ 80
FIGURA 6.2 - SCHEMATIZZAZIONE DEGLI IMPULSI DI DEFORMAZIONE ALLE INTERFACCE BARRE PROVINO
[1]. .................................................................................................................................................... 81
FIGURA 6.3 - SCHEMA DEL SISTEMA DI PROVA DELLA BARRA DI HOPKINSON A COMPRESSIONE SIMULATO
AGLI ELEMENTI FINITI. ...................................................................................................................... 83
FIGURA 6.4 - PARTICOLARE DELL'INTERFACCIA TRA BARRE E PROVINO NELLA CONFIGURAZIONE NON
DEFORMATA INIZIALE.
...................................................................................................................... 84
FIGURA 6.5 - ANDAMENTI DELLE TENSIONI DURANTE LA PROVA DI COMPRESSIONE................................... 85
FIGURA 6.6 - ANDAMENTI TEMPORALI DEGLI STRAIN RATES PER DIVERSE VELOCITÀ D’IMPATTO. ............. 87
FIGURA 6.7 - DIAGRAMMI TENSIONE-DEFORMAZIONE PER LE DIVERSE VELOCITÀ DI CARICO..................... 87
FIGURA 6.8 - DISTRIBUZIONE DELLE TENSIONI SULL’ASSE DI SIMMETRIA DELLE BARRE DI PRESSIONE A
10 µ s DALL’ISTANTE IN CUI È AVVENUTO L’IMPATTO...................................................................... 88
11
FIGURA 6.9 –ONDE DI DEFORMAZIONE INCIDENTE E RIFLESSA REGISTRATE SUGLI ESTENSIMETRI DOPO
CIRCA 320 ΜS DALL’IMPATTO. .......................................................................................................... 89
FIGURA 6.10 – CONFRONTO TRA UN PROVINO CILINDRICO NON DEFORMATO (A) ED I PROFILI FINALI
DEFORMATI CON BARRELING (B) E SENZA (C).................................................................................... 90
FIGURA 6.11 - SCHEMA FUNZIONALE DL DISPOSITIVO DI STAAB E GILAT................................................... 90
FIGURA 6.12 – PARTICOLARE DELLA MESH ADOTTATA PER LA DISCRETIZZAZIONE DEL PROVINO DI RAME
PURO. ................................................................................................................................................ 92
FIGURA 6.13 - ANDAMENTI DELLE ONDE DI TENSIONE REGISTRATE DURANTE LE PROVE............................ 93
FIGURA 6.14 - ANDAMENTI TEMPORALI DELLE VELOCITÀ DI DEFORMAZIONE RAGGIUNTE......................... 94
FIGURA 6.15 - CURVE TENSIONE-DEFORMAZIONE OTTENUTE PER DIVERSE VELOCITÀ DI CARICO............... 95
FIGURA 6.16 - MISURE DEGLI STRAIN RATES EFFETTUATE CON ESTENSIMETRI DI LUNGHEZZA VARIABILE. 96
FIGURA 6.17 – DISTRIBUZIONE DEL DANNO SULLA DEFORMATA FINALE DEL PROVINO DOPO LA ROTTURA.96
FIGURA 6.18 - SEQUENZA FOTOGRAFICA PER UNA PROVA DI TRAZIONE. .................................................... 98
FIGURA 6.19 - MESH ADOTTATA NELLE PROVE CON ARMCO-IRON........................................................... 99
FIGURA 6.20 – CONFRONTO TRA IL PROFILO TEORICO E QUELLI NUMERICI DELL’ONDA DI TRAZIONE. ....... 99
FIGURA 6.21 - CONFRONTO TRA DATI NUMERICI E SPERIMENTALI DELLA STRIZIONE. .............................. 100
FIGURA 6.22 - DEFORMATA FINALE DEL PROVINO DI ARMCO IRON UN ISTANTE DOPO LA ROTTURA......... 101
FIGURA 6.23 - CONFRONTO TRA VALORI NUMERICI E SPERIMENTALI DEGLI STRAIN RATES. ..................... 101
FIGURA 6.24 - VARIAZIONE DELLA TEMPERATURA NEL PUNTO PIÙ SOLLECITATO DEL PROVINO. ............. 103
FIGURA 6.25 - DISTRIBUZIONI DELLA TEMPERATURA DOPO LA ROTTURA DEL PROVINO........................... 103
FIGURA 7.1 – A) DIAGRAMMA LAGRANGIANO CARATTERISTICO DI UN IMPATTO PLANARE SIMMETRICO; B)
TIPICO PROFILO DI VELOCITÀ RILEVATO IN UN FLYER PLATE IMPACT TEST.................................... 107
FIGURA 7.2 - ONDA DI STRESS GENERATA IN UN FLYER PLATE IMPACT TEST A VELOCITÀ MODERATA. ... 107
FIGURA 7.3 - PROFILO DI UN'ONDA D'URTO............................................................................................... 108
FIGURA 7.4 – EFFETTO DELLO SMORZAMENTO NUMERICO SUI RISULTATI DELLE SIMULAZIONI NUMERICHE.
........................................................................................................................................................ 110
FIGURA 7.5 - CONFRONTO TRA IL PROFILO DI VELOCITÀ OTTENUTO NUMERICAMENTE ED I DATI
SPERIMENTALI................................................................................................................................. 110
FIGURA 7.6 – EVOLUZIONE NEL TEMPO DELLA DISTRIBUZIONE DELLA TRIASSIALITÀ DELLO STATO DI
SFORZO, LUNGO LO SPESSORE DEL PROVINO.
.................................................................................. 111
FIGURA 7.7 - EVOLUZIONE NEL TEMPO DEL DANNO LUNGO O SPESSORE DEL DISCO BERSAGLIO............... 111
FIGURA 7.8 – DISTRIBUZIONE DELLA POROSITÀ NEL RAME PER IMPATTI A DIVERSE PRESSIONI, [2].......... 112
FIGURA 7.9 - CONFRONTO TRA GLI SPALL SIGNALS CALCOLATO E MISURATO PER IL RAME OFHC. ......... 113
FIGURA 7.10 – DIFFERENTI MECCANISMI DI COALESCENZA DEI MICROVUOTI NELLA ROTTURA DUTTILE.. 114
FIGURA 7.11 - EFFETTO DEL COEFFICIENTE DI FORMA
α
SULLA RISPOSTA DEL SISTEMA DI MOLLE NON
LINEARE. ......................................................................................................................................... 116
FIGURA 7.12 - PROFILO DI VELOCITÀ CALCOLATO CON L’IMPIEGO DEL SISTEMA DI MOLLE NON LINEARE A
CONFRONTO CON I RISULTATI SPERIMENTALI. ................................................................................. 116
FIGURA 7.13 - SCHEMA RIASSUNTIVO DELLE CONFIGURAZIONI GEOMETRICHE ESAMINATE. .................... 118
12
FIGURA 7.14 - PROFILI DI VELOCITÀ CALCOLATI NUMERICAMENTE PER LE DIVERSE CONFIGURAZIONI
GEOMETRICHE E CON VELOCITÀ D’IMPATTO DI 185M/S: A) D/H=16; B) D/H=8; C) D/H=4; D) D/H=2.
........................................................................................................................................................ 118
FIGURA 7.15 - A) DEFORMATA E STATO DEL MATERIALE OTTENUTI CON AUTODYN PER VELOCITÀ
D’IMPATTO DI 185M/S E D/H=16; B) DEFORMATA E MAPPA DEL DANNO OTTENUTI CON MSC/MARC
PER VELOCITÀ D’IMPATTO DI 185M/S E D/H=16. ............................................................................. 119
FIGURA 7.16 - IMPULSO DI COMPRESSIONE IN DUE DIFFERENTI POSIZIONI LUNGO IL RAGGIO DEL DISCO
BERSAGLIO: SULL’ASSE DI SIMMETRIA IN BLU E IN CORRISPONDENZA DEL BORDO LIBERO DEL
PROIETTILE IN NERO. ....................................................................................................................... 120
FIGURA 7.17 - SCHEMA GEOMETRICO DELLA LOCALIZZAZIONE DELL’INNESCO DEL PROCESSO DI SPALL. 121
FIGURA 7.18 - A) DEFORMATA E STATO DEL MATERIALE OTTENUTI CON AUTODYN PER VELOCITÀ
D’IMPATTO DI 185M/S E D/H=8; B) DEFORMATA E MAPPA DEL DANNO OTTENUTI CON MSC/MARC
PER VELOCITÀ D’IMPATTO DI 185M/S E D/H=8. ............................................................................... 122
FIGURA 7.19 - A) DEFORMATA E STATO DEL MATERIALE OTTENUTI CON AUTODYN PER VELOCITÀ
D’IMPATTO DI 185M/S E D/H=4; B) DEFORMATA E MAPPA DEL DANNO OTTENUTI CON MSC/MARC
PER VELOCITÀ D’IMPATTO DI 185M/S E D/H=4. ............................................................................... 122
FIGURA 7.20 - A) DEFORMATA E STATO DEL MATERIALE OTTENUTI CON AUTODYN PER VELOCITÀ
D’IMPATTO DI 185M/S E D/H=2; B) DEFORMATA E MAPPA DEL DANNO OTTENUTI CON MSC/MARC
PER VELOCITÀ D’IMPATTO DI 185M/S E D/H=2. ............................................................................... 122
FIGURA 7.21 - A) DEFORMATA E STATO DEL MATERIALE OTTENUTI CON AUTODYN PER VELOCITÀ
D’IMPATTO DI 185M/S E D/H=1; B) DEFORMATA E MAPPA DEL DANNO OTTENUTI CON MSC/MARC
PER VELOCITÀ D’IMPATTO DI 185M/S E D/H=1. ............................................................................... 123
FIGURA 7.22 – SCHEMA DELLA CONFIGURAZIONE DEL RE-SHOCK EXPERIMENT. ...................................... 124
FIGURA 7.23 - PROFILO DI VELOCITÀ MISURATO IN UN RE-SHOCK EXPERIMENT, [3]................................. 124
FIGURA 7.24 - DISTRIBUZIONE DELLA DEFORMAZIONE PLASTICA LUNGO LO SPESSORE DEL DISCO
BERSAGLIO, A SEGUITO DELL’ONDA DI COMPRESSIONE, IN UN FLYER PLATE IMPACT TEST. ........... 126
FIGURA 7.25 - PROFILI DI SFORZO, CALCOLATI NUMERICAMENTE, A DIVERSE POSIZIONI LUNGO LO SPESSORE
DEL DISCO BERSAGLIO, IN UN FLYER PLATE IMPACT TEST STANDARD............................................ 126
FIGURA 7.26 - PROFILI DI SFORZO, CALCOLATI NUMERICAMENTE, A DIVERSE POSIZIONI LUNGO LO SPESSORE
DEL DISCO BERSAGLIO, IN UN RE-SHOCK EXPERIMENT. ................................................................... 127
FIGURA 7.27 - RAPPRESENTAZIONE DELLO STATO DI SFORZO DEL PUNTO MATERIALE CHE A SUBITO UNO
SHOCK, A DIFFERENTI POSIZIONI LUNGO LO SPESSORE DEL DISCO BERSAGLIO, NEL PIANO DEI
DEVIATORI
π . ................................................................................................................................ 128
FIGURA 7.28 - CONFRONTO TRA IL PROFILO DI VELOCITÀ SPERIMENTALE, [3], E QUELLO CALCOLATO
NUMERICAMENTE CON MSC.MARC. ............................................................................................... 129
FIGURA 7.29 - CONFRONTO TRA IL PROFILO DI VELOCITÀ SPERIMENTALE, [3], E QUELLO CALCOLATO
NUMERICAMENTE CON AUTODYN. .................................................................................................. 130
FIGURA 7.30 - RITARDO DEL GRADINO ANOMALO RISPETTO ALLA CORRISPONDENTE ONDA DI RILASCIO IN
UN FLYER PLATE IMPACT TEST STANDARD, [4]. ............................................................................. 131
13
Indice delle Tabelle
TABELLA 5.1 - PROPIETÀ MECCANICHE DEI MATERIALI INVESTIGATI E RELATIVI PARAMETRI PER IL
MODELLO DI JOHNSON E COOK, [1]. .................................................................................................. 63
TABELLA 5.2 - CONFRONTO TRA GLI ACCORCIAMENTI CALCOLATI, CON DIVERSI MODELLI DI RESISTENZA, E
I DATI SPERIMENTALI.
....................................................................................................................... 65
TABELLA 5.3 - CONFRONTO TRA I DIAMETRI CALCOLATI DELLE SUPERFICI D'IMPATTO E QUELLI MISURATI.
.......................................................................................................................................................... 69
TABELLA 5.4 - PARAMETRI DI DANNO PER IL RAME OFHC......................................................................... 73
14
1 Introduzione
La complessità dei meccanismi correlati alla dinamica dell’impatto limita l’impiego dei
modelli analitici al solo scopo di sviluppare una percezione immediata dei fenomeni
investigati. Nei casi reali, infatti, anche per le configurazioni più semplici, è molto facile
violare le ipotesi su cui è basata l’analisi teorica, la quale molto difficilmente può essere
utilizzata per effettuare previsioni.
Per la soluzione di problemi d’impatto è indispensabile ricorre agli strumenti della
simulazione numerica. Il suo impiego, ad oggi, per quanto riguarda la gestione dei
transitori e la risoluzione del moto di propagazione delle onde, è largamente diffuso e
saldamente consolidato. Ciò, comunque, non deve far pensare ad una sorta d’infallibilità
dei codici numerici, con il pericolo di compiere gravi errori di valutazione. La bontà dei
risultati delle simulazioni è direttamente legata alla qualità del modello numerico
realizzato. Qualità che, a sua volta, non può prescindere da una completa e corretta
conoscenza dei meccanismi che si desidera modellare.
In questo lavoro, gli strumenti della simulazione numerica sono stati adoperati per
analizzare tre configurazioni sperimentali classiche per la caratterizzazione della
risposta meccanica dei materiali in regime dinamico: il Taylor Test, la Hopkinson Bar e
il Flyer Plate Impact Test. L’obiettivo è stato quello di andare a studiare, nelle tre
configurazioni, già largamente investigate nel corso degli anni, i punti di maggiore
criticità: si sono indagati i limiti di modelli proposti in letteratura e le condizioni in cui
gli stessi falliscono, cercando di comprenderne il motivo e, ove possibile, di superarli.
In quest’ottica, i codici numerici, oltre che come strumenti di previsione, sono stati
impiegati come veri e propri strumenti d’investigazione. Questo modo di operare ha
permesso di risolvere alcune apparenti incongruenze derivanti dai risultati sperimentali
e di fornire nuove interpretazioni di fenomeni dinamici che non erano stati pienamente
compresi.
Inoltre, si è fornita una dimostrazione dell’enorme potenzialità degli strumenti numerici
in un campo, quale quello della dinamica dell’impatto, in cui anche le misure
sperimentali sono assai difficoltose e richiedono sempre una notevole capacità
interpretativa.
15
2 Onde di sollecitazione nei solidi
2.1 Introduzione
Una perturbazione, esercitata su una qualche grandezza fisica in una regione limitata
dello spazio, si propaga nello spazio circostante con modalità che dipendono di norma
dal tipo di perturbazione e dalle caratteristiche del mezzo che riempie lo spazio. E’ bene
puntualizzare che la propagazione, la quale avviene tramite un’onda, è del disturbo e
non della grandezza in esame e che tale propagazione ondosa comporta uno scambio di
energia. A causa di questi scambi, parte dell’energia meccanica viene convertita in
calore attraverso diversi meccanismi indicati, in genere, come attriti interni. Questi
introducono una tale complessità nei modelli matematici che descrivono il moto
ondoso, da renderli intrattabili. E’ questo il motivo per cui tali effetti vengono trascurati
nella maggior parte delle trattazioni senza tuttavia inficiare la loro validità, [1].
Nella teoria elasto-plastica dei corpi solidi, è stata diffusamente trattata, Kolsky [2],
Johnson [3], Achenbach [4], la propagazione di due tipi di onde:
1 le onde longitudinali in cui il moto delle particelle del mezzo si sviluppa
parallelamente alla direzione di propagazione;
2 le onde trasversali in cui il moto delle particelle del mezzo si sviluppa in
direzione perpendicolare alla direzione di propagazione.
Bisogna però ricordare che, in generale, quando un corpo è soggetto a carichi impulsivi,
si sviluppano in esso diverse tipologie di onde, quali quelle di torsione o quelle di
flessione; se il disturbo si propaga lungo una superficie del corpo, questa avviene per
mezzo di onde superficiali quali:
1 le onde di Rayleigh la cui intensità decade in modo esponenziale con la distanza
dalla superficie;
2 le onde di Love, onde di taglio che si formano in materiali composti da strati con
diverse caratteristiche fisiche.
Nel presente capitolo sono esposti alcuni cenni di teoria della propagazione delle onde
nei solidi, in riferimento al caso di onde longitudinali. Sono trattati la derivazione
16
dell’equazione delle onde elasto-plastiche e l’applicazione ai due casi di riferimento di
tensione ovvero deformazione uniassiale. La teoria delle onde nei solidi è trattata con
gli strumenti della “Rate Independent Theory”, per cui la risposta costitutiva del
materiale è descritta da una curva sforzo deformazione. Questo semplifica enormemente
l’analisi senza, per altro, compromettere la trattazione, in quanto il problema di
propagazione delle onde è piuttosto insensibile alla forma dell’equazione costitutiva.
2.2 Equazione delle onde
Per rendere la formulazione dell’equazione delle onde concettualmente chiara e
maggiormente intuitiva è necessario, almeno inizialmente, limitare la discussione al
caso di propagazione ondosa monodimensionale.
Si consideri il problema di un disturbo che viaggia nella direzione x rispetto ad un
sistema di riferimento fisso e si esaminino, come mostrato in figura 2.1, le forze che
agiscono su un elemento di massa dm = ρ * dx * A . Si ipotizzi di limitare l’analisi al
caso di piccole deformazioni e piccoli spostamenti.
Figura 2.1 - Forze agenti sull’elemento di massa, [5].
Sia σ =
F
la tensione, definita positiva se in trazione, per il teorema della quantità di
A
moto si può scrivere:
∂σ
∂v
=ρ
∂x
∂t
(2.1)
avendo indicato con v la velocità della particella, ovvero la derivata rispetto al tempo
17
dello spostamento nella direzione x: v =
∂u
.
∂t
Considerando che la deformazione è definita come: ε =
∂u
ne deriva che:
∂x
∂ε ∂ v
=
∂t ∂ x
(2.2)
Assumendo la tensione come funzione biunivoca della deformazione: σ = σ (ε ) si
ricava l’equazione dell’onda per moto monodimensionale:
∂ 2 u 1 ∂ 2u
=
∂ x2 c2 ∂ t 2
(2.3)
in cui c è la velocità di propagazione del fronte d’onda:
1 dσ
ρ dε
c(ε ) =
(2.4)
E’ immediato riconoscere questa equazione come il caso particolare, di propagazione
monodimensionale, della più generale equazione delle onde in forma indiciale:
∂U
1 ∂ 2ψ
= 2
∂ xi∂ xi c ∂ t 2
(2.5)
2.2.1 Tensione generata dall’impatto
Per derivare l’intensità dello sforzo generato in un impatto, si può far riferimento
all’evento, rappresentato in Figura 2.2, di un muro rigido che, al tempo t = 0 , impatta, a
velocità v = v0 , una barra o un disco in stato di quiete. Nell’intervallo di tempo dt la
barra si deformerà fino al piano B che dista v0 ⋅ dt dalla propria estremità originaria. Il
disturbo, che porta la velocità delle particelle a v0 , viaggerà, nello stesso intervallo di
tempo, fino al piano A, per una distanza pari a c ⋅ dt , in cui c indica la velocità
dell’onda. Se si indica con σ lo sforzo di compressione che si genera tra l’impattatore e
la barra di sezione A0 , l’impulso generato dell’intervallo di tempo dt è pari a σ A0 dt .
La quantità di moto della barra, inizialmente ferma, è pari a ρ A0 cdt ⋅ v0 , al prodotto,
cioè, della velocità per la massa delle particelle comprese dall’estremità iniziale della
18
barra ed il fronte d’onda A. Uguagliando l’impulso alla variazione della quantità di
moto, si ottiene:
σ = ρ cv0
(2.6)
Se lo stato iniziale di sforzo e velocità è non nullo, le quantità σ e v0 devono essere
sostituite dalle loro corrispettive variazioni ∆σ e ∆v , che portano all’espressione più
generale:
∆σ = ρ c∆v
(2.7)
v 0 dt
B
v0
A
v = v0
v=0
cdt
Figura 2.2 – Schema di un corpo rigido che impatta un cilindro a velocità
v0 , [5].
2.2.2 Riflessione di onde elastiche alle interfacce
Se, come schematicamente illustrato in Figura 2.3, un’onda longitudinale raggiunge una
superficie libera con un generico angolo d’incidenza α, dalla sua riflessione saranno
generate due onde distinte.
La prima, anch’essa longitudinale, sarà riflessa con un angolo pari a quello d’incidenza,
la seconda, di tipo distorsionale, sarà riflessa con un angolo più piccolo tale che:
sin β 2 cD
=
sin α1 cL
(2.8)
In cui cL e cD indicano rispettivamente le velocità di propagazione dell’onda
longitudinale e distorsionale.
19
onda incidente
α1
α2
β2
distorsionale
onde riflesse
longitudinale
Figura 2.3 – Riflessione di un’onda longitudinale su una superficie libera.
u, v
+σ
cL
cL
ghost
-σ
u, v
σnet
Figura 2.4 – Riflessione di un’onda longitudinale su una superficie libera ortogonale alla direzione
di propagazione dell’onda, [1].
Allo stesso modo, se un’onda distorsionale raggiungesse una superficie libera si
genererebbero due onde, quella trasversale avrebbe un angolo di riflessione pari a quello
d’incidenza, quella longitudinale sarebbe riflessa con un angolo che rispetterebbe la
relazione precedente (2.8).
In Figura 2.4, è schematicamente illustrato il caso particolare di un’onda longitudinale
che impatta normalmente una superficie libera. Poiché lo sforzo di tensione
perpendicolare alla superficie deve essere nullo, l’impulso riflesso dovrà essere di segno
inverso a quello incidente. In altre parole
20
un impulso di compressione sarà riflesso
come un impulso di trazione e viceversa. Sia u I = f ( x − ct ) lo spostamento, lungo
l’asse positivo delle ascisse, dovuto all’impulso incidente. La riflessione sulla superficie
libera genera un’onda, che si muove lungo l’asse negativo delle ascisse, che porta ad
uno spostamento u I = g ( x + ct ) . Alla superficie libera, per x = l , poiché lo sforzo netto
deve essere nullo, si ha che:
σ NET = σ I + σ R = 0
(2.9)
Se si esprime lo sforzo in campo elastico come σ = Eε = E ( ∂u ∂x ) , si ottiene:
σ NET = E ⎡⎣ f ′ ( l − ct ) + g ′ ( l + ct ) ⎤⎦ = 0
o,
(2.10)
f ′ ( l − ct ) = − g ′ ( l + ct )
in cui l’apice indica la derivazione rispetto a x . Dall’equazione (2.10) si evince che gli
impulsi incidente e riflesso hanno la stessa forma, ma segno opposto. Anche la velocità
della particella sulla superficie libera può essere ottenuta per sovrapposizione:
∂uI ∂uR
+
∂t
∂t
(2.11)
vNET = c ( − f ′ + g ′ ) = 2cg ′
(2.12)
vNET = vI + vR =
che, sulla superficie libera, x = l , porta a:
in cui, in questo caso, gli apici esprimono la derivazione rispetto al tempo. L’equazione.
(2.12) attesta che nella regione in cui gli impulsi incidente e riflesso si sovrappongono,
la velocità delle particelle e, quindi, anche lo spostamento, sono il doppio di quelli
sviluppati dagli impulsi singoli.
La tecnica, utilizzata in Figura 2.4 per visualizzare il comportamento degli impulsi di
sforzo all’interfaccia, sfrutta la linearità dell’equazione dell’onda elastica per ottenere la
soluzione come sovrapposizione di due impulsi: il primo impulso, in rosso, è costituito
dall’onda incidente “reale”, il secondo è un immaginario impulso fantasma (ghost), di
medesima forma, ma di segno opposto, che si trova inizialmente all’esterno del
materiale e viaggia nel verso contrario. Quando raggiungono l’interfaccia, il primo
21
impulso esce dal materiale, mentre il secondo, entrando, diviene progressivamente
reale. All’interno del materiale, laddove i due impulsi si sovrappongono, l’impulso netto
è nullo. Quando l’impulso incidente esce completamente dal materiale, quello che
originariamente era stato indicato come impulso fantasma dà luogo all’impulso, di
forma quadra, disegno contrario a quello d cui è stato generato.
Su una superficie fissa, invece, come illustrato in Figura 2.5, la velocità e lo
spostamento devono essere nulli. Per cui, seguendo lo stesso procedimento, si può
scrivere:
vNET = −cf ′ ( l − ct ) + cg ′ ( l + ct ) = 0
o,
(2.13)
f ′ ( l − ct ) = g ′ ( l + ct )
e per lo sforzo:
⎛ ∂u
∂u ⎞
σ NET = E ⎜ I + R ⎟ = E ⎡⎣ f ′ ( l − ct ) + g ′ ( l + ct ) ⎤⎦ = 2 Ef ′ ( l − ct )
∂x ⎠
⎝ ∂x
(2.14)
Per cui, lo sforzo, sulla superficie vincolata è il doppio, mentre la velocità e lo
spostamento sono nulli.
u, v
u, v
+σ
cL
cL
ghost
+σ
u,v = 0
Figura 2.5 – Riflessione di un’onda longitudinale su una superficie fissa ortogonale alla direzione di
propagazione dell’onda, [1].
22
2.2.3 Riflessione e trasmissione di onde elastiche in una discontinuità
meccanica
Si consideri una bara con una discontinuità meccanica dovuta ad una differenza di
materiale o ad una variazione di sezione, Figura 2.6. Sia σ I un impulso elastico di
compressione che viaggia nella barra verso destra. Alla discontinuità questo sarà in
parte riflesso, σ R , e in parte trasmesso σ T , in modo che siano verificate le seguenti
condizioni:
‰ all’interfaccia, la forza nelle due barre deve essere la medesima,
A1 (σ I + σ R ) = A2 (σ T )
(2.15)
in cui A1 e A2 sono le rispettive sezioni delle barre all’interfaccia;
‰ le velocità delle particelle all’interfaccia devono essere continue,
vI + vR = vT
(2.16)
σI
σ
σ
− R = T
ρ1c1 ρ1c1 ρ 2 c2
(2.17)
Ricordando che σ = ρ cv0 , si ottiene:
e risolvendo in funzione di σ I :
A1, ρ1, c1
σT =
2 A1 ρ 2 c2
σI
A1 ρ1c1 + A2 ρ 2 c2
(2.18)
σR =
A2 ρ 2 c2 − A1 ρ1c1
σI
A1 ρ1c1 + A2 ρ 2 c2
(2.19)
σI
σR
σT
A2 , ρ2 , c2
Figura 2.6 – Riflessione e trasmissione di un’onda ad una discontinuità meccanica.
Le equazioni (2.18) e (2.19) permettono di fare alcune considerazioni che possono
23
rivelarsi utili nell’analisi delle configurazioni sperimentali che saranno esaminate nel
prosieguo della presente trattazione. Se i materiali che costituiscono le due barre sono
identici si ha che ρ1 = ρ 2 e c1 = c2 , per cui:
σT =
2 A1
σI
A1 + A2
(2.20)
σR =
A2 − A1
σI
A1 + A2
(2.21)
σ T e σ R avranno lo stesso segno se A2 > A1 . Se, invece A2 < A1 , σ T e σ R avranno
segno opposto. Se A2 A1 → 0 , si tende alla condizione di superficie libera e quindi
σ R → −σ I . Se, al contrario, A2 A1 → ∞ , si tende alla condizione di superficie fissa, per
cui σ R → σ I . e σ T → 0 .
Non si verificano riflessioni dell’onda incidente, σ R = 0 , quando le impedenze
meccaniche delle due barre sono tra di loro uguali, A1 ρ1c1 = A2 ρ 2 c2 , da cui di ricava:
σT = σ I
E2 ρ 2
E1 ρ1
(2.22)
Nell’equazione (2.18), il coefficiente di σ I , non può mai essere negativo; questo
significa che un impulso incidente di tensione sarà sempre trasmesso come un impulso
di tensione e che un impulso incidente di compressione sarà sempre trasmesso come un
impulso di compressione.
Nell’equazione (2.19), il coefficiente di σ I , può essere positivo o negativo a seconda
che si abbia A1 ρ1c1 < A2 ρ 2 c2 o A1 ρ1c1 > A2 ρ 2 c2 rispettivamente. Se il coefficiente è
negativo, A1 ρ1c1 > A2 ρ 2 c2 , un impulso incidente di compressione sarà riflesso come un
impulso di trazione e vice versa. Se il coefficiente è positivo, A1 ρ1c1 < A2 ρ 2 c2 , gli
impulsi incidente e riflesso avranno lo stesso segno.
2.2.4 Tensione uniassiale
Adottare configurazioni che garantiscano la possibilità di effettuare alcune
semplificazioni permette di rendere matematicamente trattabili i modelli che descrivono
24
il moto delle onde. Se si studia la propagazione di un impulso di tensione in una barra
sottile, che abbia cioè una lunghezza pari o maggiore di dieci volte il suo diametro, è
possibile trascurare gli effetti dell’inerzia trasversale. Si può dunque assumere lo stato
di tensione monoassiale.
Nel caso in cui lo sforzo sia inferiore alla tensione di snervamento del materiale, lo
stesso si comporterà elasticamente e, per quanto ricavato precedentemente, nel corpo si
propagherà una perturbazione longitudinale il cui moto è descritto dall’equazione:
∂ 2 u 1 ∂ 2u
=
∂ x2 c2 ∂ t 2
(2.23)
e la cui velocità è pari a:
E
c=
(2.24)
ρ
dove si è indicato con E il modulo di Young.
Si consideri ora il caso di un materiale, il cui modello costitutivo sia descritto, secondo
le ipotesi della “Rate Independent Theory”, da una legge tensione-deformazione
bilineare, come in Figura 2.7 a, sottoposto ad un impulso di valore superiore alla sua
tensione di snervamento. Nella barra considerata si propagheranno due distinti fronti
d’onda, come illustrato nella Figura 2.7 b. Ogni fronte d’onda avrà una propria velocità
di propagazione che dipenderà dai rispettivi moduli di elasticità E ed E1.
σ
σ
E1
σy
σy
ε
ce
σy
E
⋅t
ρ
E
(a)
σ
E1
⋅t
ρ
(b)
x
cp
ce
cp
(c)
x
Figura 2.7 – Curva tensione deformazione e profili d’onda per un materiale bilineare.
Se l’impulso è di breve durata, nel solido si genereranno le onde di rilascio elastica e
plastica che, viaggiando alle velocità che le competono, porteranno alla formazione del
profilo d’onda illustrato in Figura 2.7 c.
25
Poiché la velocità di propagazione delle onde elastiche può essere anche dieci volte
maggiore della velocità delle onde plastiche, l’onda di rilascio elastica potrebbe
raggiungere l’onda plastica e scaricarla. A questo punto un’altra onda, generata dalla
riflessione sulla discontinuità rappresentata dall’onda plastica, si dirigerebbe verso la
superficie libera della barra. Si innescherebbe così un meccanismo di continue
riflessioni delle onde, tra la superficie libera e la posizione del fronte dell’onda plastica,
che porterebbe, se fossero disponibili tempi sufficientemente lunghi, al progressivo
scarico dell’onda incidente iniziale.
ε
σ
ε1
σy
ε1
εe
ε
c1
(a)
c0
ξ
(b)
Figura 2.8 – Curva tensione deformazione e profilo d’onda per un materiale elasto-plastico,
secondo la “rate independent theory”.
Nel caso in cui la curva sforzo-deformazioni possa essere rappresentata, ancora sotto le
ipotesi di “Rate Independent Theory”, dalla curva, riportata in Figura 2.8 a, con
variazione continua della pendenza della parte plastica, ne risulterà il profilo di velocità
di Figura 2.8 b, in cui si è definito ξ = x . Questo è il tipico profilo d’onda che si
t
sviluppa in un materiale elasto-plastico, quale ad esempio un metallo, in condizione di
sforzo uniassiale. Dall’equazione (2.4), infatti, si deduce che ogni livello di tensione o
deformazione propaga con una propria velocità caratteristica, che è funzione della
tangente locale alla curva sforzo-deformazione. Poiché la curva in Figura 2.8 a, assunta
come rappresentativa del comportamento del materiale in condizione di sollecitazione
uniassiale, presenta la concavità verso l’asse delle ascisse, disturbi di tensione o
deformazione più elevati sono caratterizzati da una più bassa velocità di propagazione.
L’onda, quindi, è costituita da un precursore elastico seguito da un più lento fronte
d’onda plastico disperso, a sua volta seguito da una regione a deformazione plastica
( ε1 )
costante.
26
2.2.5 Deformazione uniassiale
Un’altra importante configurazione è quella che prevede che la deformazione possa
avvenire in una sola direzione. Uno stato di deformazione monoassiale è definito come:
ε1 ≠ 0
ε 2 = ε 3 = γ 12 = γ 13 = γ 23 = 0
(2.25)
Nel derivare le equazioni per stato di deformazione monoassiale, si assume che la
deformazione totale si possa scomporre in una parte elastica ed una plastica:
ε1 = ε1e + ε1p
ε 2 = ε 2e + ε 2p
(2.26)
ε 3 = ε 3e + ε 3p
Per la seconda delle (2.25), si ha che:
ε 2p = −ε 2e
(2.27)
ε = −ε
p
3
e
3
Per l’incompressibilità del flusso plastico si può scrivere che:
ε1p + ε 2p + ε 3p = 0
(2.28)
che, sfruttando la simmetria, ε 2p = ε 3p ,porta a:
ε1p = −ε 2p − ε 3p = −2ε 2p
(2.29)
Utilizzando l’equazione (2.27) si ottiene:
ε1p = 2ε 2e
(2.30)
in modo da poter scrivere la deformazione totale in termini di sola deformazione
elastica:
ε1p = ε1e + ε1p = ε1e + 2ε 2e
(2.31)
Le deformazioni elastiche possono essere espresse, in termini di sforzi, dalle seguenti
relazioni:
27
ε1e =
ε 2e =
ε 3e =
σ1 ν
E
σ2
E
−
−
E
ν
E
σ3 ν
E
−
E
(σ 2 + σ 3 ) =
σ1
−
E
2ν
σ2
E
(σ 1 + σ 3 ) =
(1 −ν ) σ
2
−
(σ 1 + σ 2 ) =
(1 −ν ) σ
3
−
E
E
ν
E
ν
E
σ1
(2.32)
σ1
avendo posto σ 2 = σ 3 . La combinazione delle equazioni (2.32) e (2.31) permette di
ottenere:
ε1e =
σ 1 (1 − 2ν ) 2σ 2 (1 − 2ν )
E
+
E
(2.33)
Imponendo come criterio di snervamento quello di von Mises o quello di Tresca, cioè:
σ 1 − σ 2 = Y0
(2.34)
in cui Yo indica la tensione di snervamento, si ottiene:
σ1 =
E
2
2
ε1 + Y0 = K ε1 + Y0
3 (1 − 2ν )
3
3
(2.35)
in cui il “bulk modulus”, K , è definito come:
K=
E
3 (1 − 2ν )
(2.36)
Nel caso particolare di deformazione elastica unidimensionale:
ε1 = ε1e
ε 2 = ε 2e = ε 3 = ε 3e = 0
(2.37)
ε1p = ε 2p = ε 3p = 0
per cui,
ε 2e = 0 =
1 −ν
ν
σ 2 − σ1
E
E
28
(2.38)
ovvero,
σ2 =
ν
(1 −ν )
σ1
che porta alla scrittura;
ε1 =
σ1
E
−
2ν 2σ 1
E (1 −ν )
o,
σ1 =
1 −ν
Eε
(1 − 2ν )(1 +ν ) 1
(2.39)
L’ equazione (2.39) dimostra che, in caso di stato di deformazione unidimensionale, la
pendenza del tratto elastico della curva sforzo-deformazione del materiale, è, rispetto al
caso di stato di sforzo unidimensionale, più elevato di un coefficiente pari a
1 −ν
.
(1 − 2ν )(1 +ν )
Questo è chiaramente illustrato nella Figura 2.10, in cui nella parte di sinistra sono
schematicamente illustrate le curve sforzo-deformazione, per uno stato di sforzo
unidimensionale, dei due materiali elastico perfettamente plastico e elastico con
incrudimento lineare, mentre nella parte destra, sono riportate le corrispettive curve che
si ottengono, per i medesimi materiali, in caso di deformazione unidimensionale.
Figura 2.9 – Curva sforzo-deformazione per i materiali elastico perfettamente plastico ed elastico
con incrudimento lineare, in caso di stato di sforzo uniassiale ovvero di deformazione uniassiale,
[1].
29
Un altro risultato interessante è l’innalzamento del valore della σ 1 per il quale si ha il
superamento del limite elastico, dal valore dello snervamento del materiale, Y0 , per il
caso di sforzo uniassiale, allo “Hugoniot Elastic Limit”, σ HEL , per il caso di
deformazione uniassiale. Per quanto riguarda la parte plastica, l’equazione (2.35)
dimostra che lo stress, indipendentemente dall’incrudimento, continua a crescere con la
deformazione, in modo proporzionale al “bulk modulus”, e che lo scostamento dalla
parte idrostatica della curva è pari a un valore costante 2 Y0 . La curva indicata in
3
Figura 2.9 come “Hydrostat” rappresenta, quindi, il comportamento del medesimo
materiale, ma privo di capacità di resistenza a taglio, soggetto ad uno stato di
deformazione uniassiale. Per valori estremamente elevati della pressione, lo
scostamento tra le due curve diviene trascurabile e il materiale, senza compiere errori
significativi, può essere trattato come un fluido e rappresentato dalla sola parte
idrostatica.
Se al materiale elastico perfettamente plastico rappresentato in Figura 2.9 si applica uno
sforzo che supera il limite elastico di Hugoniot, si sviluppano le due onde, elastica e
plastica, che secondo le equazioni (2.4), (2.35) e (2.39) propagano, rispettivamente, con
velocità:
ce =
E (1 −ν )
ρ0 (1 − 2ν )(1 +ν )
(2.40)
e
cp =
K
ρ0
(2.41)
Contrariamente a quanto avviene in stato di sforzo uniassiale, le velocità di
propagazione delle due onde non sono significativamente differenti. Per un tipico
acciaio legato, ad esempio, l’onda elastica è più veloce di quella plastica di circa il 25%,
mentre, nel caso di sforzo uniassiale, si arriva ad un fattore pari a 10.
Per valori della pressione estremamente elevati, lo sforzo idrostatico perde il rapporto di
proporzionalità lineare con la compressione volumetrica, per salire più rapidamente,
come schematicamente illustrato in Figura 2.10. La curva, contrariamente a quanto si
30
verifica in uno stato di sforzo uniassiale, presenta la concavità rivolta verso l’asse delle
ordinate, con fortissime implicazioni sullo sviluppo e la propagazione delle onde di
sollecitazione nel materiale.
Figura 2.10 – Curva sforzo-deformazione in stato di deformazione uniassiale, per valori della
pressione estremamente elevati, [5].
Il punto A corrisponde al limite elastico di Hugoniot e, quindi, il precursore elastico
viaggerà ancora alla velocità governata dalla pendenza OA. Se il disturbo è tale da
portare lo stato del materiale oltre il valore dell’HEL, si ricade nel tratto di curva con
concavità verso l’alto, per cui, secondo la relazione (2.4), agli sforzi plastici più severi
compete una velocità di propagazione più elevata rispetto a quella degli sforzi plastici
più deboli. Questo comporta che se, per esempio, la sollecitazione porta lo stato del
materiale fino al punto B di Figura 2.10, nello stesso si genereranno un precursore
elastico, che viaggerà alla velocità definita dall’equazione (2.40), e un’onda d’urto
plastica, che viaggerà alla velocità dettata dalla pendenza del tratto AB. Se il livello di
sforzo raggiunge il punto C, che si trova sul prolungamento del tratto OA, il precursore
elastico e l’onda d’urto plastica viaggeranno alla stessa velocità. Infine, per pressioni
ancora più elevate, fino al punto D, si svilupperà una singola onda d’urto, la cui
velocità, maggiore di quella che competerebbe al precursore elastico, è determinata
dalla pendenza del tratto OD.
Un’implicazione fondamentale che deriva da quanto detto è che la curva di Figura 2.10,
generalmente indicata col termine Hugoniot, rappresenta il luogo dei punti degli stati di
31
equilibrio, ma, diversamente da quanto accade per la curva sforzo deformazione in stadi
sollecitazione uniassiale, non viene percorsa durante il processo di carico.
32
Bibliografia
[1]
Zukas, J. A, Nicholas, T., Swift, H. F., Greszczuc, L. B. and Curran, D. R.,
Impact Dynamic, John Wiley & Sons, New York, 1992.
[2]
Kolsky, H., Stress Waves in Solids, Dover, New York, 1963.
[3]
Johnson, W., Impact Strength of Materials, Crane, Russak, New York, 1972.
[4]
Achenbach, J.D., Wave Propagation in Elastic Solids, American Elsevier, 1975.
[5]
Zukas, J. A., High Velocity Impact Dynamics, John Wiley & Sons, New York,
1990.
33
3 Modellazione costitutiva
3.1 Introduzione
Lo studio dei processi associati ai fenomeni impulsivi, quali quelli conseguenti
l’impatto tra corpi, richiede necessariamente la conoscenza e la descrizione accurata del
comportamento meccanico del materiale in regime dinamico. Tali modelli devono
essere in grado di tenere in conto gli effetti, sulla resistenza del materiale, della
deformazione, della velocità di deformazione, della temperatura, del danneggiamento e,
per impatti iperveloci, della pressione idrostatica. Un approccio consolidato, e
largamente utilizzato, è quello di utilizzare un modello costitutivo composto da più
sottomodelli disaccoppiati, in grado di tenere in conto gli effetti combinati di tutte le
variabili in gioco. Tradizionalmente si descrivono gli effetti della deformazione, della
velocità di deformazione e della temperatura con un modello di resistenza, gli effetti del
danneggiamento con un modello di rottura e l’effetto della compressione volumetrica
con un’equazione di stato. Tale modo di operare semplifica enormemente la descrizione
del comportamento del materiale e, soprattutto, permette la caratterizzazione dello
stesso con un numero limitato di prove meccaniche, relativamente semplici. Anche se a
rigore si dovrebbe ricorrere ad una formulazione costitutiva che incorpori questi effetti
in maniera accoppiata, esistono numerose osservazioni sperimentali che giustificano
tale modo di operare, 0.
Nei paragrafi seguenti sono presentati i modelli di resistenza ed i modelli di
danneggiamento più importanti, con riferimento particolare ai modelli utilizzati per
l’analisi delle configurazioni sperimentali investigate.
3.2 Modelli strain rate sensitive
I metalli, generalmente, mostrano una notevole sensibilità alla velocità di deformazione
e alla temperatura. Nella maggior parte dei casi, la velocità di deformazione ha l’effetto
più rilevante sull’incremento della resistenza del materiale. In Figura 3.1 sono riportati
gli andamenti sforzo-deformazione di taglio, per un alluminio commerciale, in un
intervallo di velocità di deformazione che va da 600 a 2800 s-1, a confronto con la curva
di riferimento quasistatica ottenuta a 2,0 103 s-1.
34
Figura 3.1 – Effetto della velocità di deformazione sull’alluminio commerciale caricato a taglio.
Figura 3.2 - Effetto della temperatura sul titanio α.
La temperature, al contrario, addolcisce il materiale, come mostrato per il titanio-α in
Figura 3.2, in cui sono riportate le curve sforzo deformazione ottenute, a parità di
velocità di deformazione, in un intervallo di temperature che va dai 77 ai 288K. È noto
che la sensibilità del materiale alla velocità di deformazione e alla temperature è legata
alla struttura atomica. In particolare, i metalli con una struttura cubica a corpo centrale
35
(CCC), quali il ferro α, gli acciai ferritici, il niobio, il tantalio, etc., mostrano una forte
variazione del valore della tensione di snervamento con la temperatura, T, e la velocità
di deformazione, ε . Al contrario, i metalli con struttura cubica a facce centrate (CFC),
quali gli acciai austenitici, il nichel, l’alluminio, il rame e l’argento, non mostrano la
stessa sensibilità in modo particolare rispetto alla temperatura. I metalli a struttura
esagonale compatta (EC), infine, quali il titanio e lo zinco, esibiscono un
comportamento intermedio tra quello dei CCC e quello dei CFC.
Da tali osservazioni risulta essere evidente la necessità di superare le ipotesi
semplificative, utilizzate nel capitolo precedente, della “Rate Independent Theory”e di
utilizzare modelli basati sulla “Rate Dependent Theory”. Esiste un elevato numero di
modelli proposti in letteratura sviluppati facendo riferimento a due approcci differenti:
quelli sviluppati su base fisica, come ad esempio l’energia di attivazione o la meccanica
delle dislocazioni, o, alternativamente, gli approcci empirici. Mentre i primi descrivono
in modo più attento l’insieme dei meccanismi intimamente legati all’evoluzione della
microstruttura del materiale, i secondi sono più semplici da utilizzare in virtù d’una
maggiore maneggevolezza.
3.2.1 Modelli di resistenza formulati su basi fisiche
Descrivere la risposta inelastica di tutti i metalli con una legge generalizzata che derivi
da una teoria unificata è estremamente difficile. Nel passato sono stati compiuti diversi
tentativi con l’obiettivo di ricavare una relazione che leghi la tensione di snervamento,
σy, alla velocità di deformazione a alla temperatura:
σy
= f (ε, ε,T )
(3.1)
Un gran numero di equazioni sono state proposte da autori diversi, ampie descrizioni
sono riportate nei testi di Zukas, [2], e Zukas et al., [3]. Anche se molto differenti e, a
volte, tra di loro inconsistenti, una caratteristica fondamentale, riconoscibile in tutti i
modelli, è quella che vede una dipendenza esponenziale della tensione di snervamento
dalla temperatura e un’equivalenza di effetti tra temperatura e velocità di deformazione.
Tali caratteristiche trovano conferma nei risultati sperimentali quali quelli riportati in
Figura 3.3, in cui, per un acciaio basso legato, è riportata la variazione del valore dello
snervamento con la temperatura e la velocità di deformazione. Qui è illustrato
36
chiaramente come, in un diagramma logaritmico, lo snervamento cresce linearmente
con la velocità di deformazione fino a ε pari a 104 s-1.
Figura 3.3 - Variazione del valore dello snervamento di un acciaio basso legato con la velocità di
deformazione e la temperatura.
Tale aspetto è evidente nella relazione di Zener e Hollomon, [4]:
σy
= f (ε ⋅ eQ / RT )
(3.2)
in cui Q è l’energia di attivazione e R la costante universale dei gas.
Il modello di Zerilli-Armstrong, [5], è, invece, basato sul moto delle dislocazioni
termicamente attivato, con particolare attenzione alla differente risposta dei metalli
CCC da quelli CFC. Il modello ha un’ottima capacità di descrivere i risultati
sperimentali ed è espresso dalle seguenti relazioni:
BCC :
σy
= C 1 exp ( −C 3T + C 4T ln ε )
(3.3)
FCC :
σ flow
= C 2 ε1/ 2 exp ( −C 3T + C 4T ln ε )
Nel 1984 Hartley e Duffy, [6], proposero un modello, basato sulla dinamica delle
dislocazioni, che spiegasse il comportamento d’un materiale che manifesti una
sensibilità sia alla temperatura che alla velocità di deformazione. La legge, ricavata
37
facendo riferimento alla teoria dei meccanismi di attivazione termica, ha la forma:
1
1 p
⎧
⎫
⎡
⎤
⎛ γ0 ⎞ q ⎪
⎪
k * T *ln ⎜ ⎟ ⎥
⎪⎪ ⎢
⎝ γ ⎠ ⎥ ⎪⎪
τ = τ µ + (τ 0 − τ µ ) ⎨1 − ⎢
⎬
⎥ ⎪
F0
⎪ ⎢
⎥ ⎪
⎪ ⎢⎣
⎦ ⎪
⎩⎪
⎭
(3.4)
dove T indica la temperatura assoluta, F0 l’energia libera totale da superare, τ0 la
tensione di snervamento allo zero assoluto, τµ la componente atermica della tensione di
snervamento e p e q descrivono la forma degli ostacoli da superare.
3.2.2 Modelli di resistenza fenomenologici
Si hanno relazioni ancora più complesse se si tenta di descrivere la sensibilità del
materiale alla storia delle velocità di deformazione. Nel 1977, Campbell et al., [7],
propose un modello fenomenologico nella forma:
τ = f1 ( γ ) + f 2 ( γ , γ2 ) + f 2 ( γ − α , γ1 ) − f 2 ( γ − α , γ2 )
(3.5)
dove f1 ed f2 sono le due funzioni:
f1 ( γ ) = A * γ n
(3.6)
⎛ γ ⎞
f 2 ( γ , γ ) = m * A * γ n *ln ⎜ 1 + ⎟
⎝ B⎠
(3.7)
3.2.3 Modello di resistenza di Johnson e Cook
Il modello fenomenologico utilizzato nel presente lavoro è stato presentato da Johnson e
Cook nel 1983, [8], nella forma:
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m )
in cui ε è la deformazione plastica equivalente, ε* =
(3.8)
ε
è la velocità di deformazione
ε0
plastica adimensionalizzata per ε0 = 1.0s −1 e T* è la temperatura omologa:
T* =
T − Tro om
Tmelt − Tro om
38
(3.9)
dove T indica la temperatura assoluta, Troom la temperatura ambiente e Tmelt la
temperatura di fusione. Le costanti A, B, n, C ed m sono costanti dipendenti dal
materiale. L’espressione nel primo gruppo di parentesi esprime il valore della tensione
in funzione della deformazione, quindi la legge d’incrudimento, che si ha per una
velocità di deformazione pari a quella di riferimento ed un valore della temperatura
omologa nulla. Le espressioni nel secondo e nel terzo gruppo di parentesi esprimono,
rispettivamente, l’effetto della velocità di deformazione e quello della temperatura sulla
risposta meccanica dei materiali. Il punto di forza di tale modello è dato dalla possibilità
di trattare in modo disaccoppiato gli effetti dovuti alle tre variabili di deformazione,
velocità di deformazione e temperatura. Ciò, oltre a rendere molto semplice
l’implementazione del modello in qualunque codice numerico in commercio, permette
la caratterizzazione del materiale con un numero limitato di prove meccaniche. Bisogna
comunque sottolineare che la proporzionalità, espressa dall’equazione (3.8), della
tensione di snervamento con il logaritmo della velocità di deformazione non permette,
come illustrato in Figura 3.3, una corretta descrizione della risposta meccanica del
materiale in regimi di velocità di deformazione superiori a 104 s-1.
3.3 Modelli di Danneggiamento duttile nei metalli
La rottura duttile, se pur limitata all’ambito dei metalli, è un fenomeno estremamente
ampio e complesso. Per decenni, si è pensato alla rottura come ad un fenomeno
indipendente dalla storia dei processi di sforzo e deformazione che hanno luogo nel
materiale. Il comportamento di questo, cioè, non subiva modificazioni di sorta fino
all’improvvisa incapacità di sostenere i carichi. Le teorie di rottura, ad esempio, sono il
tentativo d’identificare il valore del carico massimo ammissibile senza interessarsi ai
meccanismi specifici di rottura. Anche se la rottura fragile ha ricevuto grande attenzione
dall’inizio del ventesimo secolo, la rottura duttile è stata studiata in dettaglio solo a
partire dagli anni sessanta.
McClintock, [9], e Rice e Tracy, [10], sono stati i primi ad identificare nel processo di
nucleazione e crescita dei microvuoti, correlate all’aumento del livello di deformazione,
il micromeccanismo responsabile della rottura duttile. Da allora, sono stati proposti un
gran numero di modelli di rottura, che, classicamente, sono suddivisi in “abrupt criteria”
e modelli “nucleation and growth (NAG)”.
39
Per i primi, a rottura improvvisa, questa avviene istantaneamente quando una variabile
interna ovvero una variabile di stato, raggiunge, in un punto, il valore critico. In tali
modelli il danno, anche se è accumulato durante la storia delle deformazioni, non è
accoppiato alle altre variabili costitutive. Questo è un modello tipico per la rottura dei
materiali fragili, per cui si ha rottura quando si raggiunge il valore critico dello sforzo
ovvero dell’intensità del campo di sforzo.
Per i modelli NAG, invece, l’attivazione dei danneggiamento è causa di una
modificazione delle proprietà meccaniche del materiale. La rottura è vista come il
risultato di un progressivo deterioramento del materiale e della sua capacità di sostenere
i carichi. La variabile, accoppiata alle altre variabili interne, che tiene in conto tale
deterioramento è comunemente indicata come danno e richiede la definizione di una
legge di evoluzione cinetica.
Gli “abrupt criteria” sono, di solito, facilmente implementabili nei codici numerici, ma,
di contro, risentono di una scarsa trasferibilità dimensionale e geometrica. Nella
dinamica dell’impatto, tali criteri sono stati largamente utilizzati con la giustificazione
del fatto che i fenomeni dinamici avvengono tanto rapidamente da confinare gli effetti
associati in volumi limitati e che, quindi, gli eventuali accoppiamenti del danno alle
altre variabili interne potessero essere trascurati. Tali modelli, però, poiché sono spesso
di natura fenomenologica, richiedono una caratterizzazione dei parametri a posteriori
che riduce fortemente l’effettiva capacità di previsione della rottura.
In questo contesto, ad esempio, il valore critico della pressione in tensione è
comunemente utilizzato per prevedere la rottura per spall in un Flyer Plate Impact Test.
La determinazione del valore critico, caratteristico del materiale che si sta investigando,
richiede l’effettuazione di un certo numero di prove a differenti velocità.
L’identificazione avviene comparando lo spall signal del profilo di velocità risultante,
con quello calcolato. Tale modo di operare non tiene in nessuna considerazione l’effetto
della triassialità dello stato di sforzo sul processo di rottura duttile.
Hancock e Mackenzie, [11], e Hancock e Brown, [12], evidenziarono che la triassialità
dello stato di sforzo (Triaxiality Factor, TF) ha un ruolo considerevole nel ridurre la
capacità di deformarsi del materiale. Proposero, allora, un modello per cui si ha rottura
quando in un punto del materiale si raggiunge un valore critico di deformazione, che
40
dipende dalla multiassialità dello stato di sforzo secondo la relazione:
⎛ 3σ
εf = α exp ⎜⎜ − m
⎜⎝ 2 σeq
⎞⎟
⎟
⎠⎟
(3.10)
in cui σm indica la pressione idrostatica, σeq la tensione equivalente di Mises, e α è una
costante dipendente dal materiale il cui valore può essere identificato in una prova a
sforzo uniassiale (TF= σm/ σeq=1/3).
Ad un’espressione simile è giunto, in modo indipendente, Manjoine, [13], interpolando
dati sperimentali per un certo numero di acciai:
εf = εuniaxial
2(
f
1−3 σm σeq )
(3.11)
In Figura 3.4 sono messe a confronto le due relazioni, che prevedono una riduzione
della duttilità al crescere della triassialità dello stato di sforzo, insieme ai dati
sperimentali relativi all’acciaio SA537 testato a differenti velocità di carico, [14].
Figura 3.4 - Riduzione della duttilità al crescere della triassialità dello stato di sforzo per l'aciaio
SA537, [14].
Johnson e Cook, [15],proposero un modello basato sul valore critico della
deformazione, in grado di considerare gli effetti della triassialità dello stato di sforzo,
41
della velocità di deformazione e della temperatura, secondo la seguente relazione:
⎛
σ
εf = ⎜⎜ D1 + D2 exp D3 m
⎜⎝
σeq
⎞⎟ ⎛
ε ⎞⎛
T − T0 ⎞⎟
⎟⎟ ⎜⎜ 1 + D4 ln ⎟⎟⎟ ⎜⎜ 1 + D5
⎟
ε0 ⎠⎝
Tmelt − T0 ⎠⎟
⎠⎝
(3.12)
Al fine di tenere in conto la storia delle deformazioni, proposero un criterio cumulativo
per cui si ha rottura quando la deformazione normalizzata, definita nell’equazione
(3.13), raggiunge il valore unitario:
D=
∑
∆εi
εif
(3.13)
Tali criteri, in genere, non sono, almeno in forma diretta, dipendenti dal tempo. Tuler e
Butcher, [16], osservando che lo sforzo in grado di causare una rottura con un impulso
di lunga durata è più basso di quello necessario con uno di breve durata, proposero la
seguente espressione:
tf
∫0
( σ − σ0 )λ dt ≥ Kc
(3.14)
in cui σ0 è il valore di soglia oltre il quale è attivato il criterio e tf è il tempo totale a
rottura. Si ha rottura quando, in un punto, l’integrale dato nell’equazione (3.14) supera il
valore di riferimento Kc. Tale modello è in grado di prevedere con buona
approssimazione la rottura per spall causata da un impulso triangolare, ma ha dei limiti
evidenti nell’incapacità di considerare gli effetti volumetrici e quelli legati alla
triassialità dello stato di sforzo.
I modelli NAG sono basati sull’assunzione che fenomeni irreversibili, che hanno luogo
durante il processo di deformazione, modificano la risposta del materiale e la sua
capacità di sostenere i carichi. Per un modello basato su questo approccio, è necessario
ridefinire le equazioni costitutive del materiale. Nel passato sono stati seguite,
prevalentemente, due strade differenti: quella dei “Porosity-based Models” e quella dei
“Contnuum Damage Models” (CDM).
Nei modelli basati sul concetto di porosità, questa è espressa con l’introduzione di una
variabile di porosità fittizia, correlata alla formazione di microvuoti nel materiale con la
deformazione plastica, che abbassa lo snervamento del materiale. Le equazioni
costitutive del materiale alla macro scala sono le equazioni elasto-plastiche standard del
materiale, ma il criterio di snervamento è modificato dalla porosità del materiale in
42
modo tale che quando questa raggiunge il valore critico, la funzione di snervamento
implode in un punto a sforzo nullo. Tale approccio è stato inizialmente formulato da
Gurson, [17],e, successivamente, è stato modificato da Tvergaard e Needleman, [18],
per tenere in conto l’effetto d’interazione tra i diversi vuoti. Needleman e Rice, [19],
modificarono il modello per la nucleazione di nuove famiglie di vuoti in fasi successive
del processo di deformazione. Anche se tale modello è largamente utilizzato in un gran
numero di applicazioni ed è disponibile nella maggior parte dei codici numerici
commerciali, è fortemente limitato da due fattori: richiede la conoscenza di un numero
eccessivo di parametri dipendenti dal materiale (fino a 9); non è trasferibile a differenti
condizioni geometriche e di vincolo, [20] e [21].
Curran et al. [22], facendo riferimento allo stesso approccio, proposero una legge di
evoluzione della porosità differente, assumendo una distribuzione esponenziale dei
vuoti rispetto alla loro dimensione. Seaman et al. [23], proposero una funzione di
distribuzione della nucleazione legata alla pressione tensile.
Nella CDM si definisce, in modo alternativo ai modelli basati sul concetto di porosità,
un set di equazioni costitutive per il materiale danneggiato. In questo approccio, il
danno è una delle variabili di stato. Assumendo l’esistenza di un potenziale di
dissipazione di danno, si può ricavare la legge cinetica di evoluzione del danno.
Lemaitre, [24], ha per primo definito il contesto costitutivo per il danno duttile nei
materiali. Successivamente, [25], sono state presentate altre forme del potenziale di
dissipazione del danno che portano a diverse leggi di evoluzione del danno con la
deformazione plastica. Nel prossimo paragrafo è presentato il modello proposto da
Bonora nel 1997, [2], che sarà utilizzato nel prosieguo della presente trattazione.
3.3.1 Modello di danno duttile non lineare
Il modello di danno utilizzato nel presente lavoro è stato sviluppato da Bonora, [2], nel
contesto della CDM, inizialmente proposta da Lemaitre, [24]. Le caratteristiche
principali del modello possono essere brevemente riassunte nei seguenti punti:
9 Il modello è derivato sotto le ipotesi dell’esistenza di un potenziale di
dissipazione del danno e dell’equivalenza delle deformazioni, che portano alla
definizione di tensione effettiva e all’accoppiamento tra danno e deformazione
plastica;
43
9 Il modello richiede un numero limitato di parametri dipendenti dal materiale,
soltanto quattro, tutti con un preciso significato fisico;
9 L’identificazione degli stessi può essere facilmente ottenuta con semplici prove
di tensione uniassiale su provini di geometria appropriata (clessidra, round
notched e prova di pura torsione);
9 I parametri sono caratterizzati da trasferibilità geometrica;
9 La formulazione del modello di danno è indipendente dal materiale, differenti
evoluzioni del danno con la deformazione plastica possono essere accuratamente
descritte con il medesimo potenziale di danno, utilizzando in modo appropriato
il set di parametri, Figura 3.5;
9 La formulazione proposta non presenta i problemi di localizzazioni tipici delle
formulazioni con softening.
1.0
Al2024-T3
0.9
Cu99.9%
0.8
AISI1045
0.7
Bonora's model
0.6
cr
D/D 0.5
0.4
0.3
0.2
0.1
0.0
-0.1
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
( ε- εth )/( εcr - εth )
Figura 3.5 - Evoluzione del danno, normalizzato rispetto al danno critico, in funzione della
deformazione plastica, per diversi tipi di metalli.
Di seguito è riportato il set base di equazioni costitutive:
Decomposizione delle deformazioni totali,
εTij = εije + εijp
44
(3.15)
Velocità delle deformazioni elastiche,
εije =
1 + ν σij
ν σkk
−
δ
E 1 − D E 1 − D ij
(3.16)
Legge evolutiva delle deformazioni plastiche,
∂ fp
3 sij
εijp = λ
= λ
∂σij
2 σeq
(3.17)
Definizione del moltiplicatore plastico
∂ fp
r = −λ
= λ = p
∂R
(3.18)
in cui le equazioni (3.17) e (3.18) sono quelle della plasticità standard, mentre la legge
cinetica di evoluzione del danno è data da:
1
α −1 p
⎛σ ⎞
( D − D0 )α
∂f
D = −λ D = α ⋅ cr
⋅ f ⎜⎜ H ⎟⎟⎟ ⋅ ( Dcr − D ) α ⋅
⎜⎝ σeq ⎠
p
ln(εf / εth )
∂Y
(3.19)
dove
f
⎛ Ρ
⎜⎜
⎜⎝ σeq
⎞⎟ 2
⎛ Ρ
⎟⎟ = ( 1 + ν ) + 3 ⋅ ( 1 − 2ν ) ⋅ ⎜⎜⎜
⎠ 3
⎝ σeq
⎞⎟2
⎟
⎠⎟
(3.20)
esprime l’effetto della triassialità degli sforzi. I parametri di danno richiesti sono: εth , la
soglia di deformazione alla quale i processi di danneggiamento hanno inizio; εf , la
deformazione teorica a rottura uniassiale; Dcr , il danno critico al quale si ha la rottura e
α , l’esponente di danno che determina la forma della curva dell’evoluzione del danno
con la deformazione plastica.
45
Bibliografia
[1]
Pirondi A. e Bonora N., Modeling Ductile Damage Under Fully Reversed
Cycling, Computational Solids Mechanics, 2002.
[2]
Zukas, J. A., High Velocity Impact Dynamics, John Wiley & Sons, New York,
1990.
[3]
Zukas, J. A, Nicholas, T., Swift, H. F., Greszczuc, L. B. and Curran, D. R., Impact
Dynamic, John Wiley & Sons, New York, 1992.
[4]
Zener C., e Hollomon J.H., J. Appl. Ph., 15, pp. 22-32, 1944.
[5]
Zerilli, F.J. e Armstrong, R.W., Dislocation-Mechanics-Based Constitutive
Relations for Material Dynamics Calculations, J. Appl. Ph., 61, pp. 1816-1825,
1987.
[6]
Hartley, K.A. e Duffy, J., Strain Rate and Temperature History Effects During
Deformation of FCC and BCC Metals, Mechanical Properties at High Rates of
Strain, Conf. Ser. No. 70, Institute of Physics, London, pp. 21-30, 1984.
[7]
Campbell, J.D., Eleiche, A.M. e Tsao, M.C., Strength of Metals and Alloys at
High Strain and Strain Rates, Fundamental Aspects of Structural Alloy Design,
Plenum, New York, pp.545-561, 1977.
[8]
Johnson G.R. e Cook W.H., A Constitutive Model and Data for Metals Subjected
to Large Strains, High Strain Rates and High Temperature, Proc. 7th Inter. Symp.
Ballistic, The Hague, Netherlands, pp. 541-547, 1983.
[9]
McClintock, F.A., A criterion for Ductile Fracture by the Growth of Holes, J.
Appl. Mech., 35, pp. 363-371, 1968.
[10] Rice, J.R., and Tracy, D.M., On the Ductile enlargement of Voids in Triaxial
Stress Fields, J. Mechanics Physics Solids, 17, pp. 210-217, 1969.
[11] Hancock, J.W. and Meckenzie, A.C., On the Mechanisms of Ductile Failure in
High Strength Steel Subjected to Multiaxial Stres-States, J. of Mech. Phy. Sol., 24,
pp. 147-169, 1976.
[12] Hancock, J.W. e Brown, D.K., On the Role of Strain and Stress State in Ductile
Failure, J. of Mech. Phy. Sol., 31, pp.1-24, 1983.
46
[13] Manjoine, M.J., Welding Research Supplement, 50s-57s, 1982.
[14] Bonora, N., e Milella, (2000) P.P., USAF/AFRL Contract no F61775-00-WE029
Final Report.
[15] Johnson G.R. e Cook W.H., Fracture Characteristics of Three Metals Subjected
to Various Strains, Strain Rates, Temperature and Pressures, Engng. Fract.
Mech., 21, pp. 31-48, 1985.
[16] Tuler, F.R., e Butcher, B.M., A Criterion for the Time Dependence of Dynamic
Fracture, Int. J. Fract. Mech., 4, pp. 431-437, 1968.
[17] Gurson, A.L., Continuum Theory of Ductile Rupture by Void Nucleation and
Growth: Part I – Yield Criterion and Flow Rules for Porous Ductile Materials, J.
Engn. Mat. Tech., 99, pp. 2-15, 1977.
[18] Tvergaard, V. e Needleman, A., , Acta Metall., 32, p.157-169, 1984.
[19] Needleman, A., e Rice, J. R., Limits to ductility by plastic flow localization in
Mechanics of Sheet Metal Forming, Koistinen and Wang Eds., Plenum Publ.
Company, New York, pp. 237-265, 1978.
[20] Gao, X., Faleskog, J., Shih, C.F., e Dodds, R. H., Eng. Fracture Mechanics, 59,
No. 6, pp. 761-777, 1998.
[21] Brocks, W., Klingbeil, D., Kunecke, G., e Sun, D.Z., in Constraint Effects in
Fracture: Theory and Application, ASTM STP 1244, Kirk and Bakker Eds.,
American Society for Testing and Materials, Philadelphia, 1994.
[22] Curran, D. R., Seaman, L. e Shockey, D.A., in Shock Waves and High Strain
Rates Phenomena in Metals, Ed. Meyers M.A. and Murr L.E., Plenum, NewYork, 1981.
[23] Seaman L., Curran, D.R., e Shockey, D.A., Computational Models for Ductile and
Brittle Fracture, J. Appl. Phys., 47, pp. 4814-4826, 1976.
[24] Lemaitre, J., A Course on Damage Mechanics, Springer-Verlag, Berlin, 1992.
[25] Chandrakanth, S. and Pandey, P.C., Int. J. Fracture, 60, R73-R76, 1993.
[26] Bonora, N., A nonlinear CDM model for ductile failure, Engineering Fracture
Mechanics,
58,
pp.
47
11-28,
1997.
4 Strumenti di simulazione numerica per l’analisi dei
fenomeni dinamici
4.1 Introduzione
Gli strumenti analitici sono molto utili per sviluppare una comprensione, che potrebbe
essere definita intuitiva, dei fenomeni fisici che si manifestano durante i processi
d’impatto e permettono di valutare, con senso critico, i risultati sperimentali.
Consentono, anche, di fare delle previsioni, di sistemi molto semplici, purché non
vengano violate le ipotesi semplificative utilizzate per la derivazione dei modelli stessi.
Se si vogliono superare tali limitazioni, però, è indispensabile ricorrere agli strumenti
della simulazione numerica. Il loro utilizzo, ad oggi, ha il limite maggiore nella carenza
di modelli costitutivi in grado di descrivere correttamente il comportamento meccanico
dei materiali in regime dinamico, nondimeno tali strumenti sono fortemente consolidati,
permettono di calcolare correttamente i fenomeni di propagazione ondosa tipici dei
fenomeni transitori di geometrie comunque complesse, di risolvere problemi accoppiati
meccanici, termici, elettromagnetici etc., [3]. Il loro utilizzo non può, però, prescindere
dall’esperienza dell’operatore e da una profonda conoscenza sia delle problematiche che
si stanno analizzando, sia degli strumenti che si stanno adoperando.
In commercio esiste un gran numero di codi numerici per trattare i processi d’impatto,
con caratteristiche diverse: codici lagrangiani, euleriani, SPH, etc. A volte lo stesso
codice integra, per garantire una maggiore flessibilità, diversi approcci, come nel caso
dei codi che permettono l’interazione, nella medesima analisi, di griglie lagrangiane ed
euleriane. Solo con un’adeguata conoscenza di tutte queste caratteristiche si può
scegliere correttamente e in modo consapevole lo strumento più adeguato alle proprie
esigenze.
Una distinzione netta, in merito al metodo di risoluzione delle equazioni differenziali,
esiste tra i codici impliciti e quelli espliciti. I primi, più tradizionali, nati per la
risoluzione di analisi quasistatiche, richiedono l’inversione della matrice di rigidezza,
questo garantisce una maggiore affidabilità dei risultati ottenuti a discapito di una
maggiore difficoltà nel raggiungere la convergenza e di una più bassa velocità di
48
calcolo. I codici espliciti sono i codici più diffusi per le analisi dinamiche, sono dedicati
alla risoluzione delle problematiche d’impatto. Sono molto “robusti” in relazione alla
loro capacità di raggiungere la convergenza e permettono di effettuare analisi, anche
molto complesse, tridimensionali, con numerosi corpi a contatto, in un tempo
relativamente contenuto. Tale facilità di convergenza, però, impone, all’operatore, una
maggiore attenzione nel controllo dei risultati ottenuti, perché l’accumulo di errori nel
processo di integrazione può portare a stime, quantitative e qualitative, del tutto
sbagliate.
Di seguito sono presentati i due codici commerciali utilizzati per le analisi delle
configurazioni sperimentali investigate: il codice implicito, agli elementi finiti
MSC.Marc e il codice esplicito Autodyn.
4.2 Analisi dinamica in MSC.Marc
Il codice agli elementi finite MSC.Marc permette di effettuare analisi dinamiche di
diverso tipo, [2]:
Analisi agli autovalori;
Analisi transiente;
Risposta armonica;
Risposta spettrale.
Il programma utilizza due metodi per l’estrazione degli autovalori e tre operatori per
integrazione temporale. Possono essere trattate non linearità, dovute al materiale, alla
geometria e alle condizioni al contorno. I problemi lineari possono essere risolti con una
sovrapposizione modale ovvero con un’integrazione diretta. I problemi non lineari,
invece, possono essere risolti esclusivamente con i metodi d’integrazione diretta. Oltre
alle masse distribuite è possibile utilizzare masse concentrate associate ad ognuno dei
gradi di libertà del sistema. Lo smorzamento numerico può essere utilizzato sia per le
analisi modali sia per quelle transienti. Si possono applicare condizioni iniziali non
uniformi di spostamento o velocità, così come forze o spostamenti dipendenti dal
tempo.
49
4.2.1.1 Integrazione Diretta
L’integrazione diretta è un metodo numerico per risolvere le equazioni del moto di un
sistema dinamico, che può essere utilizzato sia per problemi lineari sia per quelli non
lineari. Per le analisi transienti, MSC:MARC offre i tre operatori d’integrazione diretta
di seguito riportati:
‰ Newmark-β Operator;
‰ Houbolt Operator;
‰ Central Difference Operators.
Tutte le tecniche d’integrazione diretta sono imprecise e presentano almeno uno dei
seguenti inconvenienti:
conditional stability;
artificial damping;
phase errors.
Newmark-β Operator
Tale operatore è, probabilmente, il metodo d’integrazione diretta più popolare e
utilizzato per le analisi agli elementi finiti. Per i problemi lineari è incondizionatamente
stabile e non presenta smorzamenti numerici. In problemi non lineari possono nascere
delle instabilità che possono essere superate con uno smorzamento adeguato o
riducendo l’intervallo del tempo d’integrazione. La procedura supporta, infatti, la
possibilità di variazione di tale intervallo e l’utilizzo di un controllo adattativi dello
stesso.
Si considerino le equazioni del moto di un sistema strutturale scritte in forma matriciale:
Ma + Cv + K u + F = 0
(4.1)
in cui M, C, e K indicano, rispettivamente, le matrici di massa, di smorzamento e di
rigidezza, e a, v, u, e F sono i vettori di accelerazione, velocità, spostamento e forza.
La forma generalizzata del Newmark-β operator è data da:
u n +1 = u n + ∆ t ⋅ v n +
v
n +1
n
( 21 − β ) ∆t
n
2
⋅ a n + β ⋅ ∆t 2 ⋅ a n +1
= v + ( 1 − γ ) ∆t ⋅ a + γ ⋅ ∆t ⋅ a
50
n +1
(4.2)
in cui l’apice n indica l’intervallo temporale ennesimo. Le equazioni della dinamica che
corrispondono alla particolare forma della legge trapezoidale,
γ =
1
2
β =
1
4
(4.3)
risultano in:
( ∆4t
2
M +
)
(
)
2
4 n
C + K ∆u = F n + 1 − R n + M a n +
v + Cv n
∆t
∆t
(4.4)
in cui R è la forza interna data da:
R=
∫V βT σdV
(4.5)
L’equazione (4.4) permette di ottenere la soluzione implicita del problema nella forma:
u n +1 = u n + ∆u
(4.6)
È bene sottolineare che la matrice dell’operatore include il termine la matrice di
rigidezza tangente, K, per cui, ogni non linearità comporta una riformulazione della
matrice dell’operatore. Questa, inoltre, dipende dall’intervallo d’integrazione e deve
essere, di conseguenza, ricalcolata per ogni variazione dello stesso. Nell’equazione (4.2)
, γ è un parametro che, se settato ad un valore diverso da 1 2 , introduce uno
smorzamento nella risposta. Questo può permettere, ad un operatore molto esperto, di
introdurre una viscosità artificiale che può essere necessaria, ad esempio, nella
simulazione di un fenomeno caratterizzato dalla propagazione di onde d’urto. Invece β è
il parametro caratteristico della formulazione di Newmark. Al variare del suo valore
infatti, tale metodo diviene equivalente a formulazioni di volta in volta differenti, ad
esempio:
β=
1
- accelerazione costante nell’incremento temporale
4
β=
1
- accelerazione lineare
6
β=
1
- variazione a gradino dell’accelerazione
8
β = 0 - formulazione esplicita del secondo ordine.
51
4.2.2 Houbolt Operator
Tale operatore possiede la medesima stabilità incondizionata per I problemi lineari del
Newmark-β operator. Inoltre è caratterizzato da un elevato smorzamento numerico che
lo rende molto stabile anche per i problemi non lineari. La stabilità, infatti, cresce
all’aumentare del dell’intervallo d’integrazione. Di contro, però, l’elevato smorzamento
può portare, per intervalli d’integrazione molto lunghi, a soluzioni non accurate.
L’Houbolt operator è basato sull’utilizzo di un’interpolazione cubica per i valori a
quattro tempi differenti, tre determinati in precedenza e il corrente incognito. Ciò risulta
nelle equazioni:
v n +1 =
( 116 u
n +1
)
(4.7)
1
∆t 2
(4.8)
3
1
1
− 3u n + u n −1 − u n −2 ⋅
2
3
∆t
e
a n +1 = ( 2u n +1 − 5u n + 4u n −1 − u n −2 )
Sostituendo le equazioni (4.7) e (4.8) nell’equazione del moto, si ottiene:
( ∆2t
)
11
C + K ∆u =
6∆t
(4.9)
n
n −1
n −2
u
−
u
+
u
M
3
4
(
)
C
7
3
1
= F n +1 − R n +
+
u n − u n −1 + u n −2
∆t 6
2
3
∆t 2
2
M +
(
)
che fornisce uno schema di soluzione “implicito” per l’equazione (4.1), dalla quale si
ottengono vn+1 e an+1.
4.2.3 Central Difference Operator
Tale operatore esplicito è stabile solo in modo condizionale e il programma calcola,
automaticamente, il massimo intervallo ammissibile del tempo d’integrazione. Tale
metodo non è applicabile a strutture di tipo guscio o trave, perché le elevate frequenze
risultano in un limite di stabilità estremamente piccolo, è invece molto utile per l’analisi
dei fenomeni di shock.
Il Central Difference Operator assume una legge di variazione dello spostamento
rispetto al tempo, di tipo quadratico:
(
an = v
n + 12
−v
n − 12
52
) ( ∆t )
(4.10)
(
) ( ∆t )
(4.11)
a n = ( ∆u n +1 − ∆u n ) ( ∆t 2 )
(4.12)
∆u n = u n − u n −1
(4.13)
vn = u
n + 12
−u
n − 12
in modo che:
in cui:
Nella forma più generale la soluzione è data da:
M
M
n − 12
n +1
n
= F n − Rn +
2 ∆u
2 ∆u − Cv
∆t
∆t
(4.14)
4.2.4 Damping
Il damping, o smorzamento numerico, riproduce, in un’analisi dinamica transiente, la
dissipazione di energia all’interno del sistema. In Marc sono previste, per tale analisi,
due tipologie di smorzamento: il modal damping, per il metodo delle sovrapposizioni
modali, e il Rayleigh damping, per l’integrazione diretta Ad ogni incremento di tempo,
il programma associa, ad ogni modo, la frazione di damping corrispondente.
L’integrazione è basata sull’assunzione che la matrice di smorzamento numerico del
sistema è costituita da una combinazione lineare delle matrici di massa e di rigidezza, e
che quindi non modifica i modi del sistema.
Il damping è utilizzato per smorzare le eccessive oscillazioni del sistema alle alte
frequenze. Poiché, al diminuire dell’intervallo d’integrazione, la matrice di damping
può causare uno smorzamento eccessivo, è opportuno utilizzare l’opzione che vede la
matrice di damping essere dipendente dall’intervallo d’integrazione, in modo che, anche
per intervalli molto piccoli, le frequenze più elevate possano essere correttamente
rappresentate. La matrice di damping è data dalla seguente relazione:
∑ { αiM i + ( βi + γi
n
C =
i =1
In cui:
C è la matrice di damping;
Mi è la matrice di massa dell’i-esimo elemento;
Ki è la matrice di rigidezza dell’i-esimo elemento;
53
) }
∆t
Ki
π
(4.15)
αi è il coefficiente smorzamento di massa sull’i-esimo elemento;
βi è il coefficiente smorzamento di rigidezza sull’i-esimo elemento;
γi è il coefficiente che rende lo smorzamento numerico proporzionale a ∆t;
∆t è l’intervallo del tempo d’integrazione;
Se, per l’intera struttura, sono usati gli stessi valori per I coefficienti di smorzamento,
l’equazione (4.15) risulta in una formulazione del damping equivalente a quella di
Rayleigh.
4.3 Analisi dinamica in Autodyn
Il codice numerico Autodyn, della Century Dynamics, utilizza tecniche alle differenze
finite, ai volumi finiti e agli elementi finiti per risolvere una grande varietà di problemi
non lineari nella dinamica sia dei solidi sia dei fluidi; i codici di tale categoria vengono
spesso indicati col termine “Hydrocode”. Essi possono essere utilizzati per studiare
fenomeni fortemente dipendenti dal tempo e con non linearità dovute alla geometria
(grandi spostamenti e grandi deformazioni) e al materiale (plasticità, incrudimento,
softening, danneggiamento, equazioni di stato, etc.). Autodyn incorpora diversi
processori numerici ognuno dei quali è ottimizzato per risolvere il problema in
determinati domini: strutture, fluidi, gas, etc. L’accoppiamento, nello spazio e nel
tempo, dei diversi domini permette di raggiungere la soluzione ottimale al problema. I
processori numerici inclusi in Autodyn sono i seguenti, [3]:
‰ processore lagrangiano, per la modellazione dei solidi continui e delle strutture;
‰ processore euleriano, per la modellazione dei fluidi e delle grandissime
distorsioni;
‰ Arbitrary Lagrange Euler (ALE), specifico per la modellazione dei flussi;
‰ processore shell, per la modellazione di elementi strutturali sottili;
‰ Smooth Particle Hydrodynamics (SPH).
Tutti i processori elencati utilizzano un metodo di integrazione nel tempo di tipo
esplicito. In tutte le simulazioni relative al presente lavoro di tesi si è utilizzato un
processore di tipo lagrangiano, il cui schema è stato derivato dal metodo utilizzato da
54
Wilkins, , nel codice HEMP. Rispetto ad un approccio euleriano, una formulazione
lagrangiana è, dal punto di vista computazionale, più veloce, non dovendo risolvere il
calcolo del trasporto di materiale attraverso la mesh. Permette, inoltre, di trattare più
facilmente le interfacce tra i materiali, le superfici libere e l’effetto della storia sul
comportamento del materiale. Un esempio del modo di operare di una formulazione
lagrangiana è riportato in Figura 4.1.
Figura 4.1 – Distorsione di una mesh lagrangiana, [3].
Lo svantaggio maggiore è dovuto alla perdita di accuratezza che, inevitabilmente si
accompagna ad un’eccessiva distorsione della mesh. Per superare tale problema,
Autodyn permette di effettuare il “Rezoning” della mesh in modo da attenuare le
distorsioni, Figura 4.2. Laddove ciò non sia sufficiente, come ad esempio in alcuni
fenomeni di penetrazione, la tecnica dell’erosione permette ad un operatore esperto di
ottenere una soluzione sufficientemente accurata.
Figura 4.2 – Tipica procedura di “rezoning”, [3].
4.3.1 Metodo d’integrazione esplicito
Il metodo per l’integrazione delle equazioni discretizzate è detto esplicito se gli
55
spostamenti al tempo t + ∆t , nel ciclo di calcolo, sono indipendenti dalle accelerazioni
allo stesso tempo. L’algoritmo alle differenze centrali del secondo ordine è uno degli
schemi d’integrazione più utilizzati. Sia data l’equazione del moto nella forma:
Mu + Ku = F ( t, u )
(4.16)
in cui M è la matrice delle masse, K è la matrice di rigidezza, u il vettore
spostamento, u l’accelerazione e F il vettore delle forze che include i carichi
meccanici, termici e le pseudoforze dovute alle non linearità geometriche e del
materiale. Le velocità e gli spostamenti possono essere espressi, in funzione del tempo,
nella forma:
(
u t +
)
∆t
1
[ u ( t + ∆t ) − u ( t ) ]
=
2
∆t
u ( t ) =
1
∆t
(
) (
)
∆t
∆t ⎤
⎡
− u t −
⎢ u t +
⎥
2
2 ⎦⎥
⎣⎢
(4.17)
(4.18)
Combinando le equazioni (4.16), (4.17) e (4.18) si ricava la relazione:
Mu ( t + ∆t ) = ( ∆t )2 F ( t ) + ⎡⎣ 2M − ( ∆t )2 K ⎤⎦ u ( t ) − Mu ( t − ∆t )
(4.19)
Ad ogni intervallo di tempo sono noti le velocità e gli spostamenti, da cui possono
essere derivate le velocità di deformazione e le deformazioni stesse. Si ripete quindi la
procedura per determinare le accelerazioni e le velocità all’intervallo di tempo
successivo.
La risposta può diventare instabile se l’intervallo di tempo scelto non è sufficientemente
piccolo. Per i problemi non lineari non esiste un criterio di stabilità rigoroso, ma
l’esperienza ha dimostratati che si ottiene un buon risultato se:
∆t =
kl
c
(4.20)
in cui l è la dimensione minima degli elementi della mesh, c è la velocità del suono nel
mezzo, k è un coefficiente, minore di uno, generalmente compreso tra 6,0 e 9,0.
4.3.2 Viscosità artificiale
Al fine di limitare le discontinuità correlate alla comparsa di onde d’urto, si introduce
nella soluzione un termine viscoso artificiale. Von Neumann e Richtmeyer, [4],
introdussero un termine, quadratico nella velocità di deformazione, da sommare al
56
valore della pressione idrostatica, nei bilanci di energia e quantità di moto. Nel 1980,
Wilkins, [5], propose un ulteriore termine, lineare nella velocità di deformazione, per
smorzare le piccole oscillazioni ad alta frequenza che si hanno a valle dello shock. Tale
formulazione è utilizzata in gran parte dei codici espliciti in commercio, compreso
Autodyn, nella forma:
⎡⎛
⎛V ⎞
⎛V ⎞⎞⎤
q = ρ ⎢ ⎜⎜CQd ⎜⎜ ⎟⎟ − C Lc ⎜⎜ ⎟⎟
⎟⎟ ⎥
⎜⎝V ⎠⎟
⎟⎟ ⎥⎦
⎢⎣ ⎜⎝
⎝⎜V ⎠⎠
per
q =0
per
V
<0
V
(4.21)
V
<0
V
in cui CQ e C L sono costanti, ρ è la densità, d è una lunghezza caratteristica, c è la
velocità del suono nel mezzo e
V
è la variazione volumetrica.
V
4.4 Implementazione numerica del modello di danno non lineare
L’implementazione numerica del modello di danno è stata effettuata su entrambi i codi
di calcolo utilizzati in questo studio: MSC.Marc e Autodyn. La stessa è avvenuta
attraverso l’utilizzo di “user subroutines” direttamente collegate al programma
principale. La formulazione del modello permette una facile implementazione per via
del fatto che i potenziali di plasticità e di dissipazione del danno sono disaccoppiati. Le
equazioni della plasticità, di conseguenza, sono le equazioni standard già implementate
nel programma principale. Poiché la legge di evoluzione del danno, in accordo con
l’equazione (3.19), è funzione dell’ammontare di danno accumulato, della deformazione
plastica accumulata e della triassialità dello stato di sforzo, essa deve essere integrata
per l’incremento di deformazione plastica corrente. A tale scopo si utilizza lo schema
d’integrazione numerico Runge-Kutta. Il danno è calcolato ad ogni intervallo temporale,
per ogni punto di gauss. Quando, per tutti i punti di gauss di un elemento, si raggiunge il
valore del danno critico, questo viene rimosso e sforzi e deformazioni vengono
rilasciati. Tale procedura può essere causa di instabilità numeriche, che possono però
essere facilmente superate ricorrendo ad un intervallo del tempo d’integrazione
relativamente piccolo, in modo che per ogni intervallo non venga rimosso più di un
elemento. Questo permette al sistema di ristabilire gli equilibri e di evitare una
57
propagazione degli errori. Lo schema logico seguito è riportato nel diagramma di flusso
di Figura 4.3.
INPUT
EXIT
σ ij( t ) ε ijp ( t )
YES
(t)
∆σ ij( t ) , ∆ε ijp ( t ) , σ H( t ) ,σ eq
ε ijp + ( t ) = ε ijp + ( t −1 ) + ∆ε ijp ( t )
Integrate ∆D+
Runge-Kutta
D+(t)= D+(t-1)+ ∆D+
ε ijp + ( t ) = ε ijp + ( t −1 ) + ∆ε ijp ( t )
IF
NO
ε eqp + ( t ) = ε eqp + ( t −1 )
TF≤ 0
2 p(t ) p(t )
∆ε ij ∆ε ij
3
IF
⎛σ ⎞ 2
⎛σ ⎞
f ⎜⎜ H ⎟⎟ = (1 + ν ) + 3 ⋅ (1 − 2ν ) ⋅ ⎜⎜ H ⎟⎟
⎝ σ eq ⎠ 3
⎝ σ eq ⎠
2
YES
ε
p+
eq
≥ εth and Dflag( m ) = 0
NO
EXIT
Remove element.
IF
D+(t)≥ Dcr
YES Dflag(m)=1
SET
σ = 0 ,εijT (t) = 0
(t)
ij
NO
Update variables
Update stiffness matrix
~
E = E ⋅(1 − D )
Figura 4.3 - Schema logico per il calcolo numerico della variabile di danno. La procedura inizia alla
fine di ogni incremento, quando sono già state calcolate tutte le variabili globali e ed è ripetuta per
ogni punto di gauss di ogni elemento attivo.
58
Bibliografia
[1]
Zukas, J. A, Nicholas, T., Swift, H. F., Greszczuc, L. B. and Curran, D. R., Impact
Dynamic, John Wiley & Sons, New York, 1992.
[2]
MSC.Marc Volume A: Theory and User Information, Version 2005, U.S.A.,
2005.
[3]
Autodyn documentation, Theory Manual, Revision 4.3, Century Dynamics, 2003.
[4]
Von Neumann, J. e Richtmeyer, R.D., “A Method fort he Numerical Calculation
of Hydrodynamic Shocks”, J. Appl. Phys., 21, pp. 232-237, 1950.
[5]
Wilkins, M. L., “Use of Artificial Viscosity in Multidimenional Fluid Dynamic
Calculation”, J. Comp. Phys., 36, pp. 281-303, 1980.
59
5 Taylor Test
5.1 Analisi teorica del test di Taylor
Il test di Taylor è una tecnica sviluppata per determinare il valore della tensione di
snervamento di un materiale soggetto a carichi dinamici. Esso prevede che un provino
di forma cilindrica venga fatto impattare normalmente, a velocità nota, contro una
parete rigida e che si deduca la tensione di snervamento ricercata dalla velocità
d’impatto e dalla geometria iniziale e finale del provino. Taylor ha proposto [1], nel
1948, un’analisi semplificata del fenomeno, assumendo che:
‰ il materiale abbia un comportamento rigido perfettamente plastico ed
indipendente dalla velocità di deformazione, σ=σ(ε);
‰ la propagazione delle onde all’interno del cilindro sia monodimensionale;
‰ il flusso plastico sia incompressibile e la deformazione elastica sia trascurabile.
Figura 5.1 - Schematizzazione del cilindro di Taylor: (a) durante la deformazione; (b) al termine
della deformazione.
La Figura 5.1a mostra una schematizzazione del cilindro ad un certo punto durante la
prova. La regione deformata cresce con velocità pari alla velocità di propagazione
dell’onda plastica cp, mentre la porzione indeformata del cilindro, la cui lunghezza
istantanea è h, viaggia alla velocità decrescente v. Indicando con A0 l’area della sezione
iniziale del cilindro e con σy la tensione di snervamento, si possono scrivere:
l’equazione di conservazione della massa,
c p A = ( v + c p ) A0
60
(5.1)
l’equazione della conservazione della quantità di moto,
ρ (v + cp ) v = σ − σ y
(5.2)
e l’equazione del moto per la parte indeformata,
ρh
dv
= −σ y
dt
(5.3)
Assumendo, ancora, che la velocità di propagazione dell’onda plastica sia costante e che
la superficie libera del cilindro venga decelerata uniformemente, si può derivare la
formula di Taylor:
σy
(l − H )
= 0
2
ρ v0 2 ( l0 − l1 )
1
⎛l ⎞
ln ⎜ 0 ⎟
⎝H⎠
(5.4)
in cui l0 è la lunghezza iniziale del cilindro, l1 ed H sono, come mostrato in Figura 5.1b,
rispettivamente le lunghezze, rilevate al termine della prova, dell’intero cilindro e della
sua porzione indeformata.
Taylor ha introdotto un fattore correttivo nell’analisi in virtù del fatto che la
decelerazione del cilindro, in realtà, non avviene in maniera costante. Se si indica con
σ y il valore corretto di σ y determinato con la relazione precedente, si può scrivere:
⎛l ⎞
ln ⎜ 0 ⎟
σ y l0 − l1
⎝H⎠
−
=
2
σ y l0 − H ⎛
cp ⎞
⎜K − ⎟
a⎠
⎝
(5.5)
con:
a2 =
K=
2σ y
ρ
v0 + c p
a
(5.6)
(5.7)
Nel corso degli anni numerosi ricercatori hanno cercato di superare alcune delle ipotesi
semplificative, introdotte nella formulazione di Taylor, per ottenere una sua validità più
generale. Nel 1954, Lee e Tupper, [4], hanno presentato un modello che tiene in
61
considerazione la deformazione elastica. Raftopoulos e Davis [5], hanno generalizzato il
comportamento del materiale, includendo la deformazione elastica e il lavoro
d’incrudimento. Jones et al. [6], hanno proposto una nuova equazione del moto per la
parte indeformata del provino. Erlich et al. [7], nel 1981, hanno presentato una tecnica
alternative che prevede l’impatto simmetrico tra due cilindri, così da eliminare le
incertezze dovute all’indeterminatezza dovuta all’attrito sulla superficie d’impatto.
5.2 Simulazione numerica del Taylor test
Nella presentazione del test di Taylor si è visto come una delle assunzioni principali, su
cui è ricavata la relazione che permette di ricavare il valore della tensione di
snervamento del materiale in regime dinamico, è l’unidimensionalità dello stato di
sforzo. In realtà, proprio tale assunzione è stata in passato largamente contestata.
Wilkins e Guinam [8], evidenziarono la necessità di ricorrere ad un’analisi
bidimensionale per simulare correttamente il test in oggetto.
Nel presente studio, è stata effettuata un’estesa campagna numerica per verificare la
possibilità di riprodurre le caratteristiche salienti dell’esperimento e di stimare il ruolo
dei diversi parametri che caratterizzano il modello di resistenza del materiale in regime
dinamico. I risultati ottenuti sono stati sempre, per quanto possibile, messi a confronto
con risultati sperimentali disponibili in letteratura.
La modellazione agli elementi finite del Taylor Cylinder Impact Test, pur dando
l’impressione di essere estremamente semplice, nasconde, in realtà, una serie di aspetti
critici quali, il contatto tra cilindro e incudine, l’attrito tra le superfici d’impatto, le
dimensioni degli elementi, l’aspect ratio”, che rendono l’analisi molto delicata. Per
queste prime simulazioni numeriche, si è utilizzato il codice MSC.Marc facendo ricorso
alla schema di risoluzione diretto Newmark-β. Si è effettuata un’analisi parametrica
degli effetti associati al modello di resistenza utilizzato, nello specifico quello di
Johnson e Cook, [1]. Una notevole attenzione è stata dedicata alla corretta simulazione
dell’incudine, la cui influenza sulla bontà dei risultati non è, in letteratura,
opportunamente evidenziata. Le analisi hanno riguardato l’impatto di cilindri, a diverse
velocità, di rame OFHC, ferro ARMCO e acciaio AISI, le cui proprietà meccaniche
sono riportate in Tabella 5.1.
62
Tabella 5.1 - Propietà meccaniche dei materiali investigati e relativi parametri per il modello di
Johnson e Cook, [1].
Costanti per il modello di
Johnson e Cook
Proprietà Meccaniche
Materiale
Densità
[kg/m3]
Calore
specifico
[J/kg K]
Temperatur
a di fusione
(°K)
A
(MPa)
B
(MPa)
n
C
m
Acciaio
AISI 4340
7830
477
1793
792
510
.26
.014
1.03
Ferro
ARMCO
7890
452
1811
175
380
.32
.060
0.55
Rame
OFHC
8960
383
1356
90
292
.31
.025
1.09
L’analisi è stata effettuata in configurazione assialsimmetrica. Poiché gli elementi della
zona di contatto sono soggetti ad elevate deformazioni, la mesh iniziale è stata
realizzata con un aspect ratio rettangolare, in modo da prevenire schiacciamenti
eccessivi degli elementi della zona di contatto. Questa strategia consente di superare i
problemi derivanti da un’eccessiva distorsione degli elementi che possono risultare in
uno Jacobiano negativo e, quindi, all’interruzione della simulazione. In Figura 5.2 è
riportato un esempio della mesh utilizzata, con dimensioni dell’elemento di 0,5x0,25
mm2.
Figura 5.2 - Dettaglio della mesh nella zona d'impatto, per il rame OFHC.
63
Uno studio preliminare è stato effettuato modellando l’incudine come infinitamente
rigido, una strategia largamente utilizzata in letteratura, [1]. Tale assunzione si è, però,
rivelata essere inadeguata, perché porta alla generazione di disturbi ad altissima
frequenza all’interfaccia tra provino ed incudine che provocano chattering al contatto ed
altri problemi di natura numerica.
Per superare tali inconvenienti, l’incudine è stato modellato, più simile al caso reale,
come corpo deformabile a snervamento molto elevato. La mesh dell’incudine è, nella
zona del contatto, uniforme e formata da elementi quadrati, mentre nella parte
rimanente, dove è richiesta una minore accuratezza, gli elementi hanno forma
rettangolare e dimensioni via via crescenti. Per chiudere la mesh alle estremità ed
evitare che onde di riflessione possano interferire con il provino, sono stati utilizzati
elementi semi-infiniti a sei nodi e sei punti d’integrazione; l’elemento mostrato in
Figura 5.3, è formulato in modo tale da estendersi virtualmente all’infinito e da
considerare, lì, nulli gli spostamenti.
Figura 5.3 - Elemento semi-infinito usato nella modellazione dell'incudine.
Come sottolineato più volte, il fattore più critico nella simulazione dei processi dinamici
è dato dalla difficoltà di caratterizzare il materiale in modo adeguato. Per questo motivo
si sono analizzati, sfruttando la formulazione del modello Johnson e Cook, [1], gli
effetti che i diversi parametri, deformazione, velocità di deformazione e temperatura,
hanno sui risultati della simulazione numerica.
Il primo gruppo di simulazioni è stato effettuato senza tenere in considerazione il
danneggiamento del materiale ed utilizzando la legge costitutiva nella forma:
64
σ = A + Bε n
(5.8)
Si sono, in altre parole, trascurati gli effetti dovuti alla velocità di deformazione ed alla
temperatura. Le prove sono state effettuate per il rame OFHC (l0=25,4mm) a 190m/s,
per il ferro ARMCO (l0=12,6mm) a 279m/s e per l’acciaio AISI (l0=8,1mm) a 343m/s.
Le simulazioni numeriche forniscono, per tutte le configurazioni, deformazioni del
proiettile molto superiori a quelle rilevate sperimentalmente. La verifica può essere
effettuata confrontando gli accorciamenti, espressi come il rapporto tra la lunghezza del
cilindro al termine della prova e la lunghezza iniziale, per le diverse configurazioni. Per
il ferro il rapporto calcolato è pari a 0,58 contro lo 0,70 rilevato dalle sperimentazioni,
mentre per l’acciaio si ha uno 0,73 calcolato contro lo 0,8 misurato. Per il rame, inoltre,
l’eccessiva deformazione provoca l’interruzione del calcolo dopo un tempo di appena
0,2µs, a fronte di un tempo necessario alla conclusione del processo di circa 80µs,
Tabella 5.2.
Tabella 5.2 - Confronto tra gli accorciamenti calcolati, con diversi modelli di resistenza, e i dati
sperimentali.
Dati
sperimentali
Risultati numerici
σ = f (ε )
σ = f ( ε , ε ) σ = f ( ε , ε, T )
Rame OFHC
0.68
-
0.65
0.66
Ferro ARMCO
0.70
0.58
0.70
0.66
Acciaio AISI
0.73
0.80
0.77
0.61
Per il secondo gruppo di simulazioni, si è introdotto il modello di danno e si è utilizzata
la legge costitutiva nella forma:
σ = ( A + Bε n )(1 + C *ln ε* )
(5.9)
in modo da tenere in conto gli effetti dovuti alla velocità di deformazione ma non
ancora quelli legati alla temperatura. Le configurazioni per le simulazioni sono le stesse
utilizzate in precedenza ma i risultati ottenuti sono decisamente differenti. Un tipo di
legge di questo tipo è sicuramente più adeguata alla descrizione del comportamento
meccanico del materiale e ciò è confermato dal confronto tra gli accorciamenti calcolati
per i diversi provini e quelli misurati nelle sperimentazioni. Per il ferro ARMCO
65
l’accorciamento calcolato numericamente è di 0,7, pari a quello misurato; per l’acciaio
AISI è di 0,77, contro lo 0,8 misurato; per il rame OFHC, il rapporto calcolato è pari a
0,65 contro lo 0,68 misurato. Nelle Figura 5.4, Figura 5.5 e Figura 5.6 si può osservare
come i profili di deformazione per il rame e per l’acciaio corrispondano con ottima
approssimazione a quelli rilevati sperimentalmente. Questa è un’evidente conferma del
fatto che il materiale ha un comportamento meccanico differente a seconda che venga
Distanza dall'asse di simmetria [mm]
sollecitato staticamente ovvero in modo dinamico.
12
10
8
6
4
2
0
0
2
4
6
8
10
12
14
16
18
20
22
24
Posizione [mm]
Figura 5.4 - Profilo dela deformata per il rame OFHC (l0=25,4mm; V0=190m/s) ottenuta assumendo
σ = ( A + Bε n )(1 + C *ln ε* ) .
Distanza dall'asse di simmetria [mm]
12
10
8
6
4
2
0
0
2
4
6
8
10
12
14
1
Posizione [mm]
Figura 5.5- Profilo della deformata per il ferro ARMCO (l0=12,6mm; V0=279m/s) ottenuta
assumendo
σ = ( A + Bε n )(1 + C *ln ε* ) .
66
Distanza dall'asse di simmetria [mm]
12
10
8
6
4
2
0
0
2
4
6
8
10
12
1
Posizione [mm]
Figura 5.6 - Profilo della deformata per l’acciaio (l0=8,1mm; V0=343m/s) ottenuta assumendo
σ = ( A + Bε n )(1 + C *ln ε* ) .
Per l’ultimo gruppo di simulazioni si è provveduto a completare l’equazione costitutiva
aggiungendo il termine che tiene in conto gli effetti della temperatura, la relazione,
dunque, si presenta nella forma:
Distanza dall'asse di simmetria [mm]
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m )
(5.10)
12
10
8
6
4
2
0
0
2
4
6
8
10
12
14
16
18
20
22
24
Posizione [mm]
Figura 5.7 - Profilo della deformata per il rame OFHC (l0=25,4mm; V0=190m/s) ottenuta
assumendo
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) .
67
Distanza dall'asse di simmetria [mm]
12
10
8
6
4
2
0
0
2
4
6
8
10
12
14
1
Posizione [mm]
Figura 5.8 - Profilo della deformata per il ferro ARMCO (l0=12,6mm; V0=279m/s) ottenuta
assumendo
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) .
Distanza dall'asse di simmetria [mm]
12
10
8
6
4
2
0
0
2
4
6
8
10
12
1
Posizione [mm]
Figura 5.9 - Profilo della deformata per l’acciaio (l0=8,1mm; V0=343m/s) ottenuta assumendo
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) .
I risultati, per quanto riguarda i profili delle deformate, per il rame e l’acciaio
rimangono sostanzialmente gli stessi, per il ferro invece, la deformata ottenuta
utilizzando l’ultima relazione costitutiva proposta è decisamente diversa e molto più
simile a quella misurata.
68
Un tale risultato è in accordo perfetto con i dati sperimentali, reperibili in letteratura,
che hanno dimostrato, per mezzo di prove effettuate alla barra di Hopkinson, un’elevata
sensibilità del ferro alle variazioni di temperatura.
5.3 Analisi critica dei meccanismi di propagazione delle onde durante
il test di Taylor.
Come detto in precedenza nel test di Taylor, l’unidimensionalità dello stato di sforzo
non è verificata. Questo si può intuire, ad esempio, osservando la forma della deformata
del cilindro nella parte impattata, che non può essere giustificata a meno di ammettere la
presenza di deformazioni radiali.
Gli strumenti numerici sono stati allora utilizzati per investigare i meccanismi di
propagazione delle onde che si vengono a verificare durante il test. Le interpretazioni
date nel presente paragrafo derivano dalle osservazioni delle analisi numeriche già
presentate e da altre effettuate per la configurazione del RoR. In particolare si è fatto
riferimento alla configurazione adottata da Mayes et al., [10], dell’impatto simmetrico
di cilindri, di calibro 7,62 mm, di rame OFE, con due diverse dimensioni medie del
grano. Il materiale con grano medio maggiore, 75µm, è stato impattato a 300m/s e
392m/s, mentre il material con grano medio più fine, 40µm, è stato impattato a 233m/s.
Le simulazioni sono state effettuate col codice esplicito Autodyn, il modello è stato
realizzando due griglie, una per cilindro, in configurazione assialsimmetrica. In Tabella
5.3, i diametri finali delle superfici impattate sono confrontati con le misure
sperimentali, per tutte le velocità e le microstrutture, dimostrando, almeno per le
velocità più basse, un ottimo accordo.
Tabella 5.3 - Confronto tra i diametri calcolati delle superfici d'impatto e quelli misurati.
Grano
grande
300m/s
Grano fine 392m/s
Grano
fine
233m/s
Diametro calcolato
12.2mm
15.0mm
10.6mm
Risultati sperimentali, [10]
12.4mm
12.5mm
10.9mm
Nelle Figure 5.10 a, b, c e d, viene illustrato, in configurazione assialsimmetrica, il
processo di generazione, propagazione e sovrapposizione delle onde di pressione che si
69
viene a verificare nei primi 2µs del processo d’impatto. Al momento dell’impatto si
generano, all’interfaccia, nei due cilindri, onde di compressione che, secondo la
convenzione utilizzata nella dinamica dell’impatto, hanno segno positivo e sono
indicate in rosso nella prima delle Figure 5.10.
a)
b)
70
c)
d)
Figure 5.10 a, b, c e d - Generazione, propagazione e sovrapposizione delle onde di pressione in un
RoR test a diversi istanti di tempo durante il processo di deformazione.
Dopo appena 1µs, sono facilmente distinguibili le onde di rilascio che dal bordo esterno
propagano verso il centro del provino. Queste, sovrapponendosi nella regione prossima
all’asse di simmetria, dando luogo, negli istanti successivi, ad uno stato di sforzo
tensile. Tale stato di sforzo porta al distacco, per un breve intervallo di tempo, della
71
parte centrale della superficie d’impatto, appena distinguibile nella Figure 5.10 d.
Successivamente si ripristina uno stato compressivo generalizzato, che ha come effetto
il ristabilirsi del contatto, tra cilindro e incudine, lungo l’intera superficie d’impatto.
Solo dopo questa prima fase, quindi, lo stato di sforzo si avvicina molto allo stato di
unidimensionalità ipotizzato da Taylor.
La conseguenza di tale meccanica è che durante il test, la deformazione e la velocità di
deformazione nel cilindro non sono uniformi e, in ogni punto, non sono costanti, ma
variano col tempo. Questo implica che da tale test non è possibile, in nessun modo,
estrapolare il valore dello snervamento del materiale, come era stato inizialmente
proposto. Ad oggi il test è ancora largamente utilizzato perché, con un sistema
relativamente semplice, permette di raggiungere velocità di deformazione dell’ordine di
104 ÷ 105 s-1, fornendo dati che possono essere utilizzati per verificare le potenzialità
degli strumenti numerici di previsione e di nuovi modelli costitutivi.
5.4 Analisi numerica dei meccanismi di danneggiamento
In questo paragrafo viene ripresa l’analisi del test di Taylor con particolare attenzione ai
processi di danneggiamento. In modo simile a quanto proposto nei paragrafi precedenti,
attraverso un’analisi numerica delle medesime configurazioni, si sono investigati gli
effetti che, in modo indipendente, velocità di deformazione e temperatura, hanno sui
processi di danneggiamento.
In Figura 5.11 è riportata la mappa di danno ottenuta per il rame OFHC, assumendo la
legge sforzo deformazione (5.9) ed utilizzando il modello di danno non lineare
presentato, con i parametri riportati in Tabella 5.4. Si può osservare come un certo
numero di elementi, tra i quindici ed i venti, a seconda della configurazione esaminata,
sia eliminato una volta che questi abbiano raggiunto il valore di danno critico. Le
regioni maggiormente danneggiate corrispondo a quelle che risultano essere,
dall’osservazione delle foto dei cilindri impattati, maggiormente degradate e questo
conferma, qualitativamente, la validità del modello di rottura utilizzato.
La legge di Johnson e Cook nella formulazione completa, Eqn. (5.10), è stata utilizzata
per analizzare l’effetto della temperatura sul danneggiamento.
72
Tabella 5.4 - Parametri di danno per il rame OFHC.
εth
εf
0.1
3.2
α
Dcr
0.85 0.63
Figura 5.11 - Mappa di danno per il rame OFHC (l0=25,4mm; V0=190m/s) ottenuta assumendo
σ = ( A + Bε n )(1 + C *ln ε* ) .
Tensione equivalente di von Mises [MPa]
1000
σ = A + Bε n
σ = ( A + Bε n )(1 + C *ln ε* )
σ =( A+ Bεn )(1+ C*lnε*)(1−T*m)
800
600
400
200
0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Deformazione plastica equivalente
Figura 5.12 - Diagrammi tensione deformazione, al variare della legge costitutiva, per un punto
appartenente alla superficie di contatto per l’impatto di un cilindro di ferro ARMCO.
73
La Figura 5.12 esprime gli andamenti, ricavati dalle simulazioni degli impatti per il
ferro ARMCO, della tensione in funzione della deformazione per le tre diverse relazioni
costitutive utilizzate. E’ evidente il salto del valore della tensione di snervamento
causato dall’introduzione del termine legato alla velocità di deformazione, ma è altresì
interessante notare come l’aumento della temperatura, dovuto alla deformazione
plastica, abbassi, all’aumentare di quest’ultima, il livello dello sforzo. Questo comporta
una maggiore duttilità del materiale, che si manifesta in una riduzione del
danneggiamento dello stesso. Tale riduzione è messa in luce dal fatto che, per nessuna
configurazione, si ha l’eliminazione degli elementi per il raggiungimento del valore di
danno critico.
La Figura 5.13 rappresenta la mappa di danno, per il rame OFHC, rilevata dalla
simulazione numerica che si è effettuata tenendo in considerazione, nel legame
costitutivo gli effetti dovuti sia alla velocità di deformazione sia alla temperatura. E’
importante sottolineare come, per il rame e l’acciaio, la valutazione del profilo delle
deformate possa portare a ritenere ininfluente l’effetto della temperatura, mentre,
un’analisi del livello di danno raggiunto ne manifesta tutta la sua importanza.
Figura 5.13 - Mappa di danno per il rame OFHC (l0=25,4mm; V0=190m/s) ottenuta assumendo
σ = ( A + Bε n )(1 + C *ln ε* )(1 − T *m ) .
74
5.4.1 Effetto della dimensione del grano
I risultati sperimentali in [10] mostrano una particolarità che opportuno investigare in
modo approfondito. In accordo con la teoria, si osserva, per il materiale a più grande
grano medio, un aumento del danno con la velocità d’impatto. Il materiale a grana fine,
invece, manifesta, sorprendentemente, rispetto a quello a grana grande, un maggiore
ammontare del danno, a più bassa velocità d’impatto, Figura 5.15.
Tale risultato mostra chiaramente che, per quanto riguarda il danneggiamento, le due
diverse microstrutture hanno caratteristiche di danneggiamento differenti. La
dimensione del grano, come noto, influenza la resistenza del materiale. A grani più
piccoli corrispondono valori più elevati della tensione di snervamento come descritto
dalla relazione di Hall-Petch:
σ y = σ 0 + Kd −0,5
(5.11)
in cui σ 0 e K sono costanti che dipendono dal materiale e d è la dimensione media del
grano. Ad un aumento della tensione di snervamento, al diminuire delle dimensioni del
grano, corrisponde una riduzione della duttilità. Per quanto concerne il danneggiamento,
esso inizia con la nucleazione di microvuoti in prossimità delle inclusioni ovvero, per i
metalli puri come il rame in oggetto, ai bordi grano. In quest’ultimo caso, il danno è
dato dall’impossibilità di accomodare, alla mesoscala, le deformazioni imposte alla
macroscala. Di conseguenza, nelle microstrutture a grana fine, per le quali il moto delle
dislocazioni è fortemente vincolato, il danneggiamento duttile dovrebbe iniziare ad un
livello di deformazione plastica più basso rispetto al caso di grana più grossolana.
Quindi, per un dato livello di deformazione, ci si aspetta un ammontare di danno più
elevato. Si è pensato che un tale andamento possa essere governato dal valore della
deformazione di soglia e che, quindi, questa debba, in qualche modo, essere influenzata
dalla dimensione del grano.
Sono stati raccolti i valori della deformazione di soglia per diversi materiali metallici e
riportati, in funzione della dimensione media del grano, in Figura 5.14, che evidenzia
una chiara dipendenza tra εth e la dimensione media del grano, indipendentemente dal
materiale. Si è trovato, inoltre, che i dati rappresentati in Figura 5.14 sono molto bene
interpolati da una curva del tipo:
75
ε th = A ( d − d 0 )
0,5
(5.12)
in cui A è una costante e d 0 sembra indicare un limite inferiore di dimensione del
grano, al di sotto del quale i processi di danneggiamento dovrebbero essere inibiti a
causa della perdita di duttilità. Questo è in accordo con recenti osservazioni sperimentali
che indicano una variazione del tipo di rottura, da duttile a fragile, per dimensioni del
grano estremamente piccole.
0.24
0.22
Low carbon
Damage threshold strain
0.20
0.18
0.16
0.14
ARMCO
0.12
OFHC Cu
0.10
0.08
A533B
0.06
OFHC Cu
0.04
AISI4340
Al2024
W90
0.02
0.00
0
20 40 60 80 100 120 140 160 180 200 220 240
Figura 5.14 - Deformazione di soglia in funzione della dimensione media del grano.
Per quanto è stato illustrato, nelle simulazioni numeriche si sono utilizzati, per le due
diverse dimensioni medie del grano, i valori di deformazione di soglia ricavati dalla
relazione (5.12): ε th = 0,1 per d = 0, 75µ m e ε th = 0, 04 per d = 40 µ m .
La mappa di distribuzione del danno è stata analizzata e confrontata, in Figura 5.15, con
i risultati derivati dalle micrografie dei cilindri sezionati, il risultato trovato è in ottimo
accordo con i dati sperimentali disponibili, sia per la deformata finale sia per la mappa
di distribuzione del danno.
Per il materiale con grano medio maggiore, è correttamente previsto un aumento del
danno con la velocità d’impatto. Si evince, inoltre, che il danno è causato da grandi
deformazioni plastiche, che avvengono in prossimità della zona di contatto, tardi nel
processo di deformazione. Per il materiale a grana fine, caratterizzato da un valore più
76
basso della deformazione di soglia, il danno sembra essere dovuto dall’intenso impulso
di tensione che si genera nelle prime fasi del processo di deformazione, secondo il
meccanismo descritto al paragrafo 5.3. Tale processo è caratterizzato da un basso livello
di deformazione plastica e da elevata triassialità dello stato di sforzo, mostrando forti
similarità alla rottura per spall che si ha, ad esempio, nel Flyer Plate Impact Test.
Figura 5.15 - Deformate e mappe di danno calcolate a confronto con i risultati sperimentali.
77
Bibliografia
[1]
Zukas, J. A, Nicholas, T., Swift, H. F., Greszczuc, L. B. and Curran, D. R.,
Impact Dynamic, John Wiley & Sons, New York, 1992.
[2]
Zukas, J. A., High Velocity Impact Dynamics, John Wiley & Sons, New York,
1990.
[3]
Taylor, G. I., “The use of flat ended projectiles for determining dynamic yield
stress: 1. Theoretical considerations” in Proc. R. Soc. Lond. Ser. A., vol 194, pp.
289-300, 1948.
[4]
Lee, E., and Tupper, J. Appl. Mech., 63-70 (1954).
[5]
Raftopoulos, D., and Davids, N., AIAA J. 5, 2254 (1967).
[6]
Jones, S. E., Gillis, P. P., and Foster, J. C., Jr., J. of Appl. Phys. 61, 499-502
(1987).
[7]
Erlich, D. C., Shockey, D. A., and Seaman, L., “Symmetric Rod Impact
Technique for Dynamic Yield Determination”, in Shock Waves in Condensed
Matter-1981, AIP Conference Proceedings 78, Menlo Park, CA, 1981, pp. 402406.
[8]
Wilkins, M. L., and Guinam, M. W., J. Appl. Phys. 44, 1200-1206 (1973).
[9]
Johnson, G. R. and Cook, W. H., A costitutive model and data for metals
subjected to large strains, higi strain rates and high temperatures, Proc. 7° Int.
Symp. On Ballistics, pp. 541-547, Netherlands, 1983.
[10]
Mayes, J.L., Hatfield, S.L., Gillis, P.P. e House, J.W., Int. J. Impact Engng., 14,
pp.503-508, 1993.
[11]
Zener C., and Hollomon J.H., (1944), J. Appl. Ph., 15, pp. 22-32
[12]
Cottrell A.H., (1957)Proc. of Conference on Property of Materials at High Rates
of Strain, pp. 1-12.
[13]
Hollomon J. H., (1947), Trans. AIME, 171, pp. 535-545
[14]
Milella, P.P., (2001), Proc. of Shock Compression of Condensed Matter, AIP
Conf. Proc. 620, Atlanta.
78
[15]
McClintock, F.A., (1968), J. Appl. Mech., 35, pp. 363-371.
[16]
Rice, J.R., and Tracy, D.M., (1969), J. Mechanics Physics Solids, 17, 210-217.
[17]
Hancock, J.W. and Meckenzie, A.C., (1976), J. of Mech. Phy. Sol., 24, pp. 147169
[18]
Thomson, R.D., and Hancock, J.W., (1984), Int. J. Frac. 26, pp.99-112.
[19]
Manjoine, M.J., (1982), Welding Research Supplement, 50s-57s.
[20]
Bonora, N., and Milella, (2000) P.P., USAF/AFRL Contract no F61775-00-
WE029 Final Report.
[21]
Johnson G.R. and Cook W.H., (1983), Proc. 7th Inter. Symp. Ballistic, Am. Def.
Prep. Org. (ADPA), Netherlands.
[22]
Tuler, F.R., and Butcher,B.M., (1968), Int. J. Fract. Mech., 4, pp. 431-437
[23]
Gurson, A.L., J. Engn. Mat. Tech.,(1977), 99, pp. 2-15
[24]
Tvergaard, V. and Needleman, A., (1984), Acta Metall., 32, p.157-169
[25]
Needleman, A., and Rice, J. R., (1978), Limits to ductility by plastic flow
localization in Mechanics of Sheet Metal Forming, Koistinen and Wang Eds.,
Plenum Publ. Company, New York, pp. 237-265.
[26]
Gao, X., Faleskog, J., Shih, C.F., and Dodds, R. H., (1998), Eng. Fracture
Mechanics, 59, N 6, pp. 761-777
[27]
Brocks, W., Klingbeil, D., Kunecke, G., and Sun, D.Z., (1994), in Constraint
Effects in Fracture: Theory and Application, ASTM STP 1244, Kirk and Bakker
Eds., American Society for Testing and Materials, Philadelphia.
[28]
Curran, D. R., Seaman, L., Shockey, D.A., in Shock Waves and High Strain
Rates Phenomena in Metals, Ed. Meyers M.A. and Murr L.E., Plenum, NewYork
[29]
Seaman L., Curran, D.R., and Shockey, D.A., (1976) J. Appl. Phys., 47, p. 4814.
[30]
Lemaitre, J., A Course on Damage Mechanics, Springer-Verlag, Berlin, 1992
[31]
Chandrakanth, S. and Pandey, P.C., (1993), Int. J. Fracture, 60, R73-R76.
[32]
Bonora, N., (1997), Int. J. of Fracture, 88, 359-371
79
6 Hopkinson Bar
6.1 Principio di funzionamento
La barra di Hopkinson (o apparato di Kolsky) è, ad oggi, la tecnica sperimentale più
utilizzata per la caratterizzazione della risposta meccanica dei materiali, in regimi di
velocità di deformazione compresi tra 102 e 104 s-1. Il principio di funzionamento è
basato sull’assunzione che nel provino e nelle barre che costituiscono il sistema di prova
si realizzi uno stato di sforzo uniassiale. Una configurazione tipica dell’apparato di
prova e degli strumenti necessari alla rilevazione dei dati è schematicamente mostrato in
Figura 6.1.
Figura 6.1 - Schematizzazione dell'apparato e della strumentazione di una configurazione classica
della Hopkinson in compressione, 0.
Il provino è fissato tra le due barre incidente e trasmittente. Una terza barra (proiettile o
striker bar), accelerata per mezzo dell’energia trasmessagli da una molla o da una
pistola a gas, colpisce la barra incidente provocando un impulso che viaggerà in essa
fino a raggiungere il provino. All’interfaccia col provino, parte dell’impulso sarà
trasmesso e parte riflesso in rapporto alle impedenze meccaniche della barra e del
provino. L’impulso trasmesso, dopo aver attraversato il provino, all’interfaccia con la
barra trasmittente, sarà in parte riflesso e in parte trasmesso alla barra stessa. Se la
lunghezza dell’impulso è sufficientemente più lunga della lunghezza del provino, le
80
ripetute riflessioni che si realizzano garantiscono, nel provino stesso, che la
deformazione e la velocità di deformazione possano ritenersi uniformi. La lunghezza
dell’impulso è pari al doppio della lunghezza della striker bar. Le barre incidente e
trasmittente devono essere sufficientemente snelle da garantire l’instaurarsi di uno stato
di sforzo quanto più prossimo a quello uniassiale. La sezione delle barre è scelta in
modo tale che, data l’intensità dell’impulso generato, esse abbiano, durante la prova, un
comportamento elastico, mentre il provino, di sezione minore, si deforma plasticamente.
Figura 6.2 - Schematizzazione degli impulsi di deformazione alle interfacce barre provino 0.
Lo stato di sforzo e deformazione che si realizza nel provino durante la prova può essere
ricavato dalla conoscenza dei segnali di deformazione elastica sulle barre incidente e
trasmittente. Nella rappresentazione schematica del provino e delle barre di Figura 6.2,
sono riportati gli impulsi incidente, ε i , riflesso, ε r , e trasmesso, ε t . Indicando con i
pedici 1 e 2 le due estremità del provino, i loro spostamenti possono essere scritti come:
t
u1 = ∫ c0ε1dt
0
(6.1)
t
u2 = ∫ c0ε 2 dt
0
in cui c0 è la velocità dell’onda elastica nelle barre di Hopkinson. Se si scrivono le
equazioni (6.1) in termini di impulsi incidente, riflesso e trasmesso, si ottiene:
t
u1 = c0 ∫ ( ε i − ε r ) dt
0
(6.2)
t
u2 = c0 ∫ ε t dt
0
81
con la usuale convenzione, nella dinamica dell’impatto, di assumere positivi gli sforzi e
le deformazioni di compressione. La deformazione media nel provino è:
εs =
u1 − u2
L
(6.3)
o, in termini d’impulso di deformazione:
εs =
c0
L
t
∫ (ε
i
− ε r − ε t ) dt
(6.4)
0
in cui L è la lunghezza del provino. Le forze alle estremità del provino possono essere
scritte come:
P1 = EA ( ε i + ε r )
(6.5)
P2 = EAε t
in cui E e A indicano rispettivamente il modulo di Young e la sezione delle barre di
Hopkinson. La forza media è pari a:
Pav =
EA
(ε i + ε r + ε t )
2
(6.6)
Se si assume, per l’equilibrio, P1 = P2 ,si ha:
(ε i + ε r ) = ε t
(6.7)
e, quindi, dall’equazioni (6.4):
εs =
c0
L
t
∫ (ε
t
− ε r − ε r − ε t ) dt
(6.8)
0
Per un provino si sezione As , si ottengono la deformazione, lo sforzo e la velocità di
deformazione, come:
t
−2c0
εs =
ε r dt
L ∫0
σs = E
A
εt
As
82
(6.9)
(6.10)
εs =
−2c0
εr
L
(6.11)
6.2 Simulazione numerica della Hopkinson bar
L’assunzione fondamentale della teoria alla base del principio di funzionamento della
barra di Hopkinson è che lo stato di sforzo possa essere assunto uniassiale. Nel presente
lavoro si è verificato che tale ipotesi può essere considerata vera sia per la classica
configurazione a compressione della Split Hopkinson Pressure Bar, sia per una
configurazione alternativa, che permette di effettuare la prova direttamente in trazione.
6.2.1 Prova di compressione
Nelle simulazioni delle prove di compressione un provino di forma cilindrica con
diametro d = 8 mm e lunghezza l = 4 mm è stato sottoposto all’impulso generato da
una striker bar lunga 250 mm . Il diametro, D , delle barre di Hopkinson è di 10 mm e
la loro lunghezza è L = 1100 mm . Il materiale di cui è formato il provino è rame OFHC
le cui proprietà meccaniche sono riportate in Tabella 5.1. Le due barre di pressione e lo
striker sono costituite, invece, di un acciaio maraging per cui si è adottato un modello
elasto-plastico perfetto con carico di snervamento Y pari a 1764 MPa. Uno schema del
sistema di prova è riportato in Figura 6.3.
Figura 6.3 - Schema del sistema di prova della barra di Hopkinson a compressione simulato agli
elementi finiti.
83
Il rapporto l / d = 0.5 del provino ed il coefficiente d’attrito nullo tra le interfacce con le
due barre di pressione che lo trattengono, sono stati scelti allo scopo di riprodurre le
condizioni più favorevoli per minimizzare l’effetto delle inerzie e dei fenomeni
d’attrito. La simulazione numerica è stata effettuata con il codice implicito MSC.Marc:
la mesh che discretizza il campione è formata da elementi, in configurazione
assialsimmetrica, di forma rettangolare di dimensioni pari a 0.5 × 0.25 mm 2 ; l’intervallo
di tempo d’integrazione è stato scelto di ∆t = 10 −7 s che si rivela il giusto compromesso
tra la precisione richiesta e i tempi di calcolo reclamati dal fenomeno, la cui durata è
dell’ordine dei 500 µs . Per la discretizzazione delle barre di pressione sono stati
impiegati elementi rettangolari di dimensioni maggiori, 2 × 1mm 2 . In Figura 6.4 è
riportato il particolare del modello che interessa le interfacce barre provino.
Figura 6.4 - Particolare dell'interfaccia tra barre e provino nella configurazione non deformata
iniziale.
Per modellare il contatto tra i corpi si è fatto ricorso a quattro mesh indipendenti,
ciascuna associata ad un singolo corpo deformabile (il provino, lo striker, e le due barre
di pressione) avendo cura di assegnare una maggiore risoluzione alla mesh che
rappresenta il provino con un duplice scopo:
‰ seguire con precisione il moto di scorrimento sulle barre di pressione e
84
l’eventuale compenetrazione dei corpi a contatto;
‰ monitorare l’evoluzione dei gradienti di tensione indotti.
Le storie temporali delle deformazioni acquisite sulle mezzerie delle barre di Hopkinson
hanno permesso di confrontare i risultati ottenuti numericamente con quanto previsto
dalla teoria. In Figura 6.5 si riportano gli andamenti delle onde di tensione registrate
durante la prova.
Figura 6.5 - Andamenti delle tensioni durante la prova di compressione.
Sono indicate in blu ed in rosso le storie di carico che registrano rispettivamente il nodo
disposto sulla barra incidente (che rileva l’impulso incidente e in seguito quello riflesso)
e quello presente sulla barra trasmittente o d’output (che registra solo l’impulso
trasmesso), per una prova in cui lo striker viaggia inizialmente ad una velocità di
15m / s . L’onda di compressione σ i generata dall’urto con lo striker, propaga lungo la
barra fino a raggiungere l’interfaccia con il provino, dove, è parzialmente riflessa come
onda di trazione σ r e parzialmente trasmessa σ t nella barra d’output. L’impulso di
compressione σ i rilevato dal nodo che funge da estensimetro sulla barra incidente, ha
una forma pressoché rettangolare; tuttavia si distinguono in modo nitido sul suo fronte,
le tipiche oscillazioni di Pochhammer-Chree imputabili agli effetti d’inerzia radiale e
85
che si attenuano molto velocemente per via dell’elevato rapporto di snellezza
l / d = 110 scelto per le barre.
La durata temporale e l’intensità dell’onda di tensione generata dall’urto rispettano
egregiamente le previsioni teoriche secondo cui, per uno stato di sforzo elastico,
unidimensionale si dovrebbe avere:
1
2
1
2
σ = ρ C0V0 = ρ
E
ρ
⋅ V0 ≅ 292.5MPa
(6.12)
t=
2l
=
C0
2l
≅ 102.6 µ s
E ρ
Anche la contemporaneità con cui le due onde riflessa e trasmessa giungono sui due
estensimetri opposti è più che buona. Altrettanto soddisfacente è l’accordo con le
previsioni teoriche per quel che riguarda l’equilibrio delle tensioni nel provino, sebbene
le due onde non abbiano una forma simile a quella rettangolare, è facile verificare che
per ogni istante di tempo risulta vero, con buona approssimazione, che la somma degli
impulsi riflesso e trasmesso uguaglia quello incidente:
σ r (t ) + σ t (t ) = σ i (t )
(6.13)
Proprio la forma dell’onda riflessa contraddistingue in modo univoco la prova di
compressione: se si fa riferimento alla formula (6.11), si evince che gli strain rates
applicati al provino non sono perfettamente costanti, ma che i valori spesso citati nelle
prove sperimentali non sono altro che un valore "mediato" dei corrispettivi andamenti
temporali. A conferma di quanto detto, nella Figura 6.6 sono riprodotti alcuni andamenti
temporali delle velocità di deformazione subite dal provino, per diverse velocità
d’impatto.
I diagrammi sono stati ottenuti effettuando una campagna di prove con velocità
d’impatto crescenti della barra proiettile generando onde di sollecitazione la cui
ampiezza sia sempre inferiore al limite elastico del materiale costituente le barre di
Hopkinson, imponendo che:
σ max =
1
ρC 0V0 ≤ Y
2
86
(6.14)
Figura 6.6 - Andamenti temporali degli strain rates per diverse velocità d’impatto.
Figura 6.7 - Diagrammi tensione-deformazione per le diverse velocità di carico.
Sono due le cause concomitanti che provocano gli andamenti decrescenti delle velocità
di deformazione:
‰ il progressivo incrudimento del materiale;
87
‰ l’aumento della sezione resistente opposta dal campione a causa dell’espansione
radiale per effetto Poisson.
Questo spiega gli andamenti temporali decrescenti e le pendenze maggiori per i casi in
cui si genera nel provino un’onda plastica più intensa. La conferma è fornita dagli
andamenti dei diagrammi tensione-deformazione, riportati in Figura 6.7, ottenuti
secondo la teoria, dall’analisi delle onde di deformazione riflessa e trasmessa registrate
dai trasduttori, attraverso le equazioni (6.9) e (6.10).
In seguito sono state eseguite due prove con la medesima geometria per mettere a
confronto il metodo d’integrazione Newmark β con il single step Houbolt. I risultati
ottenuti hanno dimostrato la sostanziale conformità dei due algoritmi di calcolo, ma a
favore del secondo metodo sembra deporre una maggiore stabilità dimostrata da una
distribuzione delle tensioni più regolare nei primi istanti di tempo dell’analisi. Ad
esempio in Figura 6.8 si è riportato l’andamento delle tensioni all’interno delle barre di
pressione (proiettile ed incidente) dopo un tempo ∆t = 10µs , ottenuto con i due metodi
per la stessa velocità d’impatto V = 25 m / s .
Figura 6.8 - Distribuzione delle tensioni sull’asse di simmetria delle barre di pressione a
10 µ s dall’istante in cui è avvenuto l’impatto.
88
La distribuzione generata dal metodo di Houbolt modificato è certamente più regolare e
dotata di una perfetta simmetria rispetto all’asse d’impatto delle barre, contrariamente a
quanto fatto dal metodo di Newmark β che è sporcato da un rumore "di fondo"
incomprensibile. Come anticipato dagli sviluppatori del codice il primo metodo si
conferma più indicato per le analisi dinamiche che coinvolgono corpi a contatto ed in
cui è necessario attenuare l’effetto dei disturbi d’alta frequenza eccitati dall’urto. Aldilà
delle sottili divergenze iniziali però, i risultati delle analisi numeriche sul medio e lungo
termine si rivelano del tutto equivalenti, come confermato dagli andamenti delle onde di
deformazione riflessa e trasmessa misurate sugli estensimetri mostrate nella Figura 6.9.
Figura 6.9 –Onde di deformazione incidente e riflessa registrate sugli estensimetri dopo circa 320 µs
dall’impatto.
La condizione d’equivalenza dei due metodi d’integrazione numerica sul lungo periodo
è confermata anche dall’uguaglianza delle deformate finali del provino, che in ambo i
casi non hanno presentato, a causa dell’attrito nullo sulle interfacce con le barre, il
tipico profilo curvo a forma di botte. Il "barreling" è, infatti, un fenomeno molto
comune nelle prove di compressione dinamiche in cui la lubrificazione non si dimostra
del tutto efficace come mostrato nei dettagli della Figura 6.10.
89
Figura 6.10 – Confronto tra un provino cilindrico non deformato (a) ed i profili finali deformati con
barreling (b) e senza (c).
6.2.2 Prova di trazione
Per riuscire ad ottenere l’intera curva sforzo deformazioni, fino a rottura, è necessario
effettuare la prova alla barra di Hopkinson in trazione. Nel corso degli anni sono stati
concepiti numerosi dispositivi, alcuni dei quali assai ingegnosi, in grado di trasformare
l’impulso compressivo, generato dall’impatto, in un impulso tensile, [2]. Una
configurazione molto interessante è quella proposta da Staab e Gilat, perché è in grado
di generare direttamente un’onda di trazione, senza ricorrere dell’azione della barra
proiettile. Questo obiettivo si raggiunge accumulando nella parte posteriore della barra
incidente un precarico di deformazione, e lasciandolo in seguito libero di propagare
lungo la barra fino ad investire il provino, Figura 6.11.
Figura 6.11 - Schema funzionale dl dispositivo di Staab e Gilat.
90
Per realizzare un simile sistema di carico, si rendono necessari però, un servosistema
idraulico per porre in trazione la parte della barra preposta a tal fine, ed una morsa
capace di fissarne una sezione di lunghezza prefissata L . Forse è proprio la morsa a
rilascio istantaneo, l’elemento essenziale dell’impianto, e perciò la corretta
progettazione di questo assemblato ha ricadute immediate sulla qualità dell’impulso di
trazione prodotto. Un semplice schema funzionale della morsa, riprodotto nel
particolare di Figura 6.11, prevede l’azione concomitante di due bracci meccanici (con
le estremità fissate sulla loro parte inferiore) a serrare i lati della barra incidente sotto
l’azione di un bullone intagliato. Per far partire l’impulso di trazione, il bullone
sopraindicato è stretto finché non cede di schianto lasciando la barra libera di spostarsi
in senso assiale. In quel preciso istante un impulso di trazione, d’ampiezza pari alla
metà della deformazione elastica accumulata fino a quel momento, inizia a diffondere
nella barra d’input verso il provino. Contemporaneamente, in modo perfettamente
analogo, un’onda di rilascio d’uguale ampiezza inizia a muoversi dalla morsa verso
l’estremità opposta della barra che è vincolata al sistema servo-idraulico di trazione. A
causa dell’impedimento esercitato da questo dispositivo, quando l’onda di rilascio è
riflessa dall’estremità posteriore della barra, si annulla bruscamente l’ampiezza
dell’impulso che va generandosi. Il risultato del transitorio descritto è perciò un impulso
di trazione la cui durata è pari a: t = 2 L C 0 (tempo che esso impiegherebbe per coprire
il doppio della distanza esistente tra la morsa ed il sistema di trazione) che si muove
lungo la barra di pressione in direzione del provino. E’ l’intensità dell’impulso che si
desidera generare a suggerire il tipo di materiale da impiegare per il bullone e la
profondità del relativo intaglio. Comunemente, però, i materiali utilizzati a tale scopo
sono le leghe d’alluminio Al 6061-T6 o Al 2024-T6, perché dotate di una duttilità
minima ma non di un eccessivo livello di fragilità che potrebbe persino impedire di
raggiungere il livello di deformazione richiesto durante il precarico. Questo tipo
d’impianto, che ultimamente sembra essersi confermato il più compatto ed affidabile in
assoluto, ha tratto ampia ispirazione dalle tecniche sperimentali già ampiamente
collaudate per altre tipologie di prove, basti pensare a come tecniche analoghe siano
state per lungo tempo un cardine delle prove dinamiche a torsione.
La prova a trazione con il dispositivo della barra di Hopkinson su provini di rame, è
stata simulata facendo ricorso ad un’unica mesh, formata dall’unione delle barre di
91
pressione e del provino interposto. La continuità del reticolo di calcolo si è resa
necessaria per schematizzare il comportamento del campione che, nella realtà, è
assicurato alle barre per via delle estremità filettate di cui dispone. Il provino ha la
caratteristica forma ad "osso di cane" come per le normali prove di trazione, ma
dimensioni più contenute, con un diametro d = 4 mm , lunghezza utile l = 8 mm , e
raggio di raccordo r = 1 mm ; mentre le barre hanno un diametro D = 9 mm ed una
lunghezza L = 1100 mm . Per rappresentare il comportamento del sistema di carico della
barra incidente, nel modello numerico, si è fatto uso di una speciale condizione di
vincolo denominata "tying". Ai nodi disposti sull’estremità remota della barra incidente
è stato imposto uno spostamento ∆x prefissato mentre i nodi della barra disposti su una
sezione radiale ad una distanza L = 450 mm dall’estremità, sono stati bloccati per
simulare il sistema di tenuta della morsa. In seguito è stato istantaneamente rimosso il
vincolo così imposto e si è lasciato il sistema libero di evolvere. Così facendo si genera
nella barra incidente un’onda di trazione di durata t = 2 L C e di ampiezza pari alla
metà della deformazione accumulata sulla parte "in tiro" della barra fino a qualche
istante prima. Il passo adottato per l’integrazione diretta del transitorio generato, ha
dimensioni variabili ed è pari a ∆t = 2 ⋅ 10 −7 s finché l’onda di trazione non giunge sul
bordo della barra incidente, ma diventa pari a ∆t = 4.5 ⋅ 10 −8 s per i successivi 190 µs
quando l’onda elasto-plastica investe in pieno il provino.
Figura 6.12 – Particolare della mesh adottata per la discretizzazione del provino di rame puro.
92
La necessità di un time step variabile è dettata dalla durata temporale dell’intero
processo e dalla mesh di dimensioni variabili adottata per la geometria in esame. Le
dimensioni minime degli elementi approssimanti il provino sono di 0.20 × 0.25mm 2 per
quelli disposti sulla sua parte rettilinea come mostrato nella Figura 6.12.
Diversamente dalle prove di compressione, in tutte quelle di trazione realizzate, si è
avuto cura di implementare nel codice di calcolo il modello di danneggiamento non
lineare presentato.
Le peculiarità salienti della prova di trazione si sono confermate:
9 gli andamenti degli strain rates sono più regolari rispetto alla prova di
compressione;
9 le problematiche connesse alle "oscillazioni" presenti nelle curve tensione-
deformazione, tipiche di buona parte delle prove di trazione ad elevate velocità
di deformazione, costituiscono un punto critico;
9 la rottura avviene a seguito del fenomeno di strizione.
Figura 6.13 - Andamenti delle onde di tensione registrate durante le prove.
93
Figura 6.14 - Andamenti temporali delle velocità di deformazione raggiunte.
Nelle prove simulate al calcolatore si è proceduto lasciando inalterata la geometria del
sistema ed inviando sul campione onde di sollecitazione d’ampiezza crescente. Dai
profili di deformazione registrati sugli estensimetri è emerso un andamento più regolare
di quanto non risulti nelle prove di compressione. La tipica forma delle onde di trazione
e di compressione registrate sulle barre di pressione è mostrata in Figura 6.13, in cui
appaiono ben visibili gli impulsi di trazione incidenti e gli impulsi di compressione
riflessi e, più deboli, quelli trasmessi.
La forma molto regolare degli impulsi riflessi è un buon indice degli andamenti
temporali delle velocità di deformazione imposte al provino. Facendo ricorso alla
relazione (6.11) si ottengono gli andamenti mostrati in Figura 6.14, contraddistinti da
forme molto simili agli andamenti ideali "a gradino" che sarebbero indicativi di una
prova con velocità di deformazione perfettamente costante.
Ciononostante, anche la prova di trazione presenta alcuni limiti intrinseci. Innanzi tutto,
rimane confermato dall’analisi agli elementi finiti che come suggerito da diversi
ricercatori, 0, è molto difficile ottenere informazioni attendibili sul comportamento del
materiale ad elevati strain rates nel campo elastico. Sono sostanzialmente due le cause
di tale difficoltà:
94
‰ i fenomeni delle oscillazioni della risposta, che crescono con la velocità di
deformazione della prova;
‰ la mancata uniformità dello stato di deformazione interno del campione e degli
strain rates imposti nei primi istanti del test.
Ad aggravare la disomogeneità del campo di deformazione contribuiscono anche le
forme e le dimensioni del provino adottato nella prova di trazione. Tali problematiche
fanno in modo che il modulo di Young del materiale misurato in condizioni dinamiche
si riveli per diversi materiali, minore di quello ottenuto con prove di trazione
quasistatiche.
A conferma di quanto appena detto, in Figura 6.15, sono riportati gli andamenti dei
diagrammi σ − ε ottenuti con diverse velocità di carico per lo stesso materiale.
Figura 6.15 - Curve tensione-deformazione ottenute per diverse velocità di carico.
Nella realtà esiste per le prove realizzate attraverso il dispositivo della barra di
Hopkinson, anche un limite superiore d’applicabilità, determinato dalle deformazioni
plastiche accumulate. La condizione di monodimensionalità dello stato tensionale,
infatti, cessa di esistere non appena il provino inizia a subire una localizzazione del
flusso plastico a causa del fenomeno del necking. Come per gli altri test di trazione
uniassiale, quando il processo di strizione localizza sul provino, non è più possibile
convertire attraverso la semplice teoria proposta, gli spostamenti delle barre di pressione
nella curva dello strain rate subito dal
95
materiale in prova. Durate il processo di
necking la velocità di deformazione cresce localmente ben oltre i valori che si registrano
nelle condizioni di deformazione uniforme. Nella Figura 6.16 sono riportati gli
andamenti delle velocità di deformazione registrate per una prova in cui il provino è
stato deformato fin quasi a rottura ( D ≅ 0.6 ). Come si può facilmente constatare, le
misure effettuate "localmente" con degli estensimetri virtuali di lunghezza decrescente,
confermano che la velocità di deformazione cresce notevolmente a causa del necking
Figura 6.16 - Misure degli strain rates effettuate con estensimetri di lunghezza variabile.
Figura 6.17 – Distribuzione del danno sulla deformata finale del provino dopo la rottura.
96
Quando l’effetto combinato della triassialità dello sforzo, della temperatura crescente
(causata dal riscaldamento adiabatico del provino) e del danno accumulato con la
deformazione plastica, sono spinte all’eccesso si provoca il cedimento del provino che
presenta la rottura coppa-cono tipica dei materiali duttili mostrata nella Figura 6.17.
6.2.2.1 Prove di trazione con ARMCO-iron
Il campo d’applicazione del dispositivo di Kolsky per le prove di trazione può essere
esteso ben oltre l’inizio del necking facendo ricorso alle capacità della fotografia ad alta
velocità. Sono attualmente disponibili, infatti, cineprese a tamburo rotante capaci di
acquisire fino a 200000 frames in un secondo con tempi d’esposizione inferiori ai 4 µs .
Di questa tendenza affermatasi negli ultimi anni, ci si è avvalsi di recente, per effettuare
un confronto approfondito tra dati sperimentali e risultati di natura numerica. Nella
Figura 6.18 ad esempio è riportata una tipica sequenza fotografica realizzata durante
una prova di trazione con la barra di Hopkinson ad intervalli di 10 µs .
Per esaminare la validità dei modelli costitutivi adottati nelle simulazioni assistite da
calcolatore si può effettuare un confronto tra le deformate reali e numeriche accumulate
durante alcuni test come il cilindro di Taylor o i proiettili forgiati tramite esplosivi. Di
recente Noble et al., [3], hanno evidenziato che per le prove di trazione, la riduzione
dell’area nella zona di strizione e gli incrementi di temperatura ivi registrati, sono due
parametri di confronto assai più sensibili. Soprattutto le previsioni riguardanti la
temperatura, sono molto importanti quando si fa ricorso a modelli costitutivi avanzati
che implementano una dipendenza esplicita da questo parametro. Dal momento che il
legame costitutivo del materiale e la legge d’evoluzione del danno, integrate nel calcolo,
sono contraddistinti da una forma spiccatamente non lineare, si è cercata una conferma
dei parametri numerici utilizzati nelle simulazioni attraverso un confronto con alcune
prove sperimentali rinvenute in letteratura. Sono così state svolte, prove di trazione su
provini di ARMCO-iron con diametro d = 3 mm e lunghezza utile l = 8 mm lasciando
inalterata la geometria dell’intero impianto, ed utilizzando per gli elementi formanti la
mesh dimensioni minime di 0.25 × 0.25 mm 2 (come visualizzato in Figura 6.19) con
tempi inferiori d’integrazione pari a ∆t = 4 ⋅ 10 −8 s .
97
Figura 6.18 - Sequenza fotografica per una prova di trazione.
Le costanti adottate per il legame costitutivo di Johnson e Cook e per il modello di
danno di Bonora, per il materiale in esame, sono elencate rispettivamente nelle Tabella
5.1 e Tabella 5.4.
La prima serie di prove ha puntato a verificare la fondatezza del valore dello
smorzamento numerico γ i scelto per le analisi numeriche effettuate. Sono stati
confrontati, a tal fine, i profili teorici e numerici assunti dall’onda di trazione che
propaga nella barra incidente. Le prove numeriche sono state svolte adottando nei tre
98
casi valori di γ i pari rispettivamente a: 0.4, 0.8 ed 1.2 come visualizzato nella Figura
6.20.
Figura 6.19 - Mesh adottata nelle prove con ARMCO-iron.
Figura 6.20 – Confronto tra il profilo teorico e quelli numerici dell’onda di trazione.
Gli andamenti dell’impulso incidente risultano in ogni caso abbastanza vicini
all’impulso ideale a gradino che dovrebbe assumere un’ampiezza ε inc = 3.75 ⋅10 −3 ed
una durata ∆t = 185 µs . Ciò nonostante il giusto compromesso, tanto per le oscillazioni
dovute alle inerzie radiali, quanto per il livello di deformazione raggiunto, è garantito
99
dal valore di γ i = 0.8 .
In seguito si è proceduto con un confronto tra risultati numerici e valori sperimentali
registrati in tre prove distinte su provini dalla geometria identica. Sono stati così
confrontati gli andamenti temporali della riduzione percentuale dell’area del provino
causata inizialmente dalla deformazione uniforme e dopo dal necking. Com’è evidente
il trend della simulazione segue in modo soddisfacente gli andamenti delle misurazioni
sperimentali per le quali sono diagrammati il valore minimo, medio e massimo delle tre
prove, Figura 6.21.
Figura 6.21 - Confronto tra dati numerici e sperimentali della strizione.
D’altronde uno scostamento del 3 o 4% sui valori massimi della strizione è
oggettivamente una soddisfacente conferma della validità del modello di danno non
lineare anche per le situazioni (come quella in esame) che coinvolgono evoluzioni
temporali dello stato di triassialità dello sforzo. Non bisogna trascurare che prima di
giungere alla rottura del provino lo stato tensionale interno per questo ultimo, muta
notevolmente, passando da una sollecitazione pressoché monoassiale ad una
marcatamente triassiale. La deformata finale del provino, un istante dopo la rottura, è
riportata in Figura 6.22, da essa sono stati rimossi gli elementi della mesh che hanno
raggiunto il limite critico del parametro di danno. Come si nota chiaramente
100
l’asimmetria di carico, imposta nella prova, causa la rottura del campione non in
corrispondenza della sua linea di simmetria ma su una sezione radiale disposta ad una
distanza dalla prima di circa 1 mm
Figura 6.22 - Deformata finale del provino di Armco iron un istante dopo la rottura.
Figura 6.23 - Confronto tra valori numerici e sperimentali degli strain rates.
101
Questo particolare conferma il buon esito del calcolo, riproducendo l’essenza del
fenomeno reale, così come fatto anche dagli andamenti delle velocità di deformazione
sperimentate dal campione e riportate in Figura 6.23. In essa sono stati confrontati gli
andamenti sperimentali ottenuti tenendo sotto controllo due punti del provino per mezzo
della fotografia ad alta velocità, e gli analoghi andamenti temporali registrati nella
simulazione numerica.
L’altro parametro d’interesse per le simulazioni termo-meccaniche eseguite è la
temperatura. È già stato precisato che il lavoro plastico di deformazione compiuto sui
materiali di natura metallica è dissipato sotto forma di calore. Per i carichi dinamici
molto veloci s’instaurano, molto spesso, le condizioni di adiabaticità del fenomeno
perché il calore è dissipato con una velocità notevolmente inferiore a quella con cui è
generato localmente. In questi casi un’analisi semplificata, suggerisce di imporre che
un’aliquota costante k del lavoro di deformazione è convertita in calore causando
l’innalzamento della temperatura del materiale. Utilizzando un coefficiente di
conversione k = 0.95 , avvalorato dagli studi effettuati in tale ambito negli ultimi anni,
si è monitorata l’evoluzione temporale della temperatura nel punto più deformato del
campione. Come confermato dalla simulazione numerica e dalla Figura 6.24, il provino
subisce un notevole aumento della temperatura fino oltre i 250°C ma, soprattutto, la
velocità di crescita di questa ultima cambia repentinamente quando localmente si
instaura il fenomeno del necking, aumentando di circa ∆T = 180°C in meno di
∆t = 85µs .
Il notevole incremento della temperatura è stato confermato nelle prove sperimentali per
mezzo di una camera di scansione termica ad infrarossi. Con l’ausilio di questo
dispositivo si è riusciti ad ottenere, durante le prove sperimentali, la distribuzione della
temperatura sulla superficie del provino qualche istante dopo la rottura. Il confronto tra
valori sperimentali e numerici è presentato in Figura 6.25 e non deve colpire l’apparente
disuniformità dei dati, causata in realtà dai limiti del dispositivo di misura della
temperatura, incapace di acquisire immagini in un intervallo di tempo ∆t < 400µs . La
realtà è che la mappatura della temperatura sperimentale è stata eseguita circa 2 ms
dopo l’avvenuta rottura del provino, un intervallo di tempo cospicuo se si confronta con
i tempi di rottura di circa 180 µs . Questo ritardo può aver causato la ridistribuzione
della temperatura nelle immediate vicinanze della zona di frattura, che per la
102
simulazione numerica è risultata la più sollecitata come ampiamente anticipato dalla
teoria.
Figura 6.24 - Variazione della temperatura nel punto più sollecitato del provino.
Figura 6.25 - Distribuzioni della temperatura dopo la rottura del provino.
103
6.3 Conclusioni
Le simulazioni numeriche effettuate hanno permesso di verificare che nella prova alla
barra di Hopkinson, lo stato di sforzo che viene a realizzarsi sia effettivamente molto
prossimo ad uno stato unidimensionale.
Si è verificato, inoltre, che con un’attenta progettazione della prova, nello specifico
geometria del provino in relazione alla geometria delle barre e alle impedenze
meccaniche dei materiali utilizzati, la deformazione e la velocità di deformazione
possono essere ritenute, con buona approssimazione, uniformi all’interno del provino.
Questo è di importanza rilevante, in quanto, permette di accettare i risultati ottenuti con
questa tecnica sperimentale come identificativi del comportamento meccanico del
materiale in regime dinamico.
Si è riscontrato che i limiti maggiori in tale tecnica sperimentale sono dati dalla
difficoltà di mantenere costante, durante l’intera durata della prova, la velocità di
deformazione. In particolare si è osservato come nella prova di compressione la velocità
di deformazione decresce a causa dell’aumento di sezione per effetto poisson. Nella
prova a trazione, al contrario, è proprio l’effetto di strizione a consentire il
mantenimento di un’elevata velocità di deformazione ad un valore pressoché costante.
Dall’esperienza acquisita con le analisi numeriche effettuate si è partiti per la
progettazione e di una barra di Hopkinson a trazione in fase di realizzazione. Le barre
sono state dimensionate in modo da permettere la caratterizzazione del comportamento
meccanico in regime dinamico di materiali metallici quali il rame, gli acciai, alcune
leghe di nichel, come il waspaloy etc. Particolare attenzione ha richiesto la
progettazione del sistema di afferraggio. Al buon funzionamento di tale sistema, infatti,
è condizionata la limpidezza dell’impulso di trazione generato e, quindi, la pulizia della
prova effettuata.
104
Bibliografia
[1]
Zukas, J. A, Nicholas, T., Swift, H. F., Greszczuc, L. B. and Curran, D. R.,
Impact Dynamic, John Wiley & Sons, New York, 1992.
[2]
Zukas, J. A., High Velocity Impact Dynamics, John Wiley & Sons, New York,
1990.
[3]
Noble, J.P., Goldthorpe, B.D., Church, P. e Harding, J., "The use of the
Hopkinson bar to validate constitutive relations at high rates of strain", Journal
of the Mechanics and Physics of Solids, 47, pp. 1187-1206, 1999.
105
7 Flyer Plate Impact Test
L’esperimento del Flyer Plate Impact Test consiste nel realizzare un impatto planare, a
velocità nota, tra due dischi sottili. Un rapporto diametro su spessore elevato (D/h>10)
garantisce uno stato di deformazione uniassiale in prossimità dell’asse di simmetria dei
dischi. Questa configurazione sperimentale rappresenta una delle poche configurazioni
geometriche per le quali la trattazione teorica, come descritta al paragrafo 2.2.5, è
disponibile in forma esatta e può essere utilizzata per la verifica ed il confronto con i
risultati numerici. Come è già stato introdotto, è possibile, anche per impatti iperveloci,
generare un’onda di shock, perciò il Flyer Plate Impact Test è generalmente utilizzato
per determinare la curva di Hugoniot del materiale. Si ricorda che tale curva non è
percorsa durante il processo di caricamento, ma rappresenta il luogo dei punti di
equilibrio raggiunti per diverse condizioni della pressione d’impatto; in altre parole,
ogni esperimento permette di determinare un solo punto sulla curva.
Questo test è largamente utilizzato anche per il particolare tipo di rottura che è in grado
di produrre nel disco bersaglio. Tale rottura, denominata spalling, avviene per una
trazione localizzata provocata dalla sovrapposizione dell’onda di compressione, riflessa
sulla superficie libera del target, e della sopraggiungente onda di rilascio.
Nell’esperimento, la misura è effettuata mediante la rilevazione, ad esempio attraverso
tecniche d’interferometria laser, del profilo di velocità di un punto situato sulla
superficie posteriore del disco bersaglio, Figura 7.1 b. La lettura del profilo permette di
ricavare tutte le informazioni necessarie ad identificare il comportamento meccanico del
materiale. La comprensione dei tratti caratteristici del profilo di velocità, può essere
agevolata dal diagramma lagrangiano di Figura 7.1 a, in cui sull’asse delle ascisse è
riporta la distanza lungo gli spessori dei due dischi e sulle ordinate il tempo. Tale
diagramma permette di visualizzare, con tratti di retta, i fenomeni di propagazione e
riflessione delle onde durante il processo d’impatto; la pendenza del tratto di retta che
rappresenta una data onda ne indica la sua velocità di propagazione.
Al momento dell’impatto, t = 0 , le onde elastiche, generate alla superficie di contatto,
iniziano a propagare in direzione delle superfici posteriori dei due dischi. Se il materiale
che costituisce i due dischi è lo stesso, l’impatto è denominato simmetrico: allora, le
onde propagheranno alla stessa velocità e
saranno rappresentate da tratti di retta
106
simmetrici rispetto alla superficie d’impatto. Nel caso in cui il limite elastico di
Hugoniot sia superato, saranno generate due onde plastiche, che propagheranno alla
velocità determinata dall’equazione (2.4), con gli stessi meccanismi descritti per le onde
elastiche.
Figura 7.1 – a) Diagramma lagrangiano caratteristico di un impatto planare simmetrico; b) Tipico
profilo di velocità rilevato in un Flyer Plate Impact Test.
Lp
σx
Ce
Cp
Cp
Ce
Distanza
lo spessore
Distance lungo
from impact
plane
Figura 7.2 - Onda di stress generata in un Flyer Plate Impact Test a velocità moderata.
Le due onde, elastica e plastica, che viaggiano nel proiettile raggiungono la sua
superficie libera e vengono quindi riflesse come onde di trazione. Quando queste
raggiungono la superficie d’impatto, t = tc , entrano nel disco proiettile come onde di
rilascio e i due dischi, che fino a questo momento viaggiavano uniti, si separano.
A questo punto nel disco impattato si sta propagando il caratteristico impulso a gradino,
schematizzato in Figura 7.2, la cui forma è dovuta alla differenza di velocità di
107
propagazione degli impulsi elastici e plastici. La lunghezza dell’impulso, L p , è pari al
doppio dello spessore del disco proiettile.
Quando l’impulso raggiunge la superficie posteriore del disco bersaglio, questa viene
accelerata e il fenomeno può essere seguito dalla lettura del profilo di velocità. Con
l’arrivo delle onde di compressione, elastica e plastica, la velocità della particella sale,
dapprima, fino al valore che compete al limite elastico di Hugoniot, punto “A”, e poi al
valore massimo del plateau orizzontale. L’arrivo dell’onda elastica di rilascio abbassa la
velocità fino al valore corrispondente al punto “D”, mentre il processo di scaricamento è
completato dall’arrivo dell’onda di rilascio plastica, curva tratteggiata in Figura 7.1.
L’onda riflessa dalla superficie posteriore del disco bersaglio si sovrappone alla
sopraggiungente onda di rilascio, su un piano che ritrova ad una distanza, dalla
superficie libera, pari allo spessore del disco proiettile. Tale sovrapposizione genera un
impulso di trazione che, se sufficientemente elevato, provoca la rottura per spalling. In
tal caso, l’onda generata dalla separazione delle superfici di rottura, una volta raggiunta
la superficie libera, provoca la risalita della velocità e il caratteristico segnale
denominato “spall signal”.Nel caso in cui la velocità d’impatto sia tale da generare
un’onda d’urto, il suo profilo potrebbe essere schematicamente descritto dalla Figura
7.3, e sul profilo di velocità non sarà più presente lo scalino dovuto al precursore
elastico.
σx
C
Distance from impact plane
Figura 7.3 - Profilo di un'onda d'urto.
108
7.1 Simulazione numerica del Flyer Plate Impact Test
La configurazione del Flyer Plate Impact Test, come detto in precedenza, permette di
realizzare uno stato di deformazione uniassiale in prossimità dell’asse di simmetria dei
dischi. Tale condizione rende superflua, in questa prima fase del lavoro, la
discretizzazione dei due dischi nella loro interezza e, quindi, la geometria da modellare
può essere ridotta ad una semplice striscia di elementi in deformazione piana, in cui gli
spostamenti verticali siano impediti. Tale modello è indicato col termine “single strip
model”. Le dimensioni di ciascun elemento, scelte coerentemente all’RVE legato al
modello di danno utilizzato, sono di 0,1× 0,1mm 2 . L’uso del modello di danno permette,
attraverso la tecnica dell’element removal, la creazione della superficie di rottura una
volta realizzate le condizioni di spall. In condizioni d’impatto planare, l’evoluzione del
danno con la deformazione plastica è estremamente limitata, a causa della riduzione
duttilità, a valori prossimi alla deformazione di soglia, per l’elevata triassialità dello
stato di sforzo. In questo caso, il modello CDM è simile ad un criterio di rottura
improvvisa, ma la rottura è il risultato dell’accoppiamento geometria e materiale e non
richiede procedure di calibrazione post test.
La prima simulazione volta alla verifica delle capacità di previsione del modello
numerico realizzato, ha riguardato l’impatto simmetrico di due dischi di rame OFHC
secondo la configurazione riportata in [1], per la quale sono a disposizione i risultati
sperimentali. Le proprietà meccaniche del materiale sono riportate in Tabella 5.1, lo
spessore del disco proiettile è di 2mm , quello del disco bersaglio di 9mm , la velocità
d’impatto è di 185 m s . Nonostante i parametri di danno per il rame siano
sufficientemente noti, Tabella 5.4, il valore della deformazione di soglia è stato
calibrato sul tempo di risalita dello spall signal, ε th = 0, 01 .
La discordanza potrebbe essere imputata al fatto che il primo impulso di compressione,
sebbene non possa generare danno in senso stretto, in virtù del fatto che lo stato di
sforzo è compressivo, potrebbe causare delle modificazioni microstrutturali, quali, ad
esempio la rottura delle inclusioni, in grado di abbattere il valore della deformazione di
soglia. Tale speculazione andrebbe verificata con una campagna sperimentale ad hoc.
Le analisi sono state effettuate con entrambi i codici di calcolo: MSC:Marc e Autodyn.
Attenzione particolare è stata dedicata agli effetti dello smorzamento numerico sui
109
risultati delle simulazioni. Per il codice in formulazione implicita, a seguito di un’analisi
parametrica si sono scelti per i valori dei coefficienti dell’equazione (4.15), α = 0, 0 ,
β = 0, 0 e γ = 0, 4 . Nelle simulazioni effettuate con il codice esplicito si è trovato che i
valori suggeriti dei coefficienti lineare e quadratico dell’equazione (4.21), sono troppo
elevati. In Figura 7.4 sono riportati a confronto i profili di velocità calcolati, con
entrambi i codici, con i coefficienti di smorzamento numerico scelti e con i coefficienti
consigliati per Autodyn.
Figura 7.4 – Effetto dello smorzamento numerico sui risultati delle simulazioni numeriche.
Figura 7.5 - Confronto tra il profilo di velocità ottenuto numericamente ed i dati sperimentali.
110
Infine, il confronto tra il profilo di velocità calcolato e i dati sperimentali, Figura 7.5,
permette di verificare la capacità del modello di riprodurre tutte le caratteristiche chiave
dell’esperimento quali: i tempi di arrivo del precursore elastico e dell’onda plastica sulla
superficie libera del bersaglio; l’intensità e la durata del plateau di massima velocità; i
tempi di arrivo delle onde di rilascio, elastica e plastica.
Figura 7.6 – Evoluzione nel tempo della distribuzione della triassialità dello stato di sforzo, lungo lo
spessore del provino.
Spessore del piano
di spall 150 µm
Figura 7.7 - Evoluzione nel tempo del danno lungo o spessore del disco bersaglio.
In Figura 7.6, è riportata la distribuzione della triassialità dello stato di sforzo a cavallo
111
della superficie di spalling, per diversi istanti di tempo durante il processo di frattura.
Nonostante si abbiano alti valori della triassialità per un tratto considerevole dello
spessore del disco, la rottura si ha solo per la porzione di materiale interessata dal valore
di picco della stessa.
In Figura 7.7, la corrispondente evoluzione del danno è data dall’innesco del
danneggiamento all’avvenuta rottura. L’ area tratteggiata rappresenta la distribuzione
del danno al termine del processo di rimozione degli elementi. Una conferma dell’entità
delle dimensioni della regione danneggiata è data dai risultati sperimentali riportati da
Christy et al., [2], in Figura 7.8.
Figura 7.8 – Distribuzione della porosità nel rame per impatti a diverse pressioni, [2].
112
7.2 Analisi dello “Spall Signal”
Lo spall signal, come definito in Figura 7.1 b, è la parte del profilo di velocità, che ha
inizio quando l’onda di spall raggiunge la superficie di libera del disco bersaglio. Nelle
simulazioni numeriche effettuate e in quelle riportate in letteratura, per le quali sono
stati utilizzati modelli di danno differenti, si trova che il secondo picco di velocità e la
pendenza iniziale dello spall signal sono sempre superiori a quelli rilevati
sperimentalmente.
Figura 7.9 - Confronto tra gli Spall Signals calcolato e misurato per il rame OFHC.
In Figura 7.9 è riportato il particolare del confronto tra gli spall signals, per la
configurazione già analizzata in Figura 7.5. Si prende come tempo di riferimento
iniziale, il momento in cui la prima onda di spall arriva sulla superficie libera e inverte
l’andamento del profilo di velocità. Poiché, in questo diagramma, la pendenza della
curva è una misura dell’accelerazione del punto materiale, la minore pendenza della
curva sperimentale indica una perdita di quantità di moto che, per forza di cose, deve
essere imputata a fenomeni irreversibili. Tali fenomeni devono aver luogo in una fase
successiva al processo di rottura e deve essere legato al meccanismo di separazione
delle superfici di spall.
Da un punto di vista fisico, tale separazione deve essere dipendente dalla microstruttura
del materiale e dal proprio modo caratteristico di frattura. Ad esempio, sebbene sia nel
113
rame sia nell’alluminio, il danno si sviluppa, con la deformazione plastica, con la
nucleazione e la crescita di microvuoti, il processo di coalescenza può essere, nei due
materiali, considerevolmente diverso. Come schematicamente illustrato in Figura 7.10,
per l’alluminio puro, la completa separazione è dovuta alla coalescenza per “voids
sheeting” che, essendo un meccanismo sostanzialmente fragile, richiede una bassa
energia di deformazione;. per il rame, la coalescenza avviene per necking dei legamenti
tra i vuoti, attraverso un meccanismo duttile che richiede un notevole ammontare di
energia.
Ductile failure with “ductile” (high
Rottura duttile dei legamenti
strain energy) intervoid ligament
tra i vuoti - elevata energia di
rupture
deformazione
Ductile
withdei
“brittle”
(low
Rotturafailure
fragile
legamenti
strain
energy)
intervoid
ligament
tra i vuoti - bassa energia di
rupture
deformazione
Figura 7.10 – Differenti meccanismi di coalescenza dei microvuoti nella rottura duttile.
Poiché il processo di formazione dei piani di spall è analogo al processo di formazione
di una cricca duttile, utilizzando gli strumenti della meccanica della frattura, è possibile
quantificare il lavoro necessario alla generazione delle due superfici libere.
L’energia necessaria alla generazione di due superfici libere è pari a:
G = 2Γ
(7.1)
in cui Γ è l’energia libera di superficie che comprende i contributi elastico e plastico,
mentre G è il rateo di rilascio di energia di deformazione che può essere correlata al
valore della tenacità del materiale, K IC . Anche se questo è un concetto puramente
lineare elastico, può ancora essere considerato un valore di riferimento per il caso in
114
esame, in quanto la deformazione plastica lungo il piano di spall è estremamente
contenuta. Di conseguenza:
Γ=
2
1 K Ic
α
2 E
(7.2)
Assumendo la tenacità a frattura del rame pari 60MPa m , ricordando che per uno
stato di deformazione piana è α = ( 1 − ν 2 ) , si ottiene Γ 2000 J m 2 .
Poiché, nella simulazione numerica, non è stato tenuto in conto il meccanismo di
separazione descritto, l’energia dissipata, per unità di superficie, può essere ricavata
dalla differenza tra il segnale di spall calcolato e quello misurato, attraverso la relazione:
Lp ⎞
⎛
1
∆W
J
2
Lp ⎜⎜ 1 − ⎟⎟⎟ = 3842 2
= ρ∆veff
⎝
2
Lb ⎠
∆S
m
(7.3)
in cui è la densità del materiale, Lp e Lb sono gli spessori rispettivamente del disco
proiettile e di quello bersaglio e
∆veff =
1
T
T
∫0
2
[ v fem (t ) − vexp (t ) ] dt
(7.4)
Dividendo il risultato dell’equazione (7.3), per le due superfici, si ottiene
Γ 1921J m 2 , che è in ottimo accordo con il valore stimato con gli strumenti della
meccanica della frattura.
7.2.1 Modello numerico
La verifica, che il meccanismo di dissipazione descritto possa, potenzialmente, essere
responsabile della differenza tra gli spall signals calcolato e misurato, è stata effettuata
implementando, nella simulazione agli elementi finiti, un sistema costituito da molle
non lineari a cavallo del piano di spall. La rigidezza del sistema di molle è diminuita
progressivamente all’aumentare della distanza di separazione. Al raggiungimento di
un’apertura critica u0 , la forza fittizia viene annullata. La legge forza-spostamento
scelta è di forma simile a quella che descrive il legame dei piani cristallini:
⎛ ⎛ u ⎞α ⎞
f = K sin ⎜⎜⎜ π ⎜⎜ ⎟⎟ ⎟⎟⎟
⎝ ⎝ u0 ⎠ ⎠
115
(7.5)
in cui K è l’ampiezza, u0 l’apertura critica, α un esponente di forma, il cui effetto è
illustrato in Figura 7.11. I valori di K e u0 sono stati scelti imponendo che l’area sotto
la curva nel diagramma forza spostamento sia uguale al lavoro dissipato durante il
processo di separazione.
Figura 7.11 - Effetto del coefficiente di forma α sulla risposta del sistema di molle non lineare.
Figura 7.12 - Profilo di velocità calcolato con l’impiego del sistema di molle non lineare a confronto
con i risultati sperimentali.
La Figura 7.12 riporta il profilo calcolato con il modello numerico descritto a confronto
116
con il risultato sperimentale, è importante sottolineare come l’azione delle molle non
lineari nella prima fase del processo di separazione influenzi l’evoluzione dell’intero
segnale di spall.
7.3 Effetti geometrici sul processo di frattura per spalling
La condizione necessaria affinché durante l’impatto si realizzi uno stato di
deformazione uniassiale, come precisato in precedenza, è che nella regione d’interesse,
per l’intera durata del processo, non si senta l’influenza degli effetti di bordo. Anche nel
caso di dischi sottili, la condizione di deformazione uniassiale si realizza
esclusivamente in prossimità dell’asse di assialsimmetria in quanto l’onda di
deformazione radiale, che si genera al bordo libero, necessita, per raggiungere l’asse, di
un tempo maggiore a quello richiesto dall’intero fenomeno. Al diminuire del rapporto
diametro/spessore gli effetti di deformazione radiale possono intervenire direttamente
sulle modalità e localizzazione del processo di rottura per spall. A questo proposito sono
state analizzate diverse configurazioni di Flyer Plate Impact in cui il diametro del flyer è
stato progressivamente ridotto mantenendo inalterate le altre dimensioni e la velocità di
impatto.
Le simulazioni numeriche sono state effettuate utilizzando entrambi i codici numerici
presentati con lo scopo di valutare eventuali effetti dovuti alle diverse formulazioni.
Per uno spessore del flyer di 2mm sono stati esaminati i casi con un diametro di 32 ,
16 , 8 , 4 , e 2mm rispettivamente, come illustrato in Figura 7.13.
Nelle analisi agli elementi finiti è stato utilizzato un elemento a quattro nodi in
formulazione assialsimmetrica con altrettanti punti di gauss cercando di mantenere, per
quanto possibile, costante il livello di discretizzazione del modello al fine di evitare
possibili effetti di mesh. Entrambi i corpi sono considerati deformabili nel contatto. Il
criterio di rottura utilizzato nella simulazione con MSC/MARC fa riferimento al
modello di danno non lineare precedentemente descritto. Nelle simulazioni effettuate
con il codice lagrangiano AUTODYN si è adottato un criterio di rottura basato sulla
pressione massima il cui valore è stato stimato dalle prove effettuate, utilizzando il
modello di danno non lineare, con il codice implicito.
117
D/h=16
D/h=8
D/h=4
D/h=2
D/h=1
V=185 m/s
t=9 mm
Figura 7.13 - Schema riassuntivo delle configurazioni geometriche esaminate.
a)
b)
c)
d)
Figura 7.14 - Profili di velocità calcolati numericamente per le diverse configurazioni geometriche e
con velocità d’impatto di 185m/s: a) D/h=16; b) D/h=8; c) D/h=4; d) D/h=2.
118
In Figura 7.14 a-d sono riportati i profili di velocità rilevati sulla superficie posteriore
del target per le diverse configurazioni. Un valore del rapporto D/h pari a 16 è ancora in
grado di garantire che lo stato di deformazione sia, per lo meno sull’asse di simmetria,
unidimensionale. Per valori più piccoli del rapporto D/h, Figura 7.14 b, si verifica
ancora una rottura per spall, come si può rilevare dalla risalita del segnale di velocità,
anche se gli effetti associati alla deformazione radiale iniziano ad influenzare il processo
di propagazione delle onde lungo l’asse di simmetria riducendo la durata del plateau di
velocità. Per un’ulteriore diminuzione del rapporto D/h la propagazione dell’onda di
sforzo diviene, a causa dell’influenza degli effetti di bordo, sempre più complessa. Nei
profili di velocità riportati in Figura 7.14 c e d, hanno ormai perso ogni attinenza con le
soluzioni di riferimento precedentemente illustrate. L’unica caratteristica ancora
evidente è la discontinuità a cui corrisponde il limite elastico di Hugoniot. In queste
condizioni non è più possibile stabilire sulla base della sola analisi del segnale di
velocità la presenza o meno di cedimento per spall.
a)
b)
Figura 7.15 - a) deformata e stato del materiale ottenuti con AUTODYN per velocità d’impatto di
185m/s e D/h=16; b) deformata e mappa del danno ottenuti con MSC/MARC per velocità d’impatto
di 185m/s e D/h=16.
Nelle Figura 7.15 a e b vengono riportate le deformate e le mappe di danno ottenute
rispettivamente con AUTODYN e con MSC/MARC per due configurazioni simili,
caratterizzate da un elevato rapporto D/h. Per entrambe le simulazioni si ritrova il
119
cedimento per spall caratteristico di un impatto planare. Le distribuzioni di
danneggiamento ottenute con i due criteri adottati confermano la stretta correlazione
esistente tra variabile di danno e pressione idrostatica nelle condizioni di stato di
deformazione uniassiale.
Al diminuire del rapporto D/h la rottura interessa superfici del bersaglio sempre più
piccole. Si è osservato che l’innesco dei processi di rottura per spalling, nel caso di
configurazioni di diverso diametro del target e del flyer, non avviene mai in
corrispondenza dell’asse di simmetria dove invece è atteso dalla teoria. La rottura ha
luogo ad una distanza da tale asse che risulta essere in stretta correlazione con le
dimensioni del diametro del flyer, cosi come la localizzazione del piano di spall è legata
allo spessore dello stesso.
Il motivo per cui la frattura per spalling inizia fuori dall’asse di simmetria può essere
trovato nella differente forma dell’onda di compressione in corrispondenza del bordo
libero del disco proiettile, rispetto a quella sull’asse di simmetria, Figura 7.16.
Figura 7.16 - Impulso di compressione in due differenti posizioni lungo il raggio del disco bersaglio:
sull’asse di simmetria in blu e in corrispondenza del bordo libero del proiettile in nero.
Il bordo libero del proiettile è la superficie su cui l’onda di compressione, generata
nell’impatto, è immediatamente riflessa come onda di trazione. Questa, entrando come
onda di rilascio nel disco bersaglio, scarica parzialmente l’onda di compressione
modificandone il profilo. Tale profilo, di forma triangolare, riflesso dalla superficie
120
libera raggiunge più rapidamente, rispetto all’onda quadra, la condizione di massimo
sforzo di trazione. Se sufficientemente severo, l’impulso tensile porta alla rottura per
spalling il materiale prima di quanto non faccia la corrispondente onda sull’asse di
simmetria.
Una volta innescato, il processo di rottura si propaga radialmente fino ad interessare
l’asse di simmetria.
Questo fenomeno stabilisce le condizioni per la massima estensione radiale della
superficie interessata dal processo di rottura. Dal punto di vista quantitativo, si è
osservato che la localizzazione del primo innesco dei processi di rottura può essere
stimato attraverso la pendenza di una retta ideale, tracciata a partire dallo spigolo del
flyer ed incidente il piano di spall, Figura 7.17. Dalle simulazioni effettuate si è
verificato la costanza del valore di questo angolo per tutte le configurazioni esaminate.
Boundary-effect
line
Spall plane
θ
First spall
symmetry
axis
Figura 7.17 - Schema geometrico della localizzazione dell’innesco del processo di spall.
Inoltre, sempre sulla base di questo criterio, si è osservato che lo spall è impedito per
quelle configurazioni geometriche in cui l’intersezione della retta indicata, con il
termine di boudary-effect line, con la retta del piano di spall avvenga al disotto del
piano di simmetria assiale, come nel caso delle configurazioni riportate dalla Figura
7.19 alla Figura 7.21.
121
a)
b)
Figura 7.18 - a) deformata e stato del materiale ottenuti con AUTODYN per velocità d’impatto di
185m/s e D/h=8; b) deformata e mappa del danno ottenuti con MSC/MARC per velocità d’impatto
di 185m/s e D/h=8.
a)
b)
Figura 7.19 - a) deformata e stato del materiale ottenuti con AUTODYN per velocità d’impatto di
185m/s e D/h=4; b) deformata e mappa del danno ottenuti con MSC/MARC per velocità d’impatto
di 185m/s e D/h=4.
a)
b)
Figura 7.20 - a) deformata e stato del materiale ottenuti con AUTODYN per velocità d’impatto di
185m/s e D/h=2; b) deformata e mappa del danno ottenuti con MSC/MARC per velocità d’impatto
di 185m/s e D/h=2.
122
a)
b)
Figura 7.21 - a) deformata e stato del materiale ottenuti con AUTODYN per velocità d’impatto di
185m/s e D/h=1; b) deformata e mappa del danno ottenuti con MSC/MARC per velocità d’impatto
di 185m/s e D/h=1.
Gli effetti di bordo libero provocano un abbassamento del valore della triassialità dello
stato di sforzo ed un conseguente aumento del valore della deformazione a rottura. La
perdita di idrostaticità dello stato di sforzo limita la possibilità di prevedere in maniera
accurata i processi di cedimento utilizzando un criterio di massima pressione ed
evidenzia tutti i vantaggi di previsione garantiti dall’avere a disposizione il criterio di
danno non lineare presentato. Una conferma ulteriore dell’avvenuta variazione dello
stato di sforzo può essere ottenuta dall’analisi delle deformate dei flyers di più piccolo
diametro, figure 11a-b, 12a-b. Queste tendono ad assumere un profilo molto simile a
quelli ottenuti in un test di Taylor in cui si assume uno stato di sforzo unidimensionale.
7.4 Re-shock experiment
L’ultima configurazione analizzata è quella del re-shock experiment, che viene
realizzata posizionando, come illustrato in Figura 7.22, sulla parte posteriore del disco
proiettile, un disco di maggiore impedenza meccanica. Quando l’onda di compressione
generata dall’impatto raggiunge l’interfaccia con il backing, in accordo con le equazioni
(2.18) e (2.19), parte dell’impulso viene trasmesso e parte viene riflesso come impulso
di compressione. Tale impulso sovrapponendosi all’onda di compressione generata
dall’impatto causa, una ricompressione dello stato del materiale. Se la velocità
123
d’impatto è sufficientemente elevata da provocare un’onda d’urto, la ricompressione è
causa del fenomeno noto in letteratura col termine “re-shock”.
BACKING
PROIETTILE
BERSAGLIO
Figura 7.22 – Schema della configurazione del re-shock experiment.
PRECURSORE
ELASTICO
Figura 7.23 - Profilo di velocità misurato in un re-shock experiment, [3].
Tale configurazione è estremamente affascinante per la possibilità che offre di
investigare la risposta meccanica dei solidi in regime dinamico di shock ripetuto.
Secondo la teoria delle onde di sforzo, la risposta di un materiale elasto-plastico, al reshock dovrebbe essere interamente plastica, poiché lo stato del materiale si dovrebbe
trovare sulla superficie di snervamento. Dagli esperimenti si osserva invece la presenza
124
di un gradino, comunemente riconosciuto come un inaspettato precursore plastico, che
precede l’arrivo del ricaricamento plastico, Figura 7.23. I tentativi proposti in letteratura
per cercare di giustificare la presenza del gradino anomalo fanno tutti riferimento a
simili meccanismi fisici, che hanno luogo alla meso-scala, che dovrebbero portare lo
stato del materiale all’interno della superficie di snervamento. Lipkin e Asay, [4],
ritengono che le differenti orientazioni dei sistemi di scorrimento di grani contigui sono
causa di una deformazione, alla meso-scala, non uniforme. Essi proposero un modello,
una distribuzione dello stato di snervamento del materiale precompresso, in grado di
duplicare le caratteristiche chiave del profilo di velocità rilevato in un re-shock
experiment. Swegle e Grady, [5], credono che il gradino anomalo sia dovuto a fenomeni
di localizzazione delle deformazioni, causati da elevati gradienti termici, giustificati
dalla natura dinamica degli eventi, che si realizzano alla meso-scala. Tali fenomeni, a
loro volta, sarebbero responsabili di un comportamento dello stato di snervamento del
materiale dipendente dal tempo.
Nel presente lavoro di tesi si presenta una nuova interpretazione del fenomeno basata su
considerazioni alla macro-scala. In accordo con tale interpretazione, la giustificazione
della presenza del gradino anomalo va ricercata nella distribuzione non uniforme della
deformazione plastica, lungo lo spessore del disco bersaglio, dovuta a processi
dissipativi che hanno luogo durante il passaggio della prima onda di compressione.
7.4.1 Fenomenologia del re-shock
In un Flyer Plate Impact Test, il profilo di velocità della particella situata sulla
superficie libera presenta un andamento ondoso, più o meno pronunciato, all’inizio del
plateau a velocità costante. In termini di sforzo, per un punto prossimo alla superficie
libera, il profilo d’onda non è perfettamente quadrato, ma mostra un picco, più alto in
valore del susseguente plateau di sforzo. Ciò si traduce nel fatto che, quando lo stato di
sforzo in un punto raggiunge il plateau, questo non si trova esattamente sulla superficie
di snervamento.
Comunque, in base a tali considerazioni, per un punto prossimo alla superficie libera, il
divario non è grande abbastanza da giustificare la presenza del gradino anomalo,
all’arrivo della seconda onda di shock.
Un’onda di sforzo che viaggia nel disco bersaglio e soggetta a processi dissipativi, che
125
ne riducono l’intensità. Di conseguenza, la distribuzione di deformazione plastica, lungo
lo spessore, non dovrebbe essere uniforme, ma dovrebbe mostrare un massimo alla
superficie d’impatto ed un minimo alla superficie posteriore. Tale congettura può essere
verificata dall’analisi numerica di un Flyer Plate Impact Test standard. Un esempio del
risultato trovato è riportato in .
Figura 7.24 - Distribuzione della deformazione plastica lungo lo spessore del disco bersaglio, a
seguito dell’onda di compressione, in un Flyer Plate Impact Test.
Figura 7.25 - Profili di sforzo, calcolati numericamente, a diverse posizioni lungo lo spessore del
disco bersaglio, in un Flyer Plate Impact Test standard.
126
Figura 7.26 - Profili di sforzo, calcolati numericamente, a diverse posizioni lungo lo spessore del
disco bersaglio, in un re-shock experiment.
I profili di sforzo, a diverse posizioni lungo lo spessore del disco bersaglio, riportati in
Figura 7.25, mostrano come, durante la propagazione, essi si riducano progressivamente
nell’intensità del picco, mentre il plateau sembrano rimanere costanti. In prossimità
della superficie libera, la differenza in sforzo tra il picco ed il plateau, è molto piccola,
risultando nelle piccolissime oscillazione osservate sperimentalmente con le tracce
VISAR.
Come illustrato in Figura 7.26, in un re-shock test, la differenza in sforzo tra il picco e il
plateau, in prossimità della superficie libera, non è sufficiente da giustificare l’ampiezza
del gradino anomalo osservato negli esperimenti, mentre questa è consistente con la
differenza in sforzo calcolata in prossimità della superficie d’impatto.
In Figura 7.27, è riportata la rappresentazione dello stato di sforzo, al plateau, del punto
materiale, nel piano dei deviatori π , al fine di fornire una spiegazione del processo che
porta alla generazione del gradino anomalo. Il cerchio interno a tratto grigio continuo
rappresenta lo stato di sforzo corrente, che si trova in campo elastico per via del fatto
che il plateau è, rispetto al picco, ad un livello di sforzo inferiore; il cerchio a tratto
continuo nero indica la superficie di snervamento corrente; il cerchio più grande,
tratteggiato, rappresenta l’espansione della superficie di snervamento per l’arrivo
dell’onda di re-shock; il percorso di carico, in stato di deformazine uniassiale, è
rappresentato da vettore nero.
127
Distance from impact
Figura 7.27 - Rappresentazione dello stato di sforzo del punto materiale che a subito uno shock, a
differenti posizioni lungo lo spessore del disco bersaglio, nel piano dei deviatori π .
È importante sottolineare che, poiché il valore dello sforzo al plateau, può essere
considerato costante per tutti i punti, lo stato di sforzo iniziale, rappresentato dai cerchi
grigi, è lo stesso per tutti i punti lungo lo spessore del disco bersaglio. All’arrivo
dell’onda di re-shock, lo stato di sforzo, di un punto matriale prossimo alla superficie
d’impatto, dove la differenza in sforzo è più elevata, cresce dall’intervallo elastico fino
alla prima superficie di snervamento e poi cresce con essa: di conseguenza, a questo
punto, saranno generati un precursore elastico e un’onda plastica più lenta. Viaggiando
verso la superficie posteriore del bersaglio, l’onda di reshock tova punti in cui la
superficie di snervamento corrente è più piccola. Questo fa sì che il precursore elastico
generato in un punto precedentemente è forte abbastanza da snervare il materiale a
monte producendo un’onda plastica. La continua generazione, in accordo al
meccanismo descritto, di onde plastiche di diversa intensità porta alla creazione di un
profilo d’onda a “scalini” che, anche se può, in qualche modo, ricordare la caratteristica
struttura del precursore elastico e onda plastica, è quasi interamente plastico. Sulla
superficie libera, la struttura a scalino dell’onda di sforzo causa, nel profilo di velocità,
la comparsa del gradino anomalo che, a questo punto, dovrebbe essere indicato, più
128
propriamente, “precursore plastico”.
jump for
“plastic precursor”
Figura 7.28 - Confronto tra il profilo di velocità sperimentale, [3], e quello calcolato numericamente
con MSC.Marc.
In Figura 7.28, il profilo di velocità calcolato numericamente, per il re-shock
experiment è messo a confronto con i risultati sperimentali riportati da Vogler e Asay,
[3]. La configurazione si riferisce ad un impatto simmetrico di due dischi di alluminio
6061-T6, a 1, 715 km s ; il backing è costituito da un disco di rame e, per permettere
l’utilizzo dell’interferometria laser, è stata adottata una finestra di PMMA, secondo
quanto descritto in Asay e Chhabildas, [6].
La simulazione è stata effettuata con il codice agli elemeti fini MSC.Marc per mezzo
del metodo di Humbolt d’integrazione diretta. Si è assunto che la pressione idrostatica è
linearmente proporzionale alla compressione volumetrica, poiché con tale codice
numerico non è possibile utlizzare equazioni di stato differenti. Tutte le caratteristiche
principali della curva calcolata sono consistenti con gli esperimenti: velocità massima
della prima onda, durata e intensità del plateau, salto di velocità all’arrivo del precursore
plastico. Il ritardo della prima onda plastica, la presenza del precursore elastico e la più
lenta risalita del segnale delle onde plastiche, per il profilo di velocità calcolato, sono
dovuti alla formulazione dell’equazione di stato.
129
In Figura 7.29, è presentato il confronto tra i dati sperimentali e la curva calcolata con
l’hydrocode Autodyn. In questo caso, per riuscire a catturare le caratteristiche di del
segnale risultante da un impatto a 1, 7 km s , si è utilizzata l’equazione di stato MieGrunaisen. Si è osservato che lo smorzamento numerico suggerito dal codice è troppo
severo e non permette di catturare le caratteristiche fondamentali del profilo di velocità:
di
conseguenza
è
stato
appropriatamente
modificato.
L’abbassamento
dello
smorzamento numerico è causa di una risposta che presenta oscillazioni ad alta
frequenza non realistiche. Il profilo calcolato, ancora una volta è molto simile a quello
misurato, ma le caratteristiche chiave del precursore elastico rischiano di essere
nascoste dalle forti oscillazioni.
Figura 7.29 - Confronto tra il profilo di velocità sperimentale, [3], e quello calcolato numericamente
con Autodyn.
Infine, l’interpretazione del fenomeno proposta permette di giustificare il ritardo del
gradino anomalo rispetto alla corrispondente onda di rilascio in un Flyer Plate Impact
Test standard, che si osserva sperimentalmente, Figura 7.30. Per il test standard, il
primo abbassamento della velocità al termine del plateau è dovuto all’arrivo della prima
onda di rilascio, che è puramente elastica e viaggia alla velocità che le compete. Nel reshock experiment, secondo l’interpretazione data, l’onda di ricompressione porta alla
formazione dell’onda plastica a scalini, il cui fronte, inizialmente, viaggia alla velocità
130
dell’onda elastica e, successivamente, rallenta alla velocità dell’onda plastica, con
l’effetto provocato del ritardo dell’arrivo, sulla superficie libera, del precursore plastico.
Figura 7.30 - Ritardo del gradino anomalo rispetto alla corrispondente onda di rilascio in un Flyer
Plate Impact Test standard, [4].
131
Bibliografia
[1]
Rajendran, A.M., (1988), in Dynamic Constitutive/Failure Models (ed. A.M.
Rajendran,and T. Nicholas), AFWAL-TR-85-4009, Wright Patterson Afb, OH
[2]
Christy, S., Pak, H., e Mayers, M.A., in “Metallurgical Applications of Shock
Waves and High Strain Rate Phenomena” (eds., Murr et al.), Marcel Dekker,
New York, 1986.
[3]
Vogler, T.J. e Asay, J.R., “A distributional model for elastic-plastic behavior of
shock-loaded materials”, in Shock Compression of Condensed Matter, (M.D
Furnish, Y.M. Gupta, J.W. Forbes, eds.), part I, pp. 617-620, 2003.
[4]
Lipkin, J. e Asay, J.R., “Reshock and release of shock-compressed 6061-T6
aluminum”, J. of Applied Physics 48, 182, 1977.
[5]
Swegle, J. W. e Grady, D. E., “Calculation of thermal trapping in shear bands”,
Metallurgical application of shock wave and high-strain-rate phenomena”,
edited by L. E. Murr et al., New York, 1986.
[6]
Asay, J.R. e Chhabildas, L.C. “Determination of the shear strength of shock
compressed 6061-T6 aluminum”, in Shock Waves and High-Strain-Rate
Phenomena in Metals (M.A Mayers and L.E. Murr, eds.), pp. 417-431, Plenum,
New York, 1981.
132
8 Conclusioni
Nel presente lavoro di tesi, gli strumenti della simulazione numerica sono stati utilizzati
per l’analisi di tre configurazioni classiche per la caratterizzazione meccanica dei
materiali in regime dinamico: il Taylor Test, la Hopkinson Bar e il Flyer Plate Impact
Test. Gli sforzi non sono stati indirizzati alla semplice riproduzione delle caratteristiche
osservate negli esperimenti, ma con un’analisi critica si è, di volta in volta, cercato di
interpretare i processi di deformazione e rottura che si verificano in dinamica
dell’impatto, al fine di giungere ad una loro corretta modellazione.
Il contributo innovativo del presente lavoro può essere sintetizzato nei seguenti punti:
‰
è stata dimostrata la capacità di previsione di un modello di danno duttile non
lineare in regime di elevata velocità di deformazione;
‰
sono stati individuati, per il cilindro di Taylor, due diversi modi di rottura che
si realizzano a tempi diversi, durante il processo di deformazione, e per stati
della triassialità dello stato di sforzo differenti;
‰
è stata individuata una correlazione tra la deformazione di soglia, uno dei
coefficienti del modello di danno non lineare, e la dimensione media del grano;
‰
è stato individuato e quantificato, con gli strumenti della meccanica della
frattura, un processo di dissipazione nel meccanismo di separazione delle
superfici di rottura per spalling;
‰
è stata fornita una nuova interpretazione della presenza gradino anomalo che si
osserva nel re-shock experiment, identificato in letteratura come un inaspettato
“precursore elastico”;
‰
è stato dimostrato che, in realtà tale gradino è plastico;
‰
l’interpretazione proposta ha permesso di giustificare il ritardo del gradino
anomalo, rispetto alla corrispondente onda di rilascio in un Flyer Plate Impact
Test standard, che si osserva sperimentalmente.
133