UNIVERSITÀ POLITECNICA DELLE MARCHE

DIPARTIMENTO DI INGEGNERIA INDUSTRIALE E SCIENZE MATEMATICHE

CORSO DI LAUREA MAGISTRALE IN INGEGNERIA MECCANICA

Tesi di laurea

Modellazione costitutiva del comportamento dinamico di

sughero agglomerato

Constitutive modelling of the dynamic behaviour of

agglomerated cork material

Relatore: Laureando:
Prof. Marco Sasso Davide Podeschi

ANNO ACCADEMICO 2021‐2022

I
 Ringraziamenti
Ringrazio il professor Marco Sasso per l’opportunità di
tesi e la disponibilità dimostrata.

Ringrazio mio fratello Federico e mia cognata Sara per

il continuo appoggio durante questo lungo percorso di
studio.

Ringrazio la mia compagna Lara per tutto l’amore che

mi ha dato in questo difficile periodo.

Ringrazio mia mamma e mio babbo per tutto quello che

hanno fatto per me, il loro esempio sarà sempre la base
della mia vita.
Dedico a voi e ai nonni questa tesi.

II
 Indice

Indice…………………………………………………………………………………………………………………………….……Ⅲ
Introduzione………………………………………………………………………………………………………………….…….Ⅴ
Capitolo I ⸺ Il sughero…………………………………………………………………………………………………...…1
I.1 Generalità ............................................................................................................... 1
I.2 Struttura ................................................................................................................. 2
I.3 Composizione chimica ............................................................................................ 4
I.4 Applicazione ........................................................................................................... 5
I.5 Proprietà fisiche ..................................................................................................... 5
Capitolo II ⸺ Modello Iperelastico…………………………………………………………………………….…......8
II.1 Introduzione .......................................................................................................... 8
II.2 Potenziale elastico per materiali isotropi .............................................................. 9
II.3 Legame sforzi‐deformazioni ................................................................................ 12
II.4 Forme particolari di potenziale elastico .............................................................. 16
II.5 Calibrazione dei modelli di potenziale elastico per materiali elastomerici......... 21
Capitolo III ⸺ Modello Viscoelastico………………………………………………………………………………..28
III.1 Introduzione ....................................................................................................... 28
III.2 Rilassamento e scorrimento viscoso .................................................................. 30
III.3 Modelli viscoelastici elementari: modello di Maxwell e Kelvin ......................... 33
III.4 Modello di Maxwell generalizzato ..................................................................... 38
Capitolo IV ⸺ Metodo a elementi finiti e software di simulazione……………………………….....44
IV.1 Introduzione ....................................................................................................... 44
IV.2 Metodo a elementi finiti FEM ............................................................................ 48
IV.3 Analisi lineari e non lineari ................................................................................. 50
IV.4 Algoritmo esplicito ............................................................................................. 54
IV.5 Tempo di incremento stabile ............................................................................. 58
Capitolo V ⸺ Implementazione del modello e risultati delle analisi FEM……………….……..…60
V.1 Introduzione ........................................................................................................ 60

III
 V.2 Macchina per test di perforazione ...................................................................... 60
V.3 Implementazione del modello iperelastico e viscoelastico su Abaqus ............... 62
V.4 Effetto Mullins ..................................................................................................... 66
V.5 Implementazione del modello di danneggiamento ............................................ 70
V.6 Preparazione del modello FEM ........................................................................... 75
V.7 Analisi FEM e risultati .......................................................................................... 77
V.8 Mappe di deformazione FEM .............................................................................. 82
Capitolo VI ⸺ Conclusioni…………………………………………………………………………………………….....86
Riferimenti bibliografici……………………………………………………………………………………………………..88

IV
Introduzione

Negli ultimi decenni i solidi cellulari come le schiume di origine

polimerica hanno trovato larga diffusione grazie alle loro peculiarità. Tra
le caratteristiche che rendono tali materiali interessanti nel contesto
dell’ingegneria meccanica si possono sicuramente indicare la capacità di
assorbire energia da impatto, lo smorzamento di vibrazioni meccaniche, la
bassa densità relativa, l’isolamento termo-acustico, l’elevata rigidezza
specifica (intesa come rapporto tra il modulo elastico a compressione e la
densità). D’altra parte però, una crescente coscienza ambientalista e la
transizione ecologica già in atto a livello di politiche europee, sta ponendo
le basi per lo studio e l’utilizzo di materiali alternativi a minor impatto
ambientale; tutti i materiali di derivazione polimerica infatti, sebbene
versatili e di facile produzione, presentano l’inconveniente dello
smaltimento, con relativi costi. Il settore di utilizzo dei solidi cellulari
comprende i dispositivi di sicurezza anticrash usati in ambito sportivo, il
settore del packaging, l’applicazione come anima strutturale interna di
pannelli sandwich-composito (settori automotive ed aerospace), ecc. Fra i
materiali cellulari di origine naturale alternativi alle schiume strutturali
sintetiche, per quel che riguarda le caratteristiche meccaniche sopra
indicate, vi rientra indubbiamente il sughero. Visti i presupposti di utilizzo
e le normative in materia di sicurezza sempre più severe, vi è la necessità
di conoscere il comportamento meccanico del sughero qualora sia
soggetto a carichi dinamici. Con questa tesi viene studiato numericamente,
mediante software di simulazione “Abaqus/Explicit” di Dassault Systeme,
il comportamento del materiale sottoposto a un test di penetrazione con tre
diverse velocità d’impatto.

V
Com’è ben noto, anche semplicemente osservando il tappo di una
bottiglia, il sughero è un materiale che si presta a subire grandi
deformazioni e a recuperare tali deformazioni quasi integralmente; tale
comportamento non lineare viene definito iperelasticità. Partendo da un
modello iperelastico si sono aggiunti via via ulteriori modelli allo scopo di
perfezionare la risposta del materiale nel campo delle elevate velocità di
deformazione, compreso la modalità di rottura, ed avendo come riscontro
le curve sperimentali velocità-tempo di prove di penetrazione da impatto.

VI
 Capitolo I

Il sughero
I.1 Generalità

Il sughero è un materiale molto diffuso e a tutti ben noto, ma non

altrettanto ben conosciuto. L’uso che lo rende tanto comune è quello di
sigillare bottiglie, ma è molto di più: la sua struttura microscopica è
sorprendente, rendendolo campione di isolamento meccanico, acustico,
termico ed elettromagnetico. Ottenuto dalla corteccia di querce da
sughero, la sua estrazione avviene nel periodo che va dai primi di maggio
a fine agosto, quando il sughero distacca più facilmente senza causare
danni alla pianta. La prima decorticazione di una giovane quercia, la
demaschiatura, si effettua quando la pianta ha circa 25-30 anni e una
circonferenza non inferiore ai 60 cm, e se ne ottiene un sughero da macina
detto sugherone o maschio. Le successive estrazioni avvengono a
intervalli di almeno dieci anni, come previsto dalle normative, ma anche
12-13 se il sughero non ha raggiunto un calibro accettabile. Il prodotto
ottenuto è detto sughero gentile e viene utilizzato, se di buona qualità, per
la fabbricazione dei tappi. Le principali aree di diffusione della quercia da
sughero sono il Portogallo, la Spagna, la Sardegna, la Sicilia, la Maremma
grossetana, la Corsica, il sud della Francia, il Nordafrica e la Cina. (Gil et
al.,1998; Barberis et al.,2003; Pereira et al.,1984; Fialho et al.,2001;
Carrasquinho et al.,1987; Pereira et al.,1987; Costa et al.,2003).

1
I.2 Struttura

Il sughero è un materiale cellulare a celle chiuse che riempie tutti gli spazi,
senza lasciare vuoti intercellulari. Come tutti i materiali di questo tipo le
sue proprietà dipendono fortemente dalla geometria e dalla distribuzione
delle celle. Caratteristica è la presenza di una struttura regolare a celle
chiuse con una crescita regolare in direzione radiale al tronco dell’albero.
Possono esservi durante la fase di crescita piccole variazioni nella densità
e nelle dimensioni delle celle causate dal ritmo fisiologico della pianta,
portando alla creazione di diversi anelli di crescita. La regolarità della
struttura può anche essere compromessa da alterazioni locali dovute a
diversi fattori biologici ed accidentali quali per esempio fratture o
inclusioni del legno (Pereira et al.,2007). Per descrivere correttamente la
struttura del sughero sarà necessario conoscerne la corretta collocazione
spaziale all’interno della pianta, definendo con la sezione trasversale
quella perpendicolare alla direzione assiale (estensione verticale
dell’albero), la sezione tangenziale è perpendicolare a quella radiale
(direzione di sviluppo della sezione di tronco), ed infine la sezione radiale
che contiene la direzione assiale ed è perpendicolare alla direzione
tangenziale (figura 1).

Figura 1

Nella sezione tangenziale le cellule del sughero sono organizzate secondo

una struttura a nido d’ape, senza spazi vuoti tra le cellule contigue, con i
singoli poligoni costituti da elementi da quattro ad otto lati (Gibson et al.,

2
2005; Fialho et al.,2001; Pereira et al.,1987). In direzione radiale e
trasversale invece la struttura appare completamente diversa a quella
tangenziale. Le singole celle, infatti, richiamano dei mattoncini
rettangolari allineati con i tre lati che si incontrano nei singoli vertici.
Tridimensionalmente possiamo quindi pensare al sughero come costituito
da celle riconducibili a prismi esagonali impilati in righe (figura 2).

Figura 2

All’interno di ciascuna riga tutte le celle sono costituite dalla stessa base,
ma differente altezza del prisma. Le varie righe sono tra loro allineate in
direzione radiale e disposte con le basi del prisma parallelamente alla
sezione tangente. Nelle file adiacenti le basi del prisma non sono
coincidenti e giacciono tra loro in posizione sfalsata. In conclusione,
quindi, le cellule di sughero sono poliedri che entrano tra loro in contatto
in bordi facce e vertici. In particolare, due cellule possono entrare in
contatto in una faccia, tre in un bordo e quattro in un vertice, di modo che
in una struttura ideale con i prismi tutti uguali, ogni cellula entrerebbe in
contatto con ben 14 cellule. Le dimensioni sono molto piccole, per una
cella media normale parliamo di altezza del prisma di appena 40 um, bordo
di base 20 um e spessore della parete di appena 1 um. Nella realtà le
dimensioni sono molto variabili in relazione anche al ritmo di crescita
stagionale della pianta, in tal senso si possono per esempio distinguere
all’interno di uno stesso anello annuale tra: quelle prodottesi nelle prime

3
fasi di crescita da quelle nelle ultime fasi. Le prime risultano infatti più
grandi e più sottili le seconde più corte e spesse. Altro dato caratteristico
è il volume realmente occupato da materiale nel sughero, si stima infatti
che la porzione solida occupa solo l’8-9% nelle cellule del primo periodo
ed il 15-22% in quelle del secondo. (Pereira et al., 2007)
Ogni cellula si origina a partire di una cellula madre presente nel
phellogen, la direzione di crescita è quella radiale. Rispetto ad altre piante
il periderma del sughero presenta caratteristiche singolari di crescita,
regolarità e longevità. Uno strato protettivo denominato phellem,
costituito da cellule morte che fanno parte del periderma costituente la
corteccia, separa le cellule della pianta dall’ambiente esterno. Dopo la
morte del phellogen, si verifica, infatti, una rapida formazione di tessuto
filogenico traumatico che riprende la sua funzione come produttore di
strato di sughero. E’ questa caratteristica che dà alle piante da sughero la
possibilità di poter essere usate come produttrici di sughero in maniera
sostenibile.

I.3 Composizione chimica

La composizione chimica è fortemente dipendente da diversi fattori:
l'origine geografica, il clima, le condizioni del suolo, l’origine genetica, la
dimensioni dell'albero, l’età (vergine o riproduzione) e le condizioni di
crescita. La struttura della parete è costituita da lamelle centrali ricche di
lignina (un polimero reticolato, parzialmente aromatico e composto da
alcoli derivati dall'1-fenilpropano responsabile della struttura della parete)
e da una parete secondaria costituita da suberina (un copolimero aromatico
formato da monomeri di esteri con catene a 18-30 atomi di carbonio,
responsabile della comprimibilità ed elasticità del sughero) alternata a
lamelle di cera. Secondo alcuni studi la parete secondaria appare
lignificata, presentando quindi anche tracce di cellulosa. Troviamo infine

4
una terza parete di polisaccaridi. Tra i vari componenti, la suberina è la più
abbondante (circa il 40%), la lignina corrisponde al 22%, i polisaccaridi al
18% e gli estraibili al 15%. Internamente, invece, le cellule contengono
cristalli di cerina e fridelina ed una grande quantità di gas. (Pereira et
al.1979; Pereira et al., 1981; Bernards et al., 1998; Bento et al., 1992;
Cordeiro et al., 1995; Pereira et al., 1988).

I.4 Applicazione

Oggi il sughero è stato riscoperto in applicazioni tecnologicamente

avanzate grazie alla sua elaborata struttura e composizione. Viene
utilizzato in strutture leggere, nel packaging, nell’isolamento termo-
acustico, come supporto per catalizzatori, applicazioni biomediche,
sicurezza passiva in campo veicolistico, etc. Il sughero è adatto a
migliorare le capacità dissipative dei polimeri rinforzati con fibre
mediante strutture a sandwich in cui è utilizzato come core.

I.5 Proprietà fisiche

Il sughero ha una straordinaria combinazione di proprietà. È leggero ma
resistente (un tessuto vegetale che pesa appena 0.16 grammi per
centimetro cubico); è un ottimo isolante per calore e suono (Questo perché
i componenti gassosi contenuti nel sughero sono richiusi in piccoli
compartimenti stagni, isolati l’uno dall’altro da una sostanza resistente
all’umidità); ha un alto coefficiente di attrito; è impermeabile ai liquidi ed
ai gas (grazie alla suberina e ai ceroidi contenuti nelle pareti cellulari);
chimicamente stabile e resistente al fuoco(brucia senza produrre fiamma e
non emette gas tossici durante la combustione). Non assorbendo polvere è
considerabile come materiale ipoallergenico. Grazie alla sua struttura
inoltre è molto resistente alle abrasioni ed agli impatti, avendo un alto

5
coefficiente di frizione (Gibson et al., 1981). Ha inoltre mostrato una
elevata affinità per i liquidi non polari ed a polarità molto bassa ed una
buona compatibilità con polimeri acidi e basici (grazie alla sua superficie
dalle caratteristiche anfotere) (Gibson et al., 1997).
In generale il sughero è caratterizzato da un peso specifico molto basso,
ma fortemente variabile entro ampi limiti in funzione del trattamento
subito (naturale o bollito) e della sua età (vergine o riproduzione). Si
possono addirittura registrare variazioni di un fattore due passando da
circa 120 a 240 kg/m^3. (Cumber et al., 2000). Globalmente le variazioni
di densità sono riconducibili alle dimensioni della cella (corrugazione,
altezza e spessore della parete) ed alla frazione di volume occupata dai
canali lenticolari. Alte densità, per esempio, coincidono con celle spesse e
fortemente ondulate caratterizzate dalla presenza di pochi canali (Cumber
et al., 2000; Gibson et al., 1981; Fortes et al., 1988; Baptista et al.,1993).
In quest’ottica le cellule formate durante la primavera sono più sottili e
quindi meno dense di quelle autunnali. Ad influire ancora sulla bassa
densità anche la presenza di gas contenuto nelle celle. Quest’ultimo,
insieme alle ridotte dimensioni della cella sono anche alla base della
trasmissione del calore, prevalentemente per conduzione, con pareti che
presentano una conducibilità simile a quella del gas e quindi fortemente
ridotta. Allo stesso modo anche la trasmissione del suono è molto scarsa
in quanto la bassa densità e l’elevata porosità assorbono e trasformano in
energia termica la maggior parte delle onde (Mano et al., 2002).
L’attrito si aggira intorno allo 0,5 e come molte altre proprietà presenta un
comportamento anisotropo (che si attenua in presenza di sughero umido),
essendo maggiore per gli scorrimenti sul piano tangenziale rispetto ai piani
ad esso perpendicolari. Risulta influenzato da fenomeni di adesione tra le
aree di contatto e dalla deformazione per flessione delle pareti cellulari
(Vaz et al., 1998).

6
Il materiale risulta essere caratterizzato da un comportamento anisotropo,
valutabile dal rapporto tra la dimensione maggiore della cella e da quella
minore (indicato con R) ed il cui valore influenza fortemente le proprietà
del materiale. Nel caso del sughero, anche se tutt’ora non sono ancora note
tutte le cause, è possibile ricondurre l’anisotropia a due fattori principale
tra loro indipendenti: l’anisotropia della struttura, dovuta alla forma delle
celle, e l’anisotropia del materiale dovuta alla parete delle celle.
L’anisotropia strutturale in particolare deriva dalla non equidimensionalità
delle celle, per le sezioni trasversali e radiali è possibile valutare, infatti,
rapporti di anisotropia di 1,5-1,7, nella sezione tangenziale invece il
rapporto è circa 1,0-1,1. Presenta simmetria ortotropica, il che significa
che ha tre piani di simmetria ortogonali, anche se è possibile approssimarlo
ad un materiale con simmetria circolare nel piano tangenziale, ovvero
rispetto tale piano può essere visto come materiale isotropo ed
assialsimmetrico rispetto alla direzione radiale. Per l’anisotropia
materiale, risulta in gran parte sconosciuta, anche se in prima
approssimazione può essere considerato isotropo per la presenza di
cellulosa come elemento strutturale nella parete cellulare. Ad influenzare
ci sono ancora la presenza di anelli di crescita, con la presenza di sottili
strati di latercork e la formazione di canali lenticolari che non avendo le
stesse dimensioni nelle differenti direzioni, introduce un fattore di
variazione nel piano tangenziale (Pereira et al., 2007)

7
 Capitolo II

Modello Iperelastico
II.1 Introduzione

L’elasticità dipende dall’esistenza di un potenziale scalare, funzione solo

dello stato del materiale, che rappresenta l’energia elastica accumulata in
un processo deformativo ed eguaglia il lavoro di deformazione. Quando il
comportamento è elastico e al tempo stesso lineare, si dimostra, per
materiali isotropi, che la legge costitutiva dipende da due soli parametri. Il
legame costitutivo prende il nome di legge di Hooke. Nei materiali
metallici il regime elastico è caratterizzato da un comportamento
sostanzialmente lineare e termina con lo snervamento. Il campo elastico è
limitato a piccole deformazioni e la teoria elastica è generalmente
formulata utilizzando il tensore delle deformazioni infinitesime.
Alcuni materiali, quali gli elastomeri, mostrano però un comportamento
elastico anche per grandi deformazioni che è caratterizzato da un legame
niente affatto lineare. Per questi materiali la legge costitutiva è basata
sull’identificazione di un’opportuna forma per la funzione di potenziale
elastico, che dovrà poter essere espressa attraverso misure di deformazione
valide per deformazioni e rotazioni finite. Tali materiali si definiscono
iperelastici e questo capitolo presenta gli aspetti salienti teorici e pratici
della modellazione del loro comportamento costitutivo. Si presenterà in
primo luogo una forma generale del potenziale elastico per materiali
isotropi, in grado di separare il potenziale associato alla componente
deviatorica della deformazione da quella volumetrica. Tale aspetto è di
particolare interesse per gli elastomeri che hanno tipicamente un
comportamento quasi incomprimibile.

8
Saranno quindi presentate le forme corrispondenti delle leggi costitutive,
basate sull’applicazione di diverse misure di deformazione. Saranno
inoltre introdotte alcune forme specifiche del potenziale proposte in
letteratura e, infine, si discuteranno le procedure di calibrazione per le
formulazioni individuate, a partire da dati ottenuti mediante prove
sperimentali.

II.2 Potenziale elastico per materiali isotropi

Considerando la necessità di definire la legge costitutiva anche in regime
di grandi deformazioni e rotazioni, il lavoro di deformazione può essere
espresso utilizzando il tensore di deformazione di Green-Lagrange. Per
esprimere il lavoro di deformazione è necessario utilizzare una misura di
sforzo che sia energeticamente coniugata con il tensore di Green-
Lagrange, rappresentata dal secondo tensore di Piola-Kirchoff. Il tensore
di Piola-Kirchoff è legato al tensore di Cauchy, che rappresenta gli sforzi
veri sulla sezione deformata, attraverso l’espressione

II.1
𝟏 𝑻
𝜮Ⅱ 𝐽𝑭 𝝈𝑭

dove F è il gradiente di deformazione e 𝐽 è il determinante di F. 𝐽

rappresenta il rapporto fra il volume del materiale nella configurazione

deformata e quella originale indeformata 𝐽 . Il lavoro di

deformazione per unità di volume nel volume indeformato, per un

incremento di deformazione 𝑑𝑬 risulta

𝑑𝐿 II.2
𝑙 𝜮Ⅱ ∶ 𝑑𝑬
𝑑𝑉
dove il prodotto scalare rappresentato dal simbolo “:” corrisponde a
moltiplicare fra di loro i termini a componenti identici dei due tensori
rappresentati dalle matrici 𝜮Ⅱ ed E. L’assunzione di materiale con

9
comportamento elastico implica che il lavoro compiuto durante un
processo deformativo dallo stato di deformazione 𝑬 allo stato 𝑬 risulti
uguale alla variazione di potenziale elastico 𝜔 𝑬 :

𝑬 𝑬
𝛥𝑙 𝛥𝜔 𝑬 𝜮Ⅱ : 𝑑𝑬 𝑑𝜔 𝑬 II.3
𝑬𝒂 𝑬

Ciò implica che le componenti del secondo tensore di Piola-Kirchoff

possono essere ottenute derivando il potenziale elastico:

𝜕𝜔 𝑬 II.4
𝜮Ⅱ
𝜕𝑬
A differenza del caso lineare, tuttavia, il potenziale elastico non è una
forma quadratica delle deformazioni, ma può essere rappresentato da una
funzione generica. Nell’ipotesi di materiale isotropo, comunque, il
potenziale è convenientemente espresso come funzione degli invarianti del
tensore di deformazione, poiché il suo valore non può dipendere dal
sistema di riferimento considerato per descrivere lo stato di deformazione
stesso. La forma più conveniente da dare al potenziale elastico per un
materiale iperelastico ortotropo è tale da distinguere i contributi dovuti alla
componente deviatorica e volumetrica della deformazione. Il
disaccoppiamento fra le leggi costitutive che legano sforzi e deformazioni
volumetrici e deviatorici è possibile in un materiale isotropo, e comporta
importanti vantaggi nella modellazione e calibrazione del comportamento
dei materiali elastomerici. Tali materiali, infatti, hanno spesso un
comportamento quasi incomprimibile, ed il contributo al potenziale
originato dalle deformazioni volumetriche è quantitativamente e
qualitativamente diverso rispetto a quello relativo alle componenti di
deformazione deviatorica. Per ottenere la separazione fra le due
componenti il potenziale potrà essere espresso in funzione degli invarianti
del tensore sinistro di Cauchy-Green B, definito come:

10
 𝑩 𝑭 ∙ 𝑭𝑻 II.5

Il potenziale sarà espresso in funzione dei primi due invarianti di B,

espressi nell’equazione (II.6), e dal determinante del gradiente di
deformazione 𝐽

𝐼𝑩𝟏 𝑡𝑟 𝑩 𝑰∶𝑩
II.6
1
𝐼𝑩𝟐 𝐼 𝑡𝑟 𝑩: 𝑩
2 𝑩𝟏
Gli invarianti del tensore B si riferiscono alla sola componente deviatorica
dello stato di deformazione. Si consideri, infatti, un processo di
deformazione infinitesimo caratterizzato da un gradiente di spostamento
infinitesimo 𝑑𝒖, il quale è possibile mettere in relazione con l’incremento
infinitesimo del gradiente di deformazione:
𝜕𝑑𝒖
𝑑𝑳
𝜕𝑋
II.7
𝜕𝑑𝒖
𝑑𝑭 𝑰 𝑑𝑳 𝑰
𝜕𝑋
L’incremento infinitesimo di deformazione è ottenuto dalla parte
simmetrica di L (la cui derivata nel tempo è il tensore di velocità di
deformazione):
1
𝑑𝑫 𝑑𝑳 𝑑𝑳𝑻 II.8
2
L’incremento infinitesimo di deformazione volumetrica 𝑑𝜀 e
deviatorica 𝑑𝒆 risultano:

𝑑𝜀 𝑡𝑟 𝑫 𝑰∶𝑫
II.9
1
𝑑𝒆 𝑑𝑫 𝑑𝜀 𝑰
3
Nel processo l’incremento di 𝐽 dipende solo dall’incremento di
deformazione volumetrica:

11
 𝑑𝐽 𝐽𝑑𝜀 II.10

Considerando l’eq. 7, l’eq. 8, l’eq. 9 e l’eq. 10, si ottiene che l’incremento

dei due invarianti 𝐼𝑩𝟏 e 𝐼𝑩𝟐 dipende solo dall’incremento di deformazione
deviatorica 𝑑𝒆, ovvero
𝑑𝐼𝑩𝟏 2𝑩 ∶ 𝑑𝒆
II.11
𝑑𝐼𝑩𝟐 2 𝐼𝑩𝟏 𝑩 𝑩 ∙ 𝑩 ∶ 𝑑𝒆
Pertanto, definendo il potenziale elastico 𝜔 in funzione di 𝐼𝑩𝟏 , 𝐼𝑩𝟐 e 𝐽, si
riescono a separare gli effetti delle componenti di deformazione
volumetrica e deviatorica. Considerando le equazioni precedenti si ottiene:

𝜔 𝜔 𝐼𝑩𝟏 , 𝐼𝑩𝟐 , 𝐽
𝜕𝜔 𝜕𝜔 𝜕𝜔
⟹ 𝑑𝜔 𝑑𝐼 𝑑𝐼 𝑑𝐽
𝜕𝐼𝑩𝟏 𝑩𝟏 𝜕𝐼𝑩𝟐 𝑩𝟐 𝜕𝐽 II.12

𝑑𝜔 2 𝐼𝑩𝟐 𝑩 𝑩 ∙ 𝑩 ∶ 𝑑𝒆 𝐽 𝑑𝜀
𝑩𝟏 𝑩𝟐 𝑩𝟐

II.3 Legame sforzi-deformazioni

L’incremento di potenziale elastico è stato espresso in eq. 12 attraverso gli
incrementi infinitesimi di deformazione, che si ottengono dal tensore di
velocità di deformazione. Grazie a tale risultato è possibile ottenere le
espressioni degli sforzi energeticamente coniugati al tensore dD che sono,
in effetti, le componenti del tensore degli sforzi di Cauchy. Infatti il lavoro
di deformazione per un processo infinitesimo, riferito al volume nella
configurazione deformata e indeformata, si esprime come:

𝑑𝐿 𝒔: 𝑑𝒆 𝑝𝑑𝜀 𝑑𝑣 ⟹

II.13

⟹ 𝐽 𝒔: 𝑑𝒆 𝑝𝑑𝜀 𝑑𝑉

12
dove s è la componente deviatorica del tensore di Cauchy e p è la
corrispondente componente idrostatica, positiva se di compressione,
energeticamente coniugata con la deformazione volumetrica. L’eq. 13
sfrutta la possibilità di disaccoppiare il lavoro di deformazione deviatorico
e volumetrico per i materiali isotropi. Inoltre nell’equazione si è fatto uso
del legame 𝐽 per trasformare l’integrale riferito al volume deformato

in quello eseguito nel volume originale. Poiché dLD = dω per l’esistenza

del potenziale elastico, le componenti dello sforzo di Cauchy posso essere
ottenute, applicando le eq. 10, eq. 12 e eq. 13:

2 𝜕𝜔 𝜕𝜔 𝜕𝜔
𝒔 𝐼𝑩𝟐 𝑩 𝑩∙𝑩
𝐽 𝜕𝐼𝑩𝟏 𝜕𝐼𝑩𝟐 𝜕𝐼𝑩𝟐
II.14
𝜕𝜔
𝑝
𝜕𝐽
L’eq. 14 consente dunque di calcolare le componenti di sforzo di Cauchy
a partire dal gradiente di deformazione nota l’espressione del potenziale,
applicando la procedura riassunta di seguito in Figura 3.

Figura 3 – Procedura di calcolo degli sforzi per un

materiale iperelastico

13
Noti gli sforzi di Cauchy, le altre misure di sforzo (ad esempio il tensore
di sforzo nominale o il secondo tensore di Piola-Kirchoff) sono
immediatamente calcolabili noto il gradiente di deformazione. Tuttavia è
possibile mettere direttamente in relazione il potenziale con altre misure
di deformazione e ottenere, derivando, le corrispondenti misure di sforzo
coniugate. In particolare, per la calibrazione delle leggi iperelastiche è
utile mettere il potenziale in funzione dei rapporti di allungamento
principali 𝜆𝟏 , 𝜆𝟐 e 𝜆𝟑 . Tali rapporti definiscono, in assi principali di
deformazione, tutte le componenti non nulle del gradiente di
deformazione, che risulta
𝜕𝑢
⎡ 1 0 0 ⎤
𝜆 0 0 ⎢𝜕𝑋 ⎥
⎢ 0 𝜕𝑣
𝑭 0 𝜆 0 1 0 ⎥
0 0 𝜆 ⎢ 𝜕𝑌 ⎥ II.15
⎢ 𝜕𝑤 ⎥
⎣ 0 0
𝜕𝑍
1⎦

⟹ det 𝑭 𝜆 𝜆 𝜆 𝐽
Sulla base dell’eq.15 il tensore di deformazione di Green-Lagrange risulta:

1
𝑬𝜶 𝜆 1 𝛼 1; 2; 3 II.16
2
I rapporti di allungamento principali sono immediatamente correlabili alla
deformazione nominale o ingegneristica:

𝜀̃𝜶 𝜆 1 𝛼 1; 2; 3 II.17

Considerando l’esempio di una prova sperimentale di trazione uniassiale,

la deformazione nominale si ottiene dividendo la variazione di lunghezza
per la lunghezza iniziale ∆𝑙⁄𝑙 . Il rapporto di allungamento è invece pari
al rapporto fra la lunghezza deformata e quella in deformata. Sviluppando
in serie le deformazioni di Green-Lagrange è anche possibile osservare

14
che, per piccole deformazioni, le componenti di E sono approssimabili con
le deformazioni ingegneristiche.
Si può inoltre osservare che, sempre in assi principali di deformazione, la
variazione di deformazione nominale consente di calcolare le componenti
della variazione del tensore di deformazione:

𝑑𝜆 0 0 𝑑𝜀 0 0
𝑑𝑭 0 𝑑𝜆 0 0 𝑑𝜀 0 II.18
0 0 𝑑𝜆 0 0 𝑑𝜀
Dal punto di vista energetico la relazione in eq. 18 è particolarmente
significativa perché il gradiente di deformazione è energeticamente
coniugato al tensore di sforzo nominale:

𝑑𝐿 𝜮𝒏 ∶ 𝑑𝑭 𝑑𝑉 II.19

Esprimendo il potenziale in funzione dei rapporti di allungamento

principale si potrà, derivando, ottenere il legame costitutivo in termini di
sforzo nominale – deformazione ingegneristica, e questo aspetto potrà
essere vantaggiosamente sfruttato per la calibrazione dei parametri da cui
il potenziale elastico mediante la correlazione con i dati sperimentali. In
particolare, analizzando lo stato di deformazione specifico ottenuto in
prove dove le direzioni principali di deformazione siano note, sarà
possibile esprimere gli invarianti 𝐼𝑩𝟏 e 𝐼𝑩𝟐 in funzione dei rapporti di
allungamento principali 𝜆𝟏 , 𝜆𝟐 e 𝜆𝟑 . Conseguentemente, derivando il
potenziale in funzione dei rapporti di allungamento principali, si potranno
ottenere gli sforzi ingegneristici. Si consideri, infine, che gli invarianti 𝐼𝑩𝟏
e 𝐼𝑩𝟐 possono essere facilmente espressi in funzione dei valori principali
della deformazione deviatorica, 𝜆𝑩𝟏 , 𝜆𝑩𝟐 e 𝜆𝑩𝟑 :

𝐼𝑩𝟏 𝜆𝑩𝟏 𝜆𝑩𝟐 𝜆𝑩𝟑 ; 𝐼𝑩𝟐 𝜆𝑩𝟏 𝜆𝑩𝟐 𝜆𝑩𝟑 II.20

15
II.4 Forme particolari di potenziale elastico
Dal punto di vista pratico la legge costitutiva iperelastica deve poter essere
calibrata attraverso una serie di prove sperimentali, specificando
completamente la funzione potenziale elastico. La Figura 4 mostra la
risposta sperimentale a trazione di un materiale elastomerico correlata con
la predizione ottenuta da un modello analitico.

Figura 4 – Correlazione teorico-sperimentale

della risposta a trazione di un materiale
elastomerico

Per poter procedere alla calibrazione è quindi necessario definire delle

forme della funzione potenziale dipendenti da una serie di parametri, che
possano essere calibrate attraverso la correlazione con le curve
sperimentali. La letteratura fornisce delle forme particolari di riferimento,
alcune delle quali sono brevemente presentate in questo paragrafo.

 Forma polinomiale
Considerando l’eq. 20, gli invarianti del tensore sinistro di Cauchy-Green
hanno valore 3 per rapporti di allungamento principali unitari; quindi in
assenza di deformazione. Analogamente, il determinante del tensore di
deformazione 𝐽 è pari ad 1 in assenza di deformazione. Tali aspetti vanno
tenuti in considerazione per definire la forma del potenziale elastico, in
modo da annullare il potenziale per deformazioni nulle. Una forma

16
polinomiale del potenziale elastico, funzione degli invarianti 𝐼𝑩𝟏 , 𝐼𝑩𝟐 e 𝐽,
è pertanto definibile nel seguente modo:

1
𝜔 𝐶 𝐼𝑩𝟏 3 𝐼𝑩𝟐 3 𝐽 1 II.21
𝐷

Tipicamente, quando entrambi gli invarianti 𝐼𝑩𝟏 e 𝐼𝑩𝟐 sono presi in

considerazione nella forma polinomiale, approssimazioni con ordine
N > 2 sono raramente usate. I coefficienti del primo ordine definiscono
una funzione lineare negli invarianti di deformazione deviatorica che, a
loro volta, sono funzioni dei quadrati dei valori principali di deformazione
deviatorica secondo l’eq. 20. Considerando che le componenti di sforzo
sono le derivate del potenziale, la risposta individuata dai coefficienti del
primo ordine è lineare, e corrisponde al coefficiente di rigidezza a taglio
iniziale 𝐺 del materiale

II.22
𝐺 2 𝐶 𝐶

I coefficienti 𝐷 definiscono la dipendenza del potenziale dalla

deformazione volumetrica. Il coefficiente di primo ordine da luogo ad un
contributo quadratico in 𝐽 che corrisponde a un legame costitutivo lineare
fra pressione e deformazione volumetrica. Il bulk modulus iniziale 𝐾 del
materiale è definito, pertanto, dal coefficiente di primo ordine secondo la
relazione:

2 II.23
𝑘
𝐷
Minore è il valore di 𝐷 , maggiore risulta il bulk modulus 𝐾 . In generale,
bassi valori di 𝐷 configurano un materiale incomprimibile. Il coefficiente
di Poisson iniziale del materiale può essere fornito dalla seguente tabella:

17
 Tabella 1

Come si può dedurre dalla tabella 1 e dalle eq. 22 e 23, maggiore è il

rapporto fra i coefficienti ed i coefficienti 𝐶 e 𝐶 , tanto più il

materiale si avvicina al regime incomprimibile dove v=0.5. Una variante

particolare della forma polinomiale presentata si ottiene per N=1:

1
𝜔 𝐶 𝐼𝑩𝟏 3 𝐶 𝐼𝑩𝟐 3 𝐽 1 II.24
𝐷
La forma così ottenuta, che prevede la calibrazione di soli tre parametri,
corrisponde al modello di Money-Rivlin. Nei materiali elastomerici questo
tipo di modello non è in grado di rappresentare il cambio di pendenza che
si ha tipicamente oltre una data soglia di rapporto di allungamento, come
qualitativamente mostrato in Figura 5. Il rapporto di allungamento che
individua approssimativamente il cambio di pendenza viene definito
“locking stretch” 𝜆

Figura 5 – Cambio di pendenza nella risposta dei

materiali elastomerici

18
Altre forme polinomiali semplificate sono le cosiddette forme polinomiali
ridotte, nelle quali è eliminata la dipendenza dal secondo invariante delle
deformazioni deviatoriche 𝐼𝑩𝟐 , ottenendo espressioni del tipo:

1
𝜔 𝐶 𝐼𝑩𝟏 3 𝐽 1 II.25
𝐷

L’adozione di forme ridotte è motivata dal confronto fra i modelli teorici

ed i comportamenti sperimentali degli elastomeri, che risultano abbastanza
ben interpolati da modelli che non prevedono grandi sensibilità al secondo
invariante 𝐼𝑩𝟐 . Si è anche osservato che l’introduzione dei termini in 𝐼𝑩𝟐 ,
pur consentendo una più accurata riproduzione del comportamento del
materiale nelle prove usate per la calibrazione, tende a produrre risultati
non soddisfacenti nella predizione del comportamento in condizioni
diverse da quelle di calibrazione. La forma polinomiale più semplice
utilizzata è la cosiddetta neo-Hookean form, che è rappresentata da una
approssimazione polinomiale ridotta del primo ordine

1
𝜔 𝐶 𝐼𝑩𝟏 3 𝐽 1 II.26
𝐷

 Modello di Ogden
Nel modello di Ogden il potenziale elastico è espresso attraverso dei
rapporti di allungamento principali normalizzati, ed ha la forma:

2𝜇 1
𝜔 𝜆̅ 𝜆̅ 𝜆̅ 3 𝐽 1 II.27
𝛼 𝐷

Gli allungamenti principali normalizzati 𝜆̅ , 𝜆̅ e 𝜆̅ sono definiti come:

1
𝜆̅ 𝜆 ⟹ 𝜆̅ 𝜆̅ 𝜆̅ 1 II.28
𝐽

19
In base all’eq. 28 lo stato di deformazione descritto dai rapporti di
allungamento principali normalizzati è caratterizzato da una deformazione
volumetrica nulla ed è pertanto puramente deviatorico. Il modello di
Ogden preserva quindi il disaccoppiamento fra la parte volumetrica e
deviatorica della deformazione, ma fornisce direttamente l’espressione del
potenziale in funzione dei rapporti principali di allungamento consentendo
una più immediata calibrazione del modello.

 Modello di Arruda-Boyce
Il modello di Arruda-Boyce è ricavato da un’idealizzazione della struttura
molecolare di un elastomero. Viene considerato un modello di volume
rappresentativo di materiale elastomerico costituito da un cubo che
presenta al suo interno 8 elementi elastici convergenti nel centro
dell’elemento, schematizzato in Figura 6.

Figura 6 – Elemento di volume rappresentativo nel

modello di Arruda-Boyce

Il potenziale nel modello di Arruda-Boyce ha la seguente forma:

𝐶 1 𝐽 1
𝜔 𝜇 𝐼𝑩𝟏 3 𝑙𝑛 𝐽 II.29
𝜆 𝐷 2
dove i coefficienti 𝐶 non rappresentano però parametri incogniti, ma
hanno valori ottenuti applicando un modello statistico alla distribuzione di
proprietà elastiche degli otto elementi inclusi nel volume rappresentativo.
In particolare, i coefficienti 𝐶 sono ottenuti dall’espansione in serie di una
funzione matematica, ed hanno espressione:

20
 1 1 11 19 519
𝐶 ; 𝐶 ; 𝐶 ; 𝐶 ; 𝐶
2 20 1050 7000 673750

I parametri da calibrare nel modello risultano i tre coefficienti , 𝜆 e D,

dove 𝜆 rappresenta il locking stretch introdotto in Figura 3. Si osservi
anche che la parte volumetrica della deformazione dipende dal solo
parametro D, che è direttamente esprimibile in funzione del bulk modulus
del materiale secondo l’espressione riportata in eq. 23.

II.5 Calibrazionedei modelli di potenziale elastico per

materiali elastomerici
I modelli di potenziale elastico presentati in precedenza sono stati elaborati
per rappresentare la risposta di materiali elastomerici isotropi con
comportamento elastico e non lineare. Nelle procedure di calibrazione i
parametri da cui dipende il potenziale elastico devono essere fissati
basandosi sulla correlazione fra le risposte sperimentali e le predizioni
teoriche in una serie di prove condotte in condizioni semplici di sforzo.
Nelle prove sperimentali considerate le direzioni principali di sforzo e
deformazioni sono coincidenti, per l’isotropia dei materiali, e note. Le
prove individuano direttamente le risposte in termini di curve di sforzo
nominale – deformazione nominale mentre l’identificazione
dell’andamento degli sforzi di Cauchy sarebbe complicata dalla necessità
di misura le notevoli variazioni di area resistente durante il processo
deformativo.
Pertanto, è significativa la possibilità di ricavare dai modelli di potenziale
il valore teorico previsto per gli sforzi nominali principali, nin
corrispondenza di uno stato di deformazione descritto dagli allungamenti
principali 𝜆 .
Tuttavia, anche per stati di sforzo semplici e controllati, quali quelli di una
prova uniassiale di trazione, non è possibile distinguere a priori la

21
deformazione deviatorica e volumetrica senza conoscere le caratteristiche
del materiale. Per questo motivo la separazione dei contributi deviatorici
e volumetrici nel potenziale elastico è di notevole aiuto nella procedura di
calibrazione. In particolare, la bassissima comprimibilità dei materiali
elastomerici è sfruttata per assumere che, alla presenza di uno stato di
sforzo in grado di indurre deformazioni deviatoriche e volumetriche nel
materiale, la deformazione volumetrica sia, in effetti, trascurabile. Ciò
comporta che il materiale possa essere considerato incomprimibile e che i
rapporti di allungamento principali rappresentino i valori principali della
deformazione deviatorica, introdotti in eq. 20. Sotto queste ipotesi, quindi:

𝐼𝑩𝟏 𝜆𝑩𝟏 𝜆𝑩𝟐 𝜆𝑩𝟑 𝜆 𝜆 𝜆

𝐼𝑩𝟐 𝜆𝑩𝟏 𝜆𝑩𝟐 𝜆𝑩𝟑 𝜆 𝜆 𝜆 II.30

𝐽 𝜆 𝜆 𝜆

Il potenziale elastico, espresso in funzione dei rapporti di allungamento

principali, risulta:

𝜔 𝜔 𝐼𝑩𝟏 𝜆 , 𝜆 , 𝜆 ; 𝐼𝑩𝟐 𝜆 , 𝜆 , 𝜆 ; 𝐽 II.31

Gli sforzi nominali principali si ottengono derivando il potenziale e,

ricordando che 𝑑𝐽 0 per le assunzioni introdotte sul comportamento del
materiale, si ottiene:

𝜕𝜔 𝜕𝜔 𝜕𝜔
𝛴 𝑑𝐼 𝑑𝐼 𝑑𝐽
𝜕𝐼𝑩𝟏 𝑩𝟏 𝜕𝐼𝑩𝟐 𝑩𝟐 𝜕𝐽
II.32
𝜕𝜔 𝜕𝐼𝑩𝟏 𝜕𝜔 𝜕𝐼𝑩𝟐
𝜕𝐼𝑩𝟏 𝜕𝜆 𝜕𝐼𝑩𝟐 𝜕𝜆
Gli invarianti deviatorici pertanto possono essere determinati in funzione
dei rapporti di allungamento principali e, grazie all’ipotesi di
incomprimibilità, possono essere calcolati quando è noto il rapporto di
allungamento in una direzione. L’eq. 32 consente quindi di esprimere il

22
valore teorico dello sforzo nominale 𝛴 , ottenuto applicando una
determinata forma del potenziale elastico ω.
La prova più semplice per la calibrazione del potenziale elastico nei
materiali elastomerici è rappresentata dalla prova uniassiale di trazione e
compressione, schematizzata in Figura 5. Il test di trazione è generalmente
eseguito su provini tipo dog-bone, mentre il test di compressione è
eseguito comprimendo un disco di materiale fra due piastre lubrificate in
modo da lasciare libera l’espansione delle direzioni trasversali al carico
applicato.

Figura 5 – Prova di trazione uniassiale (A) e di compressione

(B) per la calibrazione del potenziale elastico.

Nelle condizioni di prova considerate la prima direzione principale dello

stato di sforzo-deformazione corrisponde alla direzione di applicazione del
carico. Sotto le assunzioni di incomprimibilità si ha, in corrispondenza di
un determinato rapporto di allungamento nella direzione del carico
applicato, 𝜆 :

1
𝜆 𝜆 ; 𝜆 𝜆 II.33
𝜆
L’applicazione dell’Eq. 32 conduce alla seguente espressione:

𝜕𝜔 𝜕𝜔 𝜕𝜔
𝛴 2 1 𝜆 𝜆 II.34
𝜕𝜆 𝜕𝐼𝑩𝟏 𝜕𝐼𝑩𝟐

23
Una seconda prova di calibrazione è la prova bi-assiale, schematizzata in
Figura 6, nella quale un carico di uguale valore è applicato secondo due
direzioni ortogonali a un provino. La prova di trazione può essere eseguita
mediante macchine di prove biassiali mentre la prova di compressione è
di più difficile esecuzione ed è raramente applicata.

Figura 6 – Prova biassiale di trazione (A) e di

compressione (B per la calibrazione del potenziale
elastico

I rapporti di allungamento principali risultano, sotto le ipotesi considerate:

1
𝜆 𝜆 𝜆 ; 𝜆 II.35
𝜆
dove 𝜆 è il rapporto di allungamento in ciascuna delle due direzioni di
applicazione del carico. Derivando il potenziale elastico in funzione di 𝜆
si ottiene il doppio del valore dello sforzo principale nelle direzioni di
carico. Pertanto il valore teorico dello sforzo nominale risulta:

1 𝜕𝜔 𝜕𝜔 𝜕𝜔
𝛴 2 𝜆 𝜆 𝜆 II.36
2 𝜕𝜆 𝜕𝐼𝑩𝟏 𝜕𝐼𝑩𝟐

Infine, incollando dei provini sottili a due piastre collegate alla macchina
di prova, è possibile effettuare test in condizioni di deformazione piana,
rappresentati in Figura 7.

24
 Figura 7 – Prova in stato di tensione piana (A) e di
compressione (B) per la calibrazione del
potenziale elastico.

I rapporti di allungamento per questo stato di sforzo deformazione, sotto

le ipotesi di incomprimibilità, risultano:

1
𝜆 𝜆 ; 𝜆 1 ; 𝜆 II.37
𝜆
Lo sforzo normale nella direzione di applicazione del carico è ottenuto
dall’espressione seguente:

1 𝜕𝜔 𝜕𝜔 𝜕𝜔
𝛴 2 𝜆 𝜆 II.38
2 𝜕𝜆 𝜕𝐼𝑩𝟏 𝜕𝐼𝑩𝟐

Le prove presentate in precedenza consentono la calibrazione della parte

deviatorica del legame sforzo-deformazione. Le prove a disposizione
possono essere anche solo alcune fra quelle indicate e consentono di
ottenere un set di dati sperimentali 𝜆 - Σ che deve essere interpolato
dal modello iperelastico. Poiché sotto l’ipotesi di incomprimibilità la
sovrapposizione di uno stato di sforzo volumetrico non influenza i valori
del potenziale elastico, alcune prove sono in realtà equivalenti fra di loro
e forniscono identiche informazioni per la calibrazione del modello.
In particolare:
- la tensione uniassiale è equivalente alla compressione biassiale;
- la compressione uniassiale è equivalente alla tensione biassiale;

25
- la tensione in stato piano di deformazione è equivalente alla
compressione in stato piano di deformazione.
Una volta raccolti i dati sperimentali, la calibrazione dei parametri del
modello di potenziale può essere condotta minimizzando l’errore
complessivo fra i valori sperimentali e quelli teorici.

1 ∑ 𝜆
𝐸𝑅𝑅 II.39
∑ 𝜆

La parte del potenziale elastico che dipende dalla deformazione

volumetrica è calibrata separatamente sebbene, in molti casi, l’assunzione
di incomprimibilità è adeguata a modellare con sufficiente accuratezza la
risposta di molti materiali elastomerici. Quando la calibrazione della
risposta volumetrica è richiesta, occorre eseguire prove volumetriche che,
tipicamente, sono condotte confinando un blocco cilindro di materiale
elastomerico in un contenitore rigido e applicando alla base uno sforzo
attraverso un pistone collegato con la macchina di prova. La prova e lo
stato di sforzo ottenuto sono schematizzati in Figura 8 e Figura 9.

Figura 8 – Prova volumetrica con confinamento

del materiale in un contenitore rigido.

26
 Figura 9 – Stato di sforzo idrostatico di trazione
(A) e di compressione (B) per la calibrazione
del potenziale elastico

In un test volumetrico puro si ha:

𝜆 𝜆 𝜆 𝜆 ⟹𝐽 𝜆

𝐼𝑩 𝜆𝑩 𝜆𝑩 𝜆𝑩 1 1 1 3 II.40

𝐼𝑩 𝜆𝑩 𝜆𝑩 𝜆𝑩 1 1 1 3

27
 Capitolo III

Modello Viscoelastico
III.1 Introduzione

Nei solidi reali possono manifestarsi deviazioni più o meno sensibili dalla
perfetta elasticità nella risposta della deformazione alla tensione, per cui è
necessario discutere la possibile forma di relazioni di tipo non elastico tra
le componenti di tensione e le componenti di deformazione. In particolare
può cadere in difetto la corrispondenza monodroma tra le dette
componenti, perché una parte della deformazione non si annulla quando
vengono rimosse le forze applicate al solido: siamo allora in presenza di
deformazioni permanenti dovute alla irreversibilità del processo di cui
esse sono conseguenza. Un comportamento perfettamente elastico
esigerebbe infatti che l’energia potenziale acquista dal solido nel corso
della deformazione venisse totalmente restituita in modo da rispettare la
perfetta reversibilità della trasformazione nel senso della termodinamica.
Evidentemente tale circostanza non può verificarsi nella realtà dovendo
essere sempre presente una dissipazione di energia nel corso della
deformazione: ad essa è dovuto ogni scostamento più o meno sensibile dal
modello ideale del solido elastico. Ammessa dunque una certa
irreversibilità di trasformazione e la presenza di dissipazione di energia
elastica, si potrà parlare generalmente di risposta anelastica; sarà quindi
conveniente pensare a modelli ideali di solidi in cui questo tipo di risposta
segua leggi particolari dove, oltre al consueto legame elastico tra le
componenti di tensione e di deformazione, dovranno comparire ulteriori
relazioni indipendenti, atte a caratterizzare il fenomeno specifico causa di
anelasticità. D’altra parte in nessun materiale reale lo stato di

28
deformazione è unicamente funzione dello stato di tensione, ma dipende
dalla temperatura, dalla distribuzione più o meno ordinata a livelli
molecolari o atomici, dal campo elettrico e magnetico. La presenza di
queste nuove variabili induce, in generale, due fenomeni tipici: la
diffusione o rilassamento delle loro fluttuazioni, cioè la tendenza ad
assumere un valore uniforme nei punti del solido, e la variazione dei loro
valori per effetto della deformazione.
La teoria classica dell’elasticità studia il comportamento meccanico dei
solidi allo stato elastico; per stato elastico si intende quello stato in cui la
tensione σ dipende dalla deformazione ɛ ed è indipendentemente dalla
velocità di deformazione ɛ : se il legame è di tipo lineare allora si parla di
solido elastico lineare (𝜎 𝐸𝜀 legge di Hooke), se il legame è di tipo non
lineare allora si parla di materiali iperelastici (figura 10).

Solido elastico lineare Solido iperelastico

σ σ

ɛ ɛ

Figura 10

In realtà non sempre si può ammettere l’ipotesi di perfetta elasticità sopra

indicata; le diverse modalità con cui può essere applicato e rimosso il
carico, sia esso deformativo o tensionale, induce un comportamento nuovo
definito viscoelasticità. Nel grafico indicato in figura 11 si vede

29
qualitativamente come un generico materiale può mostrare curve
caratteristiche diverse a seconda della velocità di deformazione ɛ imposta.

σ ɛ1 > ɛ2 ɛ1 ɛ2

ɛ
Figura 11

III.2 Rilassamento e scorrimento viscoso

Esistono due metodi per determinare il comportamento viscoelastico dei
materiali: il test a rilassamento e il test a scorrimento viscoso (creep test).

 Test a rilassamento
Nel test a rilassamento (vedi figura 12) una deformazione monoassiale di
allungamento ɛ0 viene applicata in maniera “quasi istantanea” ad un
provino e viene mantenuta per un tempo abbastanza lungo; in serie al
provino è applicata una cella di carico che registra l’andamento temporale
della tensione.

30
 ɛ

ɛ
σ
ɛ
σ t0 t
risposta
σ

rilassamento di tensione

t0 t
Figura 12

Essendo la tensione dipendente dal tempo e la deformazione costante, si

potrà definire un modulo di rilassamento E(t) anch’esso variabile nel
tempo, definito come segue:

𝜎 𝑡
𝐸 𝑡 𝑚𝑜𝑑𝑢𝑙𝑜 𝑑𝑖 𝑟𝑖𝑙𝑎𝑠𝑠𝑎𝑚𝑒𝑛𝑡𝑜 III.1
𝜀
Per sviluppi successivi è bene definire il modulo elastico iniziale (t = t0) e
il modulo all’equilibrio (t = ∞):

𝜎 𝑡 0
𝐸 𝑡 𝑡 𝐸 𝑚𝑜𝑑𝑢𝑙𝑜 𝑖𝑛𝑖𝑧𝑖𝑎𝑙𝑒 III.2
𝜀

e
𝜎 𝑡 ∞
𝐸 𝑡 ∞ 𝐸 𝑚𝑜𝑑𝑢𝑙𝑜 𝑎𝑙𝑙 𝑒𝑞𝑢𝑖𝑙𝑖𝑏𝑟𝑖𝑜 III.3
𝜀

31
  Test a scorrimento viscoso (creep)
Un altro test fondamentale per caratterizzare la viscoelasticità di un
materiale è il test a scorrimento viscoso. Una tensione uniassiale costante
𝜎 viene applicata “quasi istantaneamente” ad un provino e mantenuta per
un certo tempo; contemporaneamente si misura l’allungamento
progressivo nel tempo ottenendo delle curve qualitativamente come quelle
indicate di seguito.

σ0

t0 t1 t
recupero elastico istantaneo

ɛ deformazione anelastica (creep)

ɛ0

t0 t1 t

deformazione elastica istantanea recupero anelastico Deformazione plastica

non recuperata

Figura 13

Si può definire un modulo di deformabilità a scorrimento viscoso D(t) nel

seguente modo:

32
 𝜀 𝑡
𝐷 𝑡 𝑚𝑜𝑑𝑢𝑙𝑜 𝑑𝑖 𝑑𝑒𝑓𝑜𝑟𝑚𝑎𝑏𝑖𝑙𝑖𝑡à III.4
𝜎

III.3 Modelli viscoelastici elementari: modello di Maxwell e

Kelvin
In questo paragrafo iniziamo ad introdurre i modelli elementari di solido
che possano, per quanto possibile, rappresentare il comportamento
viscoelastico nei suoi tratti più essenziali. Come già indicato
precedentemente, un materiale viscoelastico unisce due comportamenti
diversi: un comportamento elastico e un comportamento viscoso. Con il
termine “comportamento” si intende il legame costitutivo del materiale,
cioè come il materiale si tensiona a seguito di una deformazione indotta (o
viceversa), e trovare quindi il legame analitico tra tensione e
deformazione.
Nel caso di solido elastico lineare tale legame è rappresentato dalla legge
di Hooke:

𝜎 𝐸𝜀 III.5

dove E è il modulo di rigidezza (modulo di Young); tale modello viene

rappresentato schematicamente nel seguente modo:

E
1
ɛ
Figura 14

33
La deformazione elastica, per definizione, è una deformazione che viene
totalmente recuperata in un ciclo di carico-scarico.
Il comportamento viscoso lineare è rappresentato dalla legge di Newton

𝜎 𝜂𝜀 III.6

dove 𝜂 è il coefficiente di viscosità ed 𝜀 è lo strain rate; il modello è

rappresentato schematicamente nella figura 15:

ɳ
1

ɛ
Figura 15

A sinistra della figura 15 si rappresenta schematicamente un pistone

all’interno di un cilindro riempito di fluido viscoso; in virtù del moto del
pistone il fluido si trova costretto a defluire, e quindi a dissipare in calore
l’energia meccanica posseduta dal pistone. Funzionando puramente da
dissipatore, la deformazione indotta non viene integralmente recuperata,
come mostrato nelle figure seguenti.
σ deformazione anelastica
ɛ
non recuperata

σ0

t0 t t0 t

Figura 16

34
  Modello di Maxwell:
Il modello di Maxwell consta di una molla e un dissipatore disposti in
serie, come indicato nella figura 17

Figura 17

Riscrivo il legame costitutivo della molla e del dissipatore esplicitando la

deformazione:

𝜀 𝜎 ; 𝜀 𝜎 III.7

Essendo molla e dissipatore disposti in serie la deformazione totale

all’istante t è pari alla somma delle deformazioni dei singoli elementi,
mentre la tensione è costante e pari al valore σ, quindi

𝜀 𝜀 𝜀 III.8

Differenziando rispetto al tempo si ottiene:

𝜀 𝜀 𝜀 III.9

da cui si ottiene

𝜎 𝜎 𝜂𝜀 III.10

che rappresenta l’equazione costitutiva del modello di Maxwell.

35
Eseguendo un test a scorrimento viscoso, ovvero imponendo
analiticamente una funzione 𝜎 “a gradino” si ottine, sviluppando
l’equazione, la risposta in deformazione del tipo indicato in figura 18:

1 𝑡
𝜀 𝑡 𝜎 III.11
𝐸 𝜂

𝝈𝟎
𝑬

t
Figura 18

Il modulo di deformabilità risulta essere:

𝜀 𝑡 1 𝑡
𝐷 𝑡 III.12
𝜎 𝐸 𝜂
Allo stesso modo si procede per ottenere la risposta al test di rilassamento;
imponendo analiticamente una funzione 𝜀 “a gradino” si ottiene la
relazione seguente:

𝜎 𝑡 𝜀 𝐸𝑒 III.12

dove 𝜏 è definito tempo di rilassamento; l’andamento di σ(t) mostra

una caduta esponenziale della tensione, come mostrato qualitativamente

nel grafico che segue
σ0

t
Figura 19

36
Il modulo di rilassamento è definito:

𝜎 𝑡
𝐸 𝑡 𝐸𝑒 III.13
𝜀

 Modello di Kelvin:
Il modello di Kelvin è costituito da una molla e un dissipatore disposti in
parallelo, come indicato nella figura seguente:

Figura 20

La disposizione in parallelo degli elementi fa sì che entrambi sperimentino

la stessa deformazione ma sono soggetti a tensioni diverse, quindi valgono
le seguenti relazioni:

1 1
𝜀 𝜎 𝜀 𝜎 𝜎 𝜎 𝜎 III.14
𝐸 𝜂
da cui si ricava la legge costitutiva del modello di Kelvin:

𝜎 𝐸𝜀 𝜂𝜀 III.15

Imponendo una tensione 𝜎 “a gradino”, la risposta al creep è data dalla

relazione seguente:

𝜎
𝜀 𝑡 1 𝑒 III.16
𝐸

37
 ɛ
Graficando tale relazione si
𝝈𝟎
ottiene il grafico in figura 21
𝑬

t
Figura 21

Il modulo di deformabilità 𝐷 𝑡 assume la forma seguente:

1
𝐷 𝑡 1 𝑒 III.17
𝐸
dove per valori di t elevati il modulo di deformabilità 𝐷 𝑡 si stabilizza

asintoticamente al valore .
Come si può notare dal grafico, la risposta deformativa ad uno stimolo di
tensione non avviene in maniera istantanea; questo perché sebbene la
molla possa reagire istantaneamente a tale tensione, il dissipatore non ne
consente una risposta immediata, e quindi tende a prevalere sulla molla.
La risposta al rilassamento, ottenuta fornendo una deformazione 𝜀 “a
gradino”, non fornisce risultati interessanti. Applicare una deformazione
finita in un tempo nullo implica una velocità di deformazione infinita e ciò
comporterebbe una risposta di tensione infinita da parte del dissipatore.
Quindi il modello di Kelvin non ha modulo di rilassamento.

III.4 Modello di Maxwell generalizzato

I modelli elementari di Maxwell e Kelvin hanno dei limiti intrinseci di
rappresentazione costitutiva dei materiali viscoelastici reali. Il modello di
Maxwell, pur prevedendo una risposta elastica istantanea alla
sollecitazione, non prevede nessun recupero di deformazione anelastica

38
alla rimozione del carico. Il modello di Kelvin non mostra nessuna risposta
elastica deformativa nell’istante di applicazione del carico e la
deformazione viene totalmente recuperata dalla molla come deformazione
elastica. Per i motivi appena citati i modelli elementari di Maxwell e
Kelvin non vengono utilizzati per rappresentare il comportamento
viscoelastico; verrà invece discusso e utilizzato ai fini di questa tesi il
modello di Maxwell generalizzato.
Il modello di Maxwell generalizzato (vedi figura 22) è un modello
costituito da n rami di Maxwell disposti in parallelo che rappresentano la
parte viscoelastica, più una molla (in parallelo ai rami viscoelastici) di
modulo 𝑬 che rappresenta la parte iperelastica.

σ
componente elastica Componente viscoelastica

σ
Figura 22

L’equazione alla base del modello di Maxwell generalizzato, che coniuga

equilibrio di tensione e congruenza di deformazione, è un equazione
differenziale lineare a coeffienti costanti di ordine n non omogenea del
tipo
°
σ p σ p σ ⋯ p σ
III.18
°
q ɛ q ɛ ⋯ q ɛ

39
Si assume per ora che le molle e i dissipatori del modello siano elementi
lineari, pertanto i moduli elastici 𝐸 e i coefficienti di viscosità ɳ sono
costanti durante l’intervallo di deformazione; questo rende l’equazione
differenziale di tipo lineare. Pertanto vale il principio di sovrapposizione
degli effetti (Principio di Boltzman), per cui la somma delle cause e uguale
alla somma degli effetti. Imponendo al sistema una deformazione “a
gradino”

𝜀 𝑡 𝜀 𝐻 𝑡 III.19

con H(t) funzione di Heavyside, tale deformazione si applica (per

congruenza) su tutti i rami del modello. La risposta in tensione sarà la
somma della tensione 𝜎 di ogni ramo (vedi III.12)

𝜎 𝑡 𝜀 𝐸𝑒

dove 𝜏 è il tempo di di rilassamento del ramo i-esimo. La risultante

delle tensioni di tutti i rami fornisce

𝜎 𝑡 𝜀 𝐸 𝑒 𝐸 III.20

La quantità

𝐸 𝑡 𝐸 𝑒 𝐸

viene definito modulo di rilassamento globale. La sommatoria

𝐸 𝑒 III.21

40
è una serie esponenziale decrescente definita serie di Prony. Con la serie
di Prony si rappresenta lo spettro discreto costitutivo del materiale
viscoelastico e, in generale, essa viene utilizzata anche senza far
riferimento a un modello meccanico particolare.
Qualora si abbia un materiale nuovo a presunto comportamento
viscoelastico si esegue un test a rilassamento, ovvero si impone una
deformazione “a gradino”, si collega in serie al provino una cella di carico
e si misura il valore di tensione al trascorrere del tempo. Se il materiale è
un solido viscoelastico si ottiene una curva sperimentale della tensione
con andamento ad esponenziale decrescente che tende asintoticamente ad
un certo valore diverso da zero; se il materiale è un fluido viscoelastico la
curva tende asintoticamente a zero. Ottenuta la curva sperimentale,
mediante una procedura numerica si ottengono le coppie 𝐸 e 𝜏 della serie
di Prony. La figura che segue schematizza il risultato di “fitting” di una
curva sperimentale.

Figura 23

A questo punto è importante fare alcune osservazioni sul modulo di

rilassamento del modello di Maxwell generalizzato, che riscrivo di seguito

41
 𝐸 𝑡 𝐸 𝑒 𝐸 III.20

All’istante t=0, cioè nell’istante di applicazione del carico deformativo a

“gradino”, la velocità di deformazione ɛ è molto alta (tendente
all’infinito), quindi i dissipatori si comportano come connessioni rigide (σ
tendente all’∞) in virtù del loro legame costitutivo 𝜎 ɳɛ . Sempre
all’istante t=0 gli esponenziali della serie di Prony tendono a 1. Il modello
è quindi ben approssimabile da n molle disposte in parallelo senza
dissipatori, come indicato nella figura 24.

Figura 24

Si definisce rigidezza istantanea all’istante t=0 (“short term” o

“istantaneous stiffness”) la relazione:

𝐸 𝐸 𝐸 III.21

Possiamo adimensionalizzare tale relazione rispetto 𝐸 ottenendo:

𝐸 𝐸 𝐸
1 III.22
𝐸 𝐸 𝐸

e definisco i moduli di rigidezza relativi nel seguente modo

42
 𝐸 𝐸
𝑔 ; 𝑔 III.23
𝐸 𝐸

Il modulo di rilassamento globale assume la forma:

𝐸 𝑡 𝐸 𝑔 𝑔𝑒 III.24

𝑔 e 𝜏 rappresentano i parametri necessari alla modellazione viscoelastica

e, come vedremo in seguito, tali parametri sono implementabili su Abaqus.
Per valori di tempo elevati i dissipatori esplicano la loro azione dissipando
tutta l’energia elastica accumulata nelle molle 𝐸 , gli esponenziali della
serie di Prony tendono a zero e il sistema conserva solamente l’energia
elastica accumulata nella molla di modulo 𝐸 . Quindi nel lungo periodo
il modello possiede una rigidezza rappresentata solo da 𝐸 , al quale si da
la definizione di modulo elastico di lungo periodo o “long-term stiffness”.

43
 Capitolo IV

Metodo a elementi finiti e

software di simulazione
IV.1 Introduzione

L’analisi agli elementi finiti è una tecnica di simulazione numerica

utilizzata in ambito strutturale per permettere la semplificazione della
soluzione di strutture complesse. Nel seguente capitolo sono analizzati gli
aspetti riguardanti tale tecnica di analisi ed il metodo utilizzato, facendo
un breve excursus delle varie tipologie di simulazione che si possono
realizzare. In una prima parte del seguente capitolo si fornisce un quadro
generale dell’analisi agli elementi finiti, spiegando le modalità di
interfaccia solutore-utente; successivamente vengono fornite alcune
nozioni di carattere teorico sui vari aspetti che le simulazioni possono
affrontare, con buona approssimazione della realtà. Lo sviluppo attuale
della tecnologia ha portato alla realizzazione di molti solutori numerici
dalle diverse capacità, per varie tipologie di analisi che è possibile
effettuare. Questi permettono la soluzione di innumerevoli problemi, che
vanno dalle semplici analisi di tipo statico a quelle non lineari, in cui si
coinvolgono anche zone di contatto, fino ad arrivare ad analisi dinamiche
ad alte e basse velocità. In ambito strutturale molti solutori, oltre a
effettuare analisi accurate su materiali isotropi, hanno sviluppato nel
tempo la capacità di implementare anche i materiali compositi. Il crescente
utilizzo di tali materiali in ambito industriale, ha infatti reso necessario lo
sviluppo di software dedicati. Nella presente tesi è stato utilizzato il
software ABAQUS della Simulia, uno tra i solutori commerciali

44
maggiormente diffusi sul mercato. Questo tipo di software, insieme ad altri
diffusi nel settore ingegneristico, si suddividono in tre differenti fasi:
i. pre-process.
ii. process.
iii. post-process.
La prima fase consente di costruire la struttura, definire le proprietà dei
materiali coinvolti, assegnare a ciascuna parte della struttura il materiale
costituente, definire il carico e le condizioni al contorno. Ottenute tutte le
informazioni necessarie, il solutore calcola le matrici di rigidezza del
modello, nonché le forze interne ed esterne agenti in esso. Nella fase di
process, con le informazioni fornite, vengono risolte le equazioni di
equilibrio e calcolati gli spostamenti. L’ultima fase, infine, consente di
visualizzare graficamente le soluzioni ottenute nel process.

 Pre-Process
Per effettuare il calcolo della matrice di rigidezza, il software ha bisogno
di informazioni sul tipo di materiale e sulla geometria della struttura da
simulare. Quest’ultima può essere creata attraverso i solutori CAD2
implementati nel programma FEM stesso, ma comunemente, in presenza
di geometrie complesse, si predilige creare la struttura utilizzando un buon
software CAD dedicato e solo dopo implementarlo nei solutori FEM. Si
deve però precisare che il passaggio da un software CAD al FEM può
creare problemi sulle connessioni, sui contatti e sulla geometria stessa
della struttura. Pertanto la scelta di implementare o meno un file esterno
deve essere ragionata in base alla complessità della geometria stessa. Di
solito, ogni programma FE ha una libreria che contiene una varietà di
elementi. Questi vengono identificati in categorie: elementi trave, solido,
shell. Ogni solutore commerciale identifica le diverse tipologie con un
opportuno codice. Dalla figura 25 si può notare come uno stesso elemento

45
viene identificato tramite codici diversi per due dei solutori commerciali
più diffusi, ANSYS della Dassault ed ABAQUS della Simulia.

Figura 25

Per ogni categoria precedentemente indicata si ha la possibilità di scegliere

tra differenti opzioni: su un elemento solido planare, ad esempio,
un’opzione permette di compiere analisi con deformazioni e tensioni
piane. Definiti in modo opportuno gli elementi strutturali, a questi deve
essere associato un tipo di materiale. Inoltre, in base al tipo di analisi, le
proprietà del materiale introdotte, possono essere di tipo lineare, non
lineare, isotropiche o ortotropiche, ed inoltre dipendenti o indipendenti
dalla temperatura. Altre caratteristiche del materiale possono essere
introdotte qualora le analisi necessitassero di particolari informazioni,
come, ad esempio, la tensione di rottura nel caso di failure analysis. Una
volta introdotte tutte le parti costituenti della struttura da simulare, queste
devono essere discretizzate in modo da creare una serie di nodi in cui
andranno calcolate le soluzioni.
Tale discretizzazione prende il nome di mesh e può essere realizzata con
conformazioni e dimensioni differenti in base alla precisione di soluzione
ed ai tempi di simulazione che si vogliono ottenere. Sarà compito
dell’utente trovare il giusto compromesso tra i due parametri. In sintesi, la
prima parte del pre-process, consiste nel creare la geometria del modello,

46
assegnare il materiale costituente ed i tipi di elementi che andranno a
suddividere il corpo in un insieme discreto di nodi. La seconda parte del
pre-process consiste nell’attribuire al corpo tutte le condizioni necessarie
affinché il solutore effettui i calcoli. Si indicano, a tal proposito la
definizione degli stati di carico, delle condizioni al contorno e dei vincoli.
Le condizioni al contorno sono i valori conosciuti dei gradi di libertà del
corpo sul bordo. In un’analisi strutturale i gradi di libertà consistono
essenzialmente in spostamenti e rotazioni. Esempi di condizioni al
contorno possono essere le cerniere o l’incastro. Ulteriori informazioni
necessarie al solutore sono le eventuali zone di contatto presenti, con i
rispettivi attriti; se si conoscono le superfici che interagiscono durante
l’analisi è possibile utilizzare dei contatti di tipo pair (contatto in cui viene
specificata la coppia di superfici interessate), altrimenti si utilizza il
general contact, in cui il solutore considera tutte le superfici come se
fossero coinvolte nel contatto. Il carico è un altro aspetto fondamentale per
l’analisi. Esso può essere impostato come forza concentrata o distribuita,
oppure come carico di taglio agente su una superficie. Questa seconda
parte del pre-process può essere vista come un’intersezione pre-process e
process, perché tramite i dati forniti in questa fase, in quella successiva il
solutore procede alla realizzazione dei calcoli numerici.

 Process
La fase del process è eseguita solo dopo una fase detta di check delle
informazioni inserite. Una volta terminata questa fase, parte il processo di
soluzione numerica, in cui attraverso l’uso del metodo agli elementi finiti
vengono risolte le equazioni algebriche. I tipi di analisi che si possono
effettuare sono vari, e vanno dalle analisi statiche lineari a quelle non
lineari, per finire con le analisi di tipo dinamico.

 Post- process

47
Nella fase di post-process si analizzano i risultati ottenuti dalla soluzione
che è stata calcolata. Tali risultati vengono, usualmente, visualizzati e
analizzati attraverso un’interfaccia grafica, in cui il modello appare
colorato: ad ogni colore è assegnata una fascia di valori numerici di
tensione o deformazioni, che consente di capire quali sono le zone
maggiormente sollecitate.

IV.2 Metodo a elementi finiti FEM

Le analisi agli elementi finiti sfruttano la teoria degli spostamenti per la
soluzione di strutture continue, che deriva dalla più complessa teoria delle
forze. Una rappresentazione reale, in condizioni statiche, di un corpo sotto
l’azione di un carico può essere ottenuta mediante una equazione
differenziale. Essendo, però, la sua risoluzione estremamente complessa,
se non impossibile da ottenere in forma chiusa, si ricorre alla teoria degli
spostamenti, che permette di approssimare la soluzione di un’equazione
differenziale con quella di una equazione algebrica. Tale teoria prevede la
discretizzazione della struttura in elementi, i quali permettono lo scambio
di sforzi solo tra i nodi; le informazioni sulle interfacce, quindi, non
vengono considerate. Ciò implica che ogni elemento deve garantire
l’equilibrio sui nodi. Sebbene nella realtà non esistano spostamenti nodali,
nel metodo descritto sono previsti: imporre l’equilibrio e la congruenza sui
nodi è, pertanto, condizione necessaria per ottenere la soluzione. Le forze
che agiscono su ogni nodo si valutano sulla base della legge di Hooke:

P k∗u IV.1

dove con u si indica lo spostamento nodale e con k la rigidezza di un

materiale. Come noto, la legge di Hooke è valida solo nel campo di
linearità di un materiale; la non linearità verrà discussa inseguito. Per
rendere più chiari i concetti appena esposti, ipotizziamo un caso semplice:

48
una trave incastrata, approssimata con due elementi consecutivi. Ogni
elemento possiede due nodi; si ottiene quindi una struttura formata da due
elementi e tre nodi, di cui uno dei nodi esterni vincolato, quindi con
spostamento nullo (figura 26). La rigidezza del materiale viene
approssimata tramite un elemento molla: vale cosi la legge di Hooke su
ogni nodo.

Figura 26

Imponendo l’equilibrio sui vari nodi, si otterrà un sistema matriciale in cui

le incognite risulteranno essere gli spostamenti nodali, e la cui soluzione
richiederà l’inversione della matrice di rigidezza. Dall’equilibrio delle
forze sui nodi si ottiene, infatti:

𝑃 𝑘 𝑘 0 𝑢
𝑃 𝑘 𝑘 𝑘 𝑘 ∗ 𝑢 IV.2
𝑃 0 𝑘 𝑘 𝑢

dove 𝑢 , 𝑢 e 𝑢 sono gli spostamenti dei nodi in figura 26, 𝑘 e 𝑘 sono

le rigidezze degli elementi A e B, mentre 𝑃 , 𝑃 e 𝑃 sono le forze che
agiscono sugli elementi A e B. Condizione necessaria per l’inversione
della matrice di rigidezza, e quindi per l’ottenimento degli spostamenti,
consiste nella non singolarità della matrice stessa, ottenibile con
l’imposizione di un vincolo nella struttura. Se si considera invece una
struttura continua, nonostante il principio sia molto simile a quello
descritto, intervengono effetti che complicano il processo di calcolo.
Innanzitutto, in un sistema continuo è necessario considerare il sistema di
riferimento fisso rispetto a quello relativo locale; per tener conto di ciò si
introduce una matrice di rotazione per ogni elemento, che individua la

49
posizione del sistema di riferimento locale, rispetto a quello fisso. Oltre
allo stato tensoriale sui nodi, occorre definire anche le tensioni agenti
sull’intero corpo. Per fare questo, si utilizzano le funzioni di forma che
permettono di correlare gli spostamenti nodali agli spostamenti interni.
Queste assicurano la congruenza, ma non l’equilibrio delle tensioni
dell’intero corpo. Dalle funzioni di forma si ricavano le deformazioni per
le quali è garantita, pertanto, solo la congruenza. Se si ipotizza che lo stato
di deformazione congruente equivale allo stato di deformazione
equilibrato e congruente (per cui reale), e facendo uso del principio dei
lavori virtuali, si ricava la matrice di rigidezza approssimata dell’intero
corpo. Conoscendo la matrice di rigidezza è possibile ricavare lo stato
tensoriale approssimato dell’intero corpo. Pertanto la determinazione della
matrice di rigidezza deriva dalle funzioni di forma. L’ipotesi introdotta è
fondamentale ed è alla base del metodo degli spostamenti, che permette di
ottenere risultati che approssimano bene il caso reale.
Da quanto esposto, si evince che risulta importante definire nel miglior dei
modi la discretizzazione del corpo, cioè definire la giusta quantità di
elementi in cui suddividere la struttura. Infatti più gli elementi sono di
piccola entità, maggiore sarà il numero di nodi totali posseduti dall’intero
corpo: il solutore numerico dovrà effettuare i calcoli su ogni singolo nodo
e, quindi, aumenterà la complessità di calcolo del problema. Sarà l’abilità
dell’utente a permettere di ottenere un buon compromesso tra risoluzione
della discretizzazione e complessità computazionale.

IV.3 Analisi lineari e non lineari

Nel corso degli ultimi decenni le analisi FE sono diventate uno strumento
indispensabile per la progettazione di strutture complesse. L’utilizzo
sempre maggiore di tecniche FEA ha contribuito al perfezionamento di
metodologie di soluzione e alla semplificazione dell’interfaccia CAE-

50
FEM. L’analisi agli elementi finiti deve essere supportata da una
approfondita conoscenza della fenomenologia alla base del problema da
affrontare, affinché possa essere utilizzata nel modo corretto limitando i
possibili errori. Ad esempio, un’errata definizione delle condizioni al
contorno o dei vincoli, possono causare divergenze numeriche o risultati
non veritieri. I modelli maggiormente utilizzati a supporto della fase di
progettazione sono quelli lineari, che permettono una maggiore rapidità di
soluzione. Non sempre però è sufficiente utilizzare un modello lineare; in
questi casi è necessario ricorrere a modelli di tipo non lineare, che
permettono di ottenere soluzioni con migliore approssimazione rispetto al
caso reale. Il tempo necessario ad ottenere la soluzione numerica di un
problema di tipo non lineare, è di gran lunga più elevato che nel caso
lineare. La differenza sostanziale tra le simulazioni lineari e non, risiede
nelle ipotesi che vengono fatte sulla rigidezza del materiale, che
ricordiamo è l’attitudine di un corpo a resistere alle sollecitazioni con
piccole deformazioni.
Le principali caratteristiche che influiscono sulla rigidezza del corpo sono:

 Forma: ad esempio una trave a C ed una a T, hanno differenti

rigidezze;
 Tipo di materiale costituente: ogni materiale ha una rigidezza
propria;
 Vincoli applicati: una trave vincolata in modi differenti, subisce
differenti deformazioni (vedi figura 27).

51
 Figura 27

Un elemento sottoposto ad un carico può subire una variazione di forma,

per cui la sua rigidezza può cambiare a causa di uno o più dei fattori
precedentemente elencati, fino a raggiungere i limiti di cedimento per cui
le proprietà del materiale cambiano. Ad esempio, un qualsiasi corpo
sottoposto a trazione si deforma; al crescere della deformazione, però,
aumenta la forza necessaria da applicare per riuscire a deformarlo
ulteriormente, poiché aumenta la resistenza alla deformazione del corpo.
Se si assume che le variazioni di rigidezza siano sufficientemente piccole
da non influire in maniera rilevante sulle proprietà del materiale, è
possibile considerare l’analisi come lineare. Questa è l’ipotesi alla base
delle analisi lineari. Facendo riferimento quindi alla rigidezza di un
materiale, entro i limiti di elasticità verrà utilizzato un modello non lineare,
mentre al di fuori dei limiti di elasticità dovrà essere usato un modello non
lineare.
In definitiva, le analisi FE utilizzano il metodo degli spostamenti per la
soluzione numerica del problema, ovvero seguono la legge:

𝐹 𝑘 ∗ 𝑢 IV.3

in cui [F] indica il vettore delle forze esterne applicate note, [k] la matrice
di rigidezza del materiale nota ed [u] la matrice degli spostamenti nodali

52
non noti. Ipotizzando [k] costante si imposta un’analisi lineare in cui le
equazioni sono risolte in blocco, se invece la rigidezza non può essere
considerata costante, bisognerà ricorrere ad un’analisi non lineare e la
matrice dovrà essere continuamente aggiornata, aumentando il tempo di
calcolo. Questa è la causa della lentezza computazionale delle analisi non
lineari. In base ai fattori (precedentemente elencati) che influenzano la
rigidezza, si hanno comportamenti non lineari diversi. Ci sono altri casi
in cui sarà necessario far ricorso a modelli non lineari: ad esempio il caso
dall’instabilità dovuta all’effetto di carichi di punta. Nel caso di un carico
a trazione la rigidezza del materiale, infatti, aumenta mentre se si ha
compressione, e il corpo giunge ad instabilità, si possono avere due
possibili scenari: cedimento o nuova rigidezza strutturale.
L’analisi lineare è utile per il calcolo del valore del carico critico di una
struttura caricata con carico di punta; in alcuni casi tale carico può essere
maggiore del caso fisico reale: pertanto i risultati devono essere utilizzati
con cautela. Una volta raggiunto il carico di punta, non è detto che esso
causi una rottura catastrofica. È, quindi, necessaria una analisi non lineare
per prevedere il comportamento della struttura dopo il raggiungimento del
carico critico. Anche il contatto tra superfici può richiedere l’uso di analisi
non lineari, perché nella zona di contatto ci sono delle variazioni di
rigidezza (anche per zone di contatto non molto estese). Infine un altro
motivo che può rendere necessaria un’analisi non lineare è la presenza di
vincoli non lineari

Figura 28

53
La soluzione di tali metodi non lineari, sfrutta un codice di calcolo
numerico di tipo iterativo, che prende il nome di “metodo della rigidezza
variabile”. Spesso l’approssimazione lineare è ritenuta accettabile e
fornisce utili informazioni sul comportamento del corpo; in altri casi, un
modello lineare non è sufficiente a descrivere il comportamento reale di
una struttura che può essere meglio approssimato da un modello non
lineare. In generale, quindi, la scelta del tipo di analisi da effettuare,
dipende dalla necessità del caso e sarà frutto di un compromesso tra la
qualità dei risultati in termini di previsione del comportamento reale della
struttura e la complessità di calcolo.

IV.4 Algoritmo esplicito

Le analisi agli elementi finiti possono essere di tipo statico (descritte nei
paragrafi precedenti) o di tipo dinamico; in questi ultimi si tiene conto
degli effetti di inerzia, dei carichi dipendenti dal tempo e degli
smorzamenti. Anche l’analisi dinamica può essere di tipo lineare o non
lineare, e le regole che determinano il loro diverso utilizzo sono le stesse
discusse nel precedente paragrafo. Se la rigidezza del modello non muta
significativamente sotto il carico applicato, è sufficiente un approccio
analitico dinamico di tipo lineare. Ad esempio, un motore che vibra
subisce piccole deformazioni al punto di equilibrio. Altri problemi come
la simulazione di un impatto, richiede analisi dinamiche non lineari, a
causa delle elevate deformazioni in gioco. In questo lavoro di tesi, sono
stati costruiti dei modelli, sui quali sono state eseguite delle analisi
dinamiche, che meglio rappresentano il test di compressione che si vuole
simulare. Come anticipato, le analisi dinamiche sono utilizzate per
modellare prove di impatto o test in cui le forze siano dipendenti dal
tempo, permettendo una migliore comprensione dei fenomeni fisici che le
caratterizzano. Le analisi di tipo dinamico sono di grande complessità, ma

54
permettono di ottenere dei risultati sullo stato di sforzo tridimensionale
coerenti con la fenomenologia reale. Una corretta modellazione, quindi,
permette di ottenere risultati, che consentono di approssimare bene il caso
reale. Di contro, la dipendenza dal tempo dei modelli dinamici, può
causare, nel caso di errori nell’implementazione di vincoli o contatti,
eccessive non linearità ed elevato tempo computazionale. I metodi di
risoluzione delle equazioni differenziali che governano la meccanica
strutturale, si suddividono in metodi che utilizzano codici impliciti od
espliciti. I codici impliciti utilizzano il metodo di Newmark, il quale, però,
non permette di ottenere informazioni sulla rottura del materiale. Non è,
infatti, possibile effettuare l’analisi del danno, perché nessun criterio di
rottura è implementato nei solutori impliciti. I metodi espliciti, invece,
consentono di ottenere informazioni sulla rottura del materiale. A tal fine
si è fatto uso di un solutore esplicito, vista la necessità di rilevare lo stato
tensoriale a rottura del materiale, durante un test di compressione. La
differenza sostanziale tra i due metodi consiste nel fatto che quelli espliciti
determinano gli spostamenti al passo successivo, imponendo l’equilibrio
al passo precedente, mentre, i metodi impliciti impongono l’equilibrio ad
ogni passo. Ne deriva che ad ogni step di calcolo il metodo esplicito risulta
essere computazionalmente più rapido, ma richiede più intervalli
temporali per avere delle soluzioni accettabili. Gli impliciti richiedono
meno passi temporali, ma per ognuno è necessaria l’inversione della
matrice di rigidezza dinamica.
Il modello utilizzato da Abaqus/Explicit usa un metodo alle differenze
centrali per integrare nel tempo le equazioni del moto, partendo dalle
condizioni cinematiche (accelerazioni sugli spostamenti) al tempo i per
ottenere quelle al tempo i+1. Allo step iniziale il programma calcola
l’equilibrio dinamico utilizzando la formula derivante dal secondo
principio della dinamica:

55
 𝐹 𝐼 𝑀∗𝑢 IV.4

dove [F] è il vettore delle forze esterne applicate, [I] è il vettore delle forze
interne, [M] la matrice diagonale della massa e [𝑢] è il vettore delle
accelerazioni. Le incognite del problema saranno le accelerazioni nodali,
per cui è necessario effettuare l’inversione di tale equazione matriciale.
Tutti i termini dell’equazione precedente sono considerati all’i-esimo step.
Se indichiamo con t il tempo corrispondente all’i-esimo incremento si
possono calcolare le accelerazioni all’inizio dell’intervallo:

𝑢 𝑡 𝐹 𝐼 𝑡 ∗ 𝑀 IV.5

Questa equazione esprime l’equilibrio dinamico del sistema. Si noti che

l’accelerazione su ciascun nodo è determinata dalla massa che ad esso
compete e dalle forze che agiscono su di esso.
La velocità viene ottenuta come integrazione della accelerazione che è
considerata costante. Il metodo alle differenze centrali è considerato
esplicito perché utilizza gli stati noti 𝑢 𝑡 , accelerazione al tempo t, e
∆ ∆
𝑢 𝑡 velocità al tempo 𝑡 . Pertanto la velocità allo stato 𝑡
∆
sarà data dalla somma della velocità alla metà dell’incremento
∆
temporale precedente 𝑡 , che è noto, e della variazione temporale per

l’accelerazione al tempo t. Si effettua, cioè, una interpolazione lineare

delle velocità medie secondo l’equazione seguente:

∆𝑡 ∆ ∆𝑡
𝑢 ∆ 𝑢 ∆ ∗𝑢 IV.6
2

Le velocità così calcolate vengono integrate, ottenendo gli spostamenti.

Inserendole nell’equazione seguente e sommandole agli spostamenti

56
nodali ottenuti nell’intervallo precedente t, si ottengono gli spostamenti
all’intervallo 𝑡 ∆𝑡 :

𝑢 ∆ 𝑢 ∆𝑡 ∆ ∗𝑢 ∆
IV.7

Per come è strutturato, questo metodo integra solo le accelerazioni

costanti: per non avere errori numerici si devono considerare intervalli di
tempo non elevati. Lo schema di calcolo esplicito descritto può essere
riassunto nella figura 29:

Figura 29

Il punto di partenza di questo metodo numerico è rappresentato dalle

condizioni iniziali. Se si considera l’intervallo di tempo iniziale ∆𝑡, per
∆
come è stato definito il metodo, sarà necessario specificare 𝑢 , che

viene considerato dal programma pari a:

∆𝑡 ∆𝑡 0
𝑢 𝑢0 ∗𝑢 0 IV.8
2 2

Una volta calcolati gli spostamenti, ad ogni incremento si valutano le

deformazioni e successivamente le tensioni, che determinano le forze
interne.

57
IV.5 Tempo di incremento stabile
Il tempo di incremento stabile determina la stabilità del metodo. Dal punto
di vista fenomenologico, si può dire che la stabilità del processo di
integrazione diretta implica che gli errori numerici non si amplificano ad
ogni passo di integrazione. Una procedura viene definita
“incondizionatamente stabile” se la soluzione particolare dell’equazione
differenziale rimane limitata nel tempo per qualsiasi condizione iniziale.
Invece si parlerà di processo “condizionatamente stabile” quando la
soluzione è limitata nel tempo per Δ𝑡 di integrazione inferiori a valori
critici oltre i quali si può avere instabilità numerica. Non si deve
confondere la stabilità con l’accuratezza del metodo: il primo parametro
consente di definire la divergenza numerica del metodo, mentre il secondo
parametro indica l’errore che si commette rispetto ad un valore di
riferimento. Inoltre, il tempo di incremento stabile quantifica la durata di
una simulazione. Vista l’importanza di tale parametro è necessario
conoscere cosa lo influenzi, così da poter ottimizzare la durata della
simulazione. Il modello utilizzato da Abaqus/Explicit è un solutore alle
differenze centrali condizionatamente stabile. Il limite di stabilità, per
questa classe di metodi risulta di difficile determinazione, ma, essendo un
parametro fondamentale, se ne effettua una stima conservativa. Si può
dimostrare che per un sistema senza smorzamento, il limite di stabilità
sarà:

2
Δt IV.9
ω
dove 𝜔 è la massima frequenza naturale del corpo.
Invece, nel caso di un sistema smorzato, risulterà essere

2
Δt 1 𝜉 𝜉 IV.10
ω

58
dove 𝜔 è la frequenza naturale massima del corpo e ξ lo smorzamento
della struttura.
La massima frequenza naturale del sistema non è calcolabile in modo
esatto, poiché influenzata da molti fattori; se ne valuta, quindi, una stima
approssimativa. Si calcola la 𝜔 del singolo elemento, che si può
dimostrare essere maggiore della più alta frequenza dell’intero modello.
Si definisce, allora, il limite di stabilità come nella seguente relazione:

𝐿
∆𝑡 𝑚𝑖𝑛 IV.11
𝑐

dove 𝐿 è la lunghezza caratteristica dell’elemento di discretizzazione e

𝑐 la velocità di propagazione delle onde dilatazionali del modello. Si nota
che maggiori sono le dimensioni degli elementi della mesh tanto maggiore
sarà il tempo di incremento stabile. Inoltre, tale parametro aumenta al
diminuire della velocità di propagazione delle onde dilatazionali, che a sua
volta è inversamente proporzionale alla rigidezza specifica del materiale.
Conoscendo quindi le dimensioni del più piccolo elemento in cui è stata
discretizzata la struttura e la velocità di propagazione delle onde
dilatazionali del materiale, si può stimare il tempo di incremento stabile.
Il solutore ha la capacità di controllare la discretizzazione imposta,
manualmente o utilizzando le opzioni di default in modo da verificare
preventivamente se possono incorrere in instabilità durante l’analisi.
Anche il materiale può influire sul tempo di incremento; infatti, in base
alle sue caratteristiche, varia la velocità di propagazione delle onde
dilatazionali. Se il materiale è nella fase di comportamento elastico, ad
esempio, la velocità di propagazione è costante, e quindi il tempo di
incremento stabile dipenderà dal valore di rigidezza maggiore presente
nella struttura; se invece, si entra in campo plastico, la rigidezza e la
velocità di propagazione diminuiscono aumentando il limite di stabilità.

59
 Capitolo V

Implementazione del modello e

risultati delle analisi FEM
V.1 Introduzione

La tesi qui sviluppata si inserisce all’interno di una collaborazione fra i

gruppi di ricerca di costruzioni di macchine delle facoltà di ingegneria
dell’Università Politecnica delle Marche e dell’Università La Sapienza.
Gli aspetti generali della ricerca sono contenuti nell’articolo AIAS 2020–
1496 intitolato “Modellazione Costitutiva e Analisi Numerico-
Sperimentale del comportamento a impatto di schiume di PVC”, al
quale si è fatto ampio riferimento. Occorre dire che PVC e sughero sono
sì materiali diversi, ma il loro comportamento è descritto dallo stesso
modello costitutivo: modello hyperfoam, che rappresenta un modello di
Ogden modificato (vedi capitolo 2). Come già accennato, lo scopo della
tesi è la validazione di parametri costitutivi esistenti del materiale
mediante il confronto tra le curve velocità-tempo ottenute
sperimentalmente dalle prove di penetrazione a impatto e quelle prodotte
numericamente mediante Abaqus. Si sono eseguite simulazioni per tre
velocità di impatto diverse e, modellando opportunamente il materiale, si
è poi verificata l’indipendenza dei suddetti parametri dalla velocità di
impatto.

V.2 Macchina per test di perforazione

Le prove di perforazione a impatto vengono eseguite su campioni circolari
di sughero del diametro di 50 mm e 15 mm di spessore, in apposita torre

60
di caduta (modello Instron/Ceast9340) dotata di sistema di acquisizione
dati. Il sistema di acquisizione permette di registrare la forza di reazione
del materiale sull’impattatore durante l’urto, e di misurare la velocità
iniziale dell’impattatore stesso grazie ad un sensore ottico posto appena
sopra l’area di impatto.

Figura 30

La macchina consta essenzialmente di un penetratore a testa emisferica

(figura 30A) e di un supporto a base circolare (figura 30B).
Dalla registrazione della forza viene estrapolato l’andamento temporale
della velocità, dello spostamento e dell’energia attraverso un processo di
integrazione eseguito direttamente dal sistema. I campioni vengono posti
su un supporto di base circolare cavo con un diametro interno di 40 mm e
vengono fissati ad esso grazie ad un sistema pneumatico. Il penetratore
(figura 30A) ha un diametro di 12.7 mm e una massa di 3.055 kg. I
campioni di sughero vengono testati a 1.0 J, 2.5 J e 5 J, corrispondenti
rispettivamente alle velocità di impatto di 0,79 𝑚⁄𝑠 , 1,27 𝑚⁄𝑠 e 1,80
𝑚⁄𝑠, come già precedentemente indicato. Velocità ed energie di impatto
crescenti vengono ottenute aumentando l’altezza iniziale dell’impattatore.

61
La piena penetrazione dei campioni è stata registrata per la velocità di
impatto di 1,80 𝑚⁄𝑠.

V.3 Implementazione del modello iperelastico e

viscoelastico su Abaqus
 Comportamento iperelastico
Nello studio dei materiali solidi una delle prime prove che viene eseguita
è la prova uniassiale di trazione o compressione. Tipicamente, i materiali
cellulari come il sughero sopportano deformazioni fino 90% e trovano
applicazione laddove i carichi esterni producono stati tensionali di
compressione (vedi infatti le applicazioni nel settore packaging). Va da sé,
quindi, che la prova più importante per la caratterizzazione statica è la
prova di compressione.

Figura 31

Il grafico in figura 31 rappresenta un tipico andamento del legame

tensione-deformazione in una prova di compressione quasi statica, cioè
per basse velocità di deformazione 𝜀. Tale curva, nella quale il materiale
mostra deformazioni molto elevate, rappresenta un tipico comportamento
iperelastico ed è percorsa in fase di carico e in fase di scarico. Per materiali

62
cellulari si possono distinguere 3 fasi diverse: una prima fase per piccole
deformazioni a comportamento puramente elastico in cui avviene una
flessione delle pareti delle celle; una seconda fase cosiddetta di “plateau”
nella quale inizia il collasso delle pareti delle celle di materiale e la
deformazione è indipendente dalla tensione. Infine la fase di
densificazione delle celle di materiale nella quale la tensione riprende ad
aumentare.
Il comportamento iperelastico, in generale, viene agevolmente descritto
modellando la funzione potenziale di deformazione, come descritto nel
capitolo 2. Una variante al modello di Ogden, valida per materiali isotropi
non lineari e ad elevata porosità, è il modello hyperfoam. Il potenziale di
deformazione di tale modello è espresso dalla relazione

2𝜇 1
𝜔 𝜆 𝜆 𝜆 3 𝐽 1 V.1
𝛼 𝛽

dove ricordo che:

 N (numero intero) rappresenta l’ordine del polinomio;
 𝐽 rappresenta lo Jacobiano di deformazione;
 𝜆 sono i coefficienti di stretch principali nelle direzioni 1,2 e3;
 𝛼 e 𝛽 sono parametri utili al fitting della curva reale, il loro valore
è scelto in modo da avere la migliore approssimazione possibile.
Tali parametri mostrano una dipendenza dalla temperatura, in
particolare per ogni termine della funzione energia, il coefficiente
𝛽 determina il grado di compressibilità ed è correlato al
coefficiente di Poisson dalla relazione:

𝜈
𝛽 V.2
1 2𝜈

63
I parametri 𝜇 sono correlati al modulo di taglio iniziale 𝜇 dalla seguente
relazione

𝜇 𝜇 V.3

mentre i parametri 𝜇 e 𝛽 sono correlati al modulo di comprimibiltà

iniziale 𝐾 in base alla seguente relazione

1
𝐾 2𝜇 𝛽 V.3
3

La figura 32 seguente rappresenta la scheda di compilazione dei parametri

hyperfoam su Abaqus. I valori 𝜈 0 attestano la comprimibilità del
sughero, come risulta dall’evidenza sperimentale.

Figura 32

 Comportamento viscoelastico
Come descritto nel terzo capitolo, il comportamento meccanico di un
materiale generico è sempre intermedio tra un comportamento
perfettamente elastico ed uno perfettamente viscoso. Tale comportamento

64
viscoelastico lo si riscontra anche nel sughero; se si immagina di
schiacciare molto velocemente un provino di sughero e di rilasciare la
forza, si nota che il recupero alla forma originale avviene in un determinato
tempo. La viscoelasticità, essendo una proprietà (intrinseca) del materiale
e non attinente a un certo tipo di sollecitazione, la si riscontra anche nel
caso della prova di penetrazione.

Figura 33

65
La figura 33 mostra l’inserimento dei coefficienti della serie di Prony su
Abaqus; è stato considerato un modello di Maxwell generalizzato a 10
rami viscoelastici.

V.4 Effetto Mullins

Allo scopo di modellare correttamente la dissipazione di energia nel
sughero, che non è solo dovuta alla viscoelasticità, si deve introdurre un
nuovo concetto. Ricordo che stiamo modellando un provino di sughero
soggetto all’urto di un corpo penetrante e, nel caso di bassa velocità di
impatto, il provino non si rompe e fa rimbalzare il penetratore. Ricordo,
inoltre, che gli aspetti teorici fin qui sviluppati sono relativi a basse
velocità di deformazione 𝜀. Il provino di sughero subisce un ciclo di carico
(punzone che sta rilasciando energia cinetica e il provino che sta
accumulando energia potenziale elastica) e un ciclo di scarico (l’energia
di deformazione elastica rilasciata dal provino e recuperata dal punzone
sotto forma di energia cinetica). Questo trasferimento di energia, per
rispetto del secondo principio della termodinamica, avviene a meno di una
certa dissipazione.
Si consideri una prova di compressione quasi statica su un provino di
sughero; si esegue la lettura degli strumenti sia in fase di carico che in fase
di scarico del provino e si impone una storia di carico (deformazione-
tempo o tensione-tempo) del seguente tipo (figura 34)
𝜀
c’

c
b’ b’
C
b B B

a Figura 34
𝑡

66
l’andamento della curva tensione-deformazione è indicato
qualitativamente nella figura 35 seguente

Figura 35

Se consideriamo il primo ciclo di carico-scarico (a-b-b’-B) notiamo che le

curve non collimano, ovvero racchiudono un’area (zona colorata grigio)
che è indicativa di una avvenuta dissipazione di energia; l’area sottesa alla
curva (zona bianca) rappresenta l’energia recuperabile. Si può notare che,
a parità di deformazione raggiunta, la tensione al generico ciclo è inferiore
a quella del ciclo precedente; anche la rigidezza del materiale (data dalla
pendenza locale della curva) è inferiore. Questo fenomeno va sotto il nome
di stress softening o “Effetto Mullins”. Si può notare, inoltre, che al
generico ciclo di carico la curva collima con la curva iperelastica
“primaria” solo al superamento della deformazione massima raggiunta al
ciclo precedente.

 Aspetti teorici
I materiali iperelastici, come indicato nel secondo capitolo, sono descritti
in termini della funzione energia di potenziale di deformazione ω(F), la

67
quale rappresenta l’energia di deformazione accumulata nel materiale per
unità di volume, ed F rappresenta il gradiente di deformazione. Per tener
conto dell’Effetto Mullins, i ricercatori Ogden e Roxburgh proposero una
descrizione basata su una funzione energia potenziale di deformazione
modificata, del tipo ω(F,ɳ) dove ɳ è una variabile scalare detta “variabile
di danno”. Richiamando l’equazione rappresentativa del modello
hyperfoam di seguito,

2𝜇 1
𝜔 𝜆 𝜆 𝜆 3 𝐽 1 V.1
𝛼 𝛽

il potenziale di deformazione è costituito da una parte responsabile della

variazione di forma e una responsabile della variazione di volume. Ma nel
caso specifico del sughero, essendo questo un materiale altamente
comprimibile, il coefficiente di Poisson è tendente a zero. Quindi, per
l’equazione

𝜈
𝛽 V.2
1 2𝜈

i coefficienti 𝛽 risultano nulli; consegue quindi che la parte volumetrica

dell’energia di deformazione è nulla.
Per computare l’effetto Mullins si aggiunge all’energia di deformazione
una apposita funzione danno 𝜙 𝜂 , nel seguente modo

𝜔 𝜆 ,𝜂 𝜂𝜔 𝜆 𝜙 𝜂 V.4

dove 𝜔 𝜆𝑖 rappresenta la parte deviatorica dell’energia di deformazione.

La funzione danno 𝜙 𝜂 è una funzione continua di 𝜂 dove 𝜂 rappresenta
la variabile di danno; vale la condizione 0 𝜂 1. La condizione 𝜂 1
significa assenza di danno, cioe 𝜙 1 0, ovvero la deformazione sta
avvenendo sulla curva iperelastica primaria.

68
Poiché per definizione di iperelasticità la tensione è esprimibile come
gradiente della funzione energia potenziale di deformazione, derivando la
V.4 rispetto 𝜆 otteniamo

𝝈 𝜂, 𝜆 𝜂𝝈 𝜆 V.5

dove 𝝈 𝜆𝑖 è la tensione corrispondente alla curva iperelastica primaria e

relativa alla deformazione corrente 𝜆𝑖 . La tensione che si genera a causa
dell’effetto Mullins si ottiene moltiplicando 𝝈 𝜆𝑖 per lo scalare 𝜂 . Il
valore della variabile 𝜂 dipende da 3 parametri m, β ed r attraverso la
relazione

1 𝜔 𝜔
𝜂 1 𝑒𝑟𝑓 V.6
𝑟 𝑚 𝛽𝜔

dove erf 𝑥 è la funzione degli errori, 𝜔 è il massimo valore di 𝜔

in un punto durante la storia deformativa.

Figura 36

Nella figura 36 viene mostrata la scheda di compilazione dei parametri di

Mullins su Abaqus; i valori di m, β e r indicati sono quelli che danno il
miglior fitting delle curve velocità-tempo.

69
V.5 Implementazione del modello di danneggiamento
Il comportamento costitutivo di un materiale semplice, tipo acciaio, viene
determinato con prove meccaniche che servono a produrre stati di tensione
e deformazione ben definiti all’interno del materiale; tali prove sono per
esempio la prova di trazione, compressione, taglio ecc. Queste prove,
come visto precedentemente, vanno eseguite imponendo una certa velocità
di deformazione, più o meno elevata, e si può osservare quanto il materiale
sia sensibile alla velocita di deformazione. E’ possibile, inoltre,
quantificare la resistenza ultima del materiale e avere indicazione di
quando avverrà la rottura. Tutte queste informazioni si possono acquisire
fin quando sono misurabili; se, come nel caso in discussione, la velocità
di deformazione è elevatissima (si consideri che l’intervallo temporale di
acquisizione dei dati nella prova di penetrazione a impatto è di 0.03
secondi) diventa problematico procedere. Con la prova di penetrazione il
provino di sughero mostra segni di danneggiamento permanenti già a
velocita di 1,27 𝑚⁄𝑠; l’energia trasferita dal punzone al provino porta il
materiale in crisi di resistenza nell’arco di centesimi di secondo, fenomeno
che accade in maniera del tutto simile in molti dispositivi di sicurezza
antiurto. Il problema che ci siamo posti è stato innanzitutto quello di
assumere un criterio di rottura (o criterio di crisi) adatto al tipo di materiale
in prova e, in secondo luogo, implementare tale criterio su Abaqus.

 Criterio di rottura di Mohr

La scelta di un criterio di rottura dipende dalla modalità stessa di rottura,
se duttile o fragile; altro fattore è l’asimmetria di comportamento a
trazione e compressione. Il sughero mostra indubbiamente una asimmetria
di comportamento fra trazione e compressione, a compressione si presta a
deformazioni elevatissime senza rompersi mentre a trazione si trancia con
poco sforzo. Inoltre, osservando il provino di sughero perforato a 1,80
𝑚⁄𝑠 si nota una zona di rottura netta, pari circa all’aerea della sezione del

70
penetratore. Questo fa pensare che a certe velocità di penetrazione, dove
la rottura si verifica, il sughero tende a comportarsi come un materiale
fragile in cui le azioni di taglio producono decoesione del materiale. Visti
i presupposti sperimentali si è deciso di adottare il criterio di rottura di
Mohr. Nella figura 37 è mostrata la foto del provino impattato a 1,80 𝑚⁄𝑠.

Figura 37

Nella figura 38 è rappresentato graficamente il criterio di Mohr. Il cerchio

rosso a destra rappresenta il cerchio di Mohr di uno stato di trazione, il
cerchio rosso grande a sinistra rappresenta il cerchio di Mohr di uno stato
di compressione. Se quei cerchi sono massimi, ovvero rappresentativi
della rottura a trazione e compressione, le due rette tangenti a tali cerchi
(colore blu) rappresentano la “curva di danno” del materiale. Nel caso
mostrato in figura 38 si sono usati solo due cerchi di Mohr per ottenere le
curve di danno, in realtà occorrono diversi cerchi corrispondenti a diverse
prove poiché la curva di danno non è, a priori, approssimabile ad una retta.

71
 𝜏

𝑐 𝜎 𝜎 𝜎

𝛼
𝜎 𝜎 0 𝜎
𝜎 𝜎
𝜎

𝜎 𝜎
Figura 38

La retta di danno, ricavata mediante tangenza ai cerchi che rappresentano

la crisi del materiale, interseca l’asse delle 𝜏 e l’asse delle 𝜎. Nel punto di
intersezione con l’asse delle 𝜎 le tre tensioni principali sono tutte tre
uguali e rappresentano una trazione triassiale idrostatica; in questo caso il
materiale va subito in crisi senza l’apporto di alcuna tensione tangenziale.
E’ interessante osservare come la presenza di una compressione migliori
nei materiali fragili la resistenza al taglio rispetto a quella relativa al taglio
puro (cerchio di Mohr centrato nell’origine degl’assi). La presenza di una
compressione idrostatica, infatti, aumenta il raggio del cerchio ma allo
stesso tempo ne sposta il centro verso sinistra allontanando così la
condizione di crisi che, ripeto, avviene quando il cerchio dello stato
tensionale corrente è tangente (o addirittura fuoriesce) alla curva di danno.

 Implementazione del criterio di rottura di Mohr: subroutine

VUSDFLD
Come già anticipato in precedenza, è stato utilizzato Abaqus per le
simulazioni agli elementi finiti. Abaqus consente, grazie alle sue librerie,
la modellazione della rottura di diversi materiali ma, purtroppo, non è

72
implementato il criterio di rottura di Mohr per materiali fragili. In questi
casi è comunque possibile implementare un criterio di rottura andando a
richiamare una subroutine in linguaggio Fortran opportunamente
compilata. Tra le tante subroutine che Abaqus mette a disposizione per
risolvere problemi specifici, quella relativa alla modellazione di
danneggiamento di un materiale si chiama VUSDFLD, che è l’acronimo
di “User Defined Field” per risolutore esplicito. In figura 39 è
rappresentata la finestra di dialogo Abaqus per l’implementazione della
subroutine.

Figura 39

In figura 40 è riportata, per conoscenza, una subroutine compilata in

linguaggio Fortran.

73
 cerchio di Mohr

condizione di rottura

Figura 40

74
Operativamente la subroutine VUSDFLD esegue il calcolo dello stato di
deformazione, ovvero calcola il tensore di deformazione (o di tensione, a
seconda di come vogliamo impostare l’analisi) in ogni elemento del
modello. Una volta acquisito il tensore di deformazione, la subroutine
calcola il cerchio di Mohr nei parametri che lo rappresentano, cioè il raggio
e il centro del cerchio. A questo punto si impone la condizione di rottura
di Mohr scrivendo l’equazione della “retta di danno” e successivamente si
esegue la verifica di avvenuta rottura con un ciclo di IF. Come discusso
nel paragrafo precedente la rottura avviene se il raggio del cerchio supera
un valore limite (epsFail) e, in tal caso, la subroutine procede con
l’istruzione StateNew(K,1)=0. Operativamente quest’ultima istruzione
annulla (spegne) gli elementi che non hanno superato la verifica di rottura,
tali elementi sono quelli che nella sperimentazione vengono strappati via
dal punzone e non collaborano più alla resistenza.
I parametri della retta di danno, ovvero l’intercetta con l’asse verticale e
la pendenza, si sono determinati mediante una serie di simulazioni tese a
trovare il miglior fitting fra le curve sperimentali e quelle ottenute
numericamente,

V.6 Preparazione del modello FEM

E’ stato costruito un modello assialsimmetrico 2D del blocco cilindrico di
sughero e del penetratore come mostrato nella figura 41. Il blocco di
sughero è stato discretizzato in 1066 elementi e modellato con elementi
assialsimmetrici 2D ad integrazione ridotta (CAX4R). Il penetratore è
stato discretizzato come una superficie rigida, con una massa di 3.055kg e
una velocità iniziale pari ai valori sperimentali. Il modello di sughero è
stato vincolato sulla superficie laterale esterna mediante un vincolo tipo
incastro. Considerando che il provino di sughero presenta una rugosità
elevata, ci si deve aspettare un attrito elevato; per questo motivo si è

75
modellato il contatto tangenziale fra il provino e il penetratore inserendo
un coefficiente di attrito di 0,5. Si è considerato un valore della densità
pari a 140 𝐾𝑔⁄𝑚 .

Figura 41

Ritengo utile, anche ai fini di eventuali future consultazioni, elencare

sommariamente la sequenza logica di costruzione del modello:
1. Si costruiscono i modelli 2d del provino di sughero e del
penetratore;
2. Si procede con l’inserimento dei parametri costitutivi e fisici del
modello di sughero, quindi densità, parametri hyperfoam,
parametri viscoelastici (serie di Prony del modello di Maxwell
generalizato), parametri di Mullins, attivazione della subroutine
VUSDFLD;

76
 3. Si imposta lo step ovvero l’intervallo temporale di simulazione che,
come nel caso sperimentale, è pari a 0,03 secondi;
4. Nella sezione Field Output si vanno a settare tutti i valori che
vogliamo ottenere al termine della simulazione; in particolare a noi
interessano i valori del raggio del cerchio di Mohr (ovvero la
deformazione tangenziale massima) e la deformazione di rottura
epsfail;
5. Nella sezione History Output si va a richiedere l’andamento
temporale di grandezze di varia natura; per i nostri scopi sono state
calcolate le curve velocità-tempo di un punto rappresentativo del
penetratore per le tre diverse velocità di impatto;
6. Vengono impostate le condizioni cinematiche di vincolo del
modello di sughero, le condizioni cinematiche e dinamiche del
penetratore (massa e velocità) e le condizioni di attrito fra modello
di sughero e penetratore;
7. Infine viene creato il Job di simulazione impostando la chiamata
della subroutine VUSDFLD;

V.7 Analisi FEM e risultati

Una volta preparato il modello e compilata la subroutine VUSDFLD, si
sono eseguite una serie di simulazioni per ogni velocità di impatto. Si deve
a questo punto fare un’osservazione riguardante i dati che costituiscono
l’intero modello: come già descritto precedentemente, i dati relativi alla
parte hyperfoam, alla serie di Prony della parte viscoelastica e i parametri
di Mullins si sono ritenuti sufficientemente attendibili in quanto testati in
precedenti lavori di ricerca (vedi bibliografia AIAS2020-1496). Pertanto,
essendoci posti l’obiettivo di modellare l’intero fenomeno fino a rottura
ed avendo come unico riscontro le curve velocità-tempo sperimentali,

77
abbiamo impostato le simulazioni partendo da valori fittizi dei parametri
di rottura Mohr 𝑐; 𝛼 (figura 42).

𝒄 𝜶

Figura 42

Procedendo per approssimazioni successive, abbiamo determinato le

diverse coppie 𝑐; 𝛼 e valutato gli andamenti qualitativi delle curve
velocità-tempo del penetratore, a parità di velocità di impatto. Si è
proceduto, quindi, a scegliere la coppia 𝑐; 𝛼 che fornisse il miglior
fitting fra la curva sperimentale e la curva numerica. Questa operazione di
scelta dei parametri ha necessariamente richiesto un compromesso fra
fitting globale, intercettazione della velocità di uscita e il tempo di
rimbalzo. Le curve, infatti, seppur simili si presentavano traslate. Trovate
le coppie 𝑐; 𝛼 che meglio rappresentano il comportamento del materiale
in studio, le abbiamo riportate e messe in evidenza nelle tabelle seguenti,
insieme ai rispettivi grafici velocità-tempo.
Come ci aspettavamo c’è una leggera differenza fra i valori dei parametri
nei tre casi di velocità, ma tale differenza è stata al di sotto delle nostre
aspettative. Possiamo affermare che il modello impostato coglie in
maniera sufficientemente accurata il comportamento globale del
materiale, dalla fase di schiacciamento fino al rimbalzo del penetratore.
Dopo il rimbalzo si ha una leggera perdita della curva, specialmente nel
caso di velocità 1,27 m/s, con una previsione più ottimistica di recupero di
energia cinetica. Questo è probabilmente dovuto ad una non ottimale
calibrazione dei parametri di Mullins che, come descritto nei capitoli
precedenti, gioca un ruolo importante nella fase di scarico della
deformazione.

78
 Tabella simulazioni 1 𝒗𝒊𝒏 𝟎, 𝟕𝟗 𝒎/𝒔
𝑐 𝛼 danno 𝑣 FEM (m/s) 𝑣 EXP (m/s)
0,74 0,20 no 0,476 rimbal. 0,32
0,60 0,20 no 0,476 rimbal. 0,32
0,40 0,20 no 0,476 rimbal. 0,32
0,30 0,20 parziale 0,296 rimbal. 0,32
0,35 0,20 parziale 0,372 rimbal. 0,32
0,33 0,20 parziale 0,329 rimbal. 0,32
0,34 0,20 parziale 0,329 rimbal. 0,32
0,34 0,19 parziale 0,358 rimbal. 0,32
0,36 0,20 parziale 0,370 rimbal. 0,32
0,36 0,22 parziale 0,425 rimbal. 0,32
0,35 0,22 parziale 0,432 rimbal. 0,32
0,33 0,22 parziale 0,408 rimbal. 0,32
0,37 0,20 parziale 0,370 rimbal. 0,32

v (m/s)
1

0,5

0
0 5 10 15 20 25 30 35
t (ms)
‐0,5

EXP FEM
‐1

‐1,5

79
 Tabella simulazioni 2 𝒗𝒊𝒏 𝟏, 𝟐𝟕 𝒎/𝒔
𝑐 𝛼 danno 𝑣 FEM (m/s) 𝑣 EXP (m/s)
0,74 0,20 no 0,730 rimbal. 0,32
0,60 0,20 no 0,730 rimbal. 0,32
0,50 0,20 no 0,450 rimbal. 0,32
0,40 0,20 parziale 0,480 rimbal. 0,32
0,20 0,20 penetrato -0,553 0,32
0,45 0,20 parziale 0,435 rimbal. 0,32
0,50 0,22 no 0,435 rimbal. 0,32
0,48 0,18 parziale 0,435 rimbal. 0,32
0,35 0,18 parziale 0,076 rimbal. 0,32
0,40 0,19 parziale 0,432 rimbal. 0,32
0,38 0,20 parziale 0,235 rimbal. 0,32
0,30 0,20 penetrato -0,098 0,32

v (m/s)
0,3

‐0,2 0 5 10 15 20 25 30 35
t (ms)

‐0,7 EXP FEM

‐1,2

‐1,7

80
 Tabella simulazioni 3 𝒗𝒊𝒏 𝟏, 𝟖𝟎 𝒎⁄𝒔
𝑐 𝛼 danno 𝑣 FEM (m/s) 𝑣 EXP (m/s)
0,60 0,15 parziale 0,682 rimbal. -1,1
0,55 0,15 parziale 0,262 rimbal. -1,1
0,55 0,10 parziale 0,380 rimbal. -1,1
0,50 0,15 parziale 0,268 rimbal. -1,1
0,45 0,15 penetrato -0.431 -1,1
0,45 0,10 parziale 0,132 -1,1
0,45 0,18 penetrato -0,403 -1,1
0,40 0,18 penetrato -0,707 -1,1
0,35 0,18 penetrato -1.055 -1,1
0,30 0,18 penetrato -1.196 -1,1
0,38 0,18 penetrato -0,781 -1,1
0,35 0,195 penetrato -0,970 -1,1

v (m/s)
0,3

‐0,2 0 5 10 15 20 25
t (ms)

‐0,7

‐1,2

‐1,7 EXP FEM

‐2,2

81
V.8 Mappe di deformazione FEM
I risultati grafici delle simulazioni sono riportati nelle figure seguenti, in
termini di deformazione a taglio massima e deformazione a taglio
ammissibile. Nella tabella 3 si hanno le immagini del provino un istante
prima della rottura e, come si può notare, il materiale va in crisi allo stesso
valore di deformazione nei tre casi di velocità di impatto. Nella zona
centrale del provino, dove si risente maggiormente dell’urto, nonostante si
possa ipotizzare una componente idrostatica di deformazione che gioca a
favore della resistenza, il materiale è prossimo alla rottura. In tabella 4
sono riportate le immagini relative a diversi istanti dopo la rottura. Come
ci si poteva aspettare, essendo il materiale arrivato a rottura non è più in
grado di disporre di una componente di deformazione idrostatica che gli
favoriva la capacità di resistenza, ci ricordiamo infatti che il criterio di
Mohr utilizzato si basa proprio su questo concetto. Ossevando la tabella 4
notiamo che già a velocità di 0,79 𝑚⁄𝑠 interviene la rottura nei primi strati
di sughero a contatto col penetratore; lateralmente al penetratore il
materiale ha ancora un’ultima resistenza indicata dalle zone rosse
prossime a rottura. Alla velocità di 1,27 𝑚⁄𝑠 il fenomeno è ancor più
accentuato e la zona di rottura è maggiormente estesa, ma non si verifica
il trapasso del punzone in quanto gli strati sottostanti permangono a stati
deformativi lontani dalla crisi. Nel caso delle velocità di 0,790 𝑚⁄𝑠 e
1,270 𝑚⁄𝑠, non essendo tutto il materiale arrivato a rottura, questo rilascia
l’energia elestica accumulata facendo rimbalzare il punzone con una
velocità che, in fase sperimentale, è stato possibile determinare.
Nell’ultimo caso, la velocità del punzone di 1,80 𝑚⁄𝑠 produce la rottura
completa del modello di sughero; se osserviamo le due immagini possiamo
affermare che la rottura avviene in maniera fragile poiché sembra
danneggiata solo la zona coinvolta dal punzone mentre le zone limitrofi
restano a valori di deformazione molto bassa. In questo ultimo caso,

82
l’avvenuta perforazione completa del provino fa si che l’impattatore non
possa rimbalzare e prosegua la sua corsa verso il basso ad una velocità
inferiore rispetto al valore iniziale.

83
 1.0 J – 0,79 𝑚⁄𝑠
2.5 J ‐ 1,27 𝑚⁄𝑠
5.0 J ‐ 1,80 𝑚⁄𝑠 Deformazione a taglio massima Deformazione a taglio ammissibile

Tabella 3: Mappe deformative pre-rottura

84
 1.0 J – 0,79 𝑚⁄𝑠
2.5 J ‐ 1,27 𝑚⁄𝑠
5.0 J ‐ 1,80 𝑚⁄𝑠 Deformazione a taglio massima Deformazione a taglio ammissibile

Tabella 4: Mappe deformative post-rottura

85
 Capitolo VI

Conclusioni
In questo lavoro è stato riprodotto numericamente, sfruttando la
simulazione FEM, la prova di urto eseguita su dei provini di sughero
mediante apparecchiatura con penetratore a caduta Instron Ceast. Il
sughero è un materiale che ben si presta in applicazioni dove esistano
elevate deformazioni ed elevate velocità di deformazione; si pensi per
esempio ai dispositivi di sicurezza passiva come caschi da moto,
smorzatori di vibrazioni, ma anche come materiale d’anima per strutture
in composito. Appurato il comportamento iperelastico e viscoelastico in
regime di deformazione quasi-statico, per completare lo spettro costitutivo
del materiale restava da analizzare il comportamento a rottura per velocità
di deformazione tipiche di un crash. Partendo da tali presupposti, mediante
l’utilizzo di Abaqus CAE, abbiamo rimodellato il materiale aggiungendo
due ulteriori contributi: l’effetto Mullin, per tener conto del
danneggiamento a seguito del ciclo di carico-scarico, e la modalità di
rottura del materiale ipotizzando valido il criterio di Mohr. Grazie alle
curve sperimentali velocità-tempo rese disponibili per tre diverse velocità
di impatto, dopo una serie di simulazioni siamo risaliti a ritroso ai
parametri di rottura di Mohr cercando il miglior fitting fra la curve
velocità-tempo prodotte numericamente e quelle sperimentali. Dall’analisi
delle curve possiamo affermare che il modello proposto coglie, in maniera
soddisfacente, la fase di schiacciamento del materiale. La fase di rimbalzo
del penetratore, in tutti tre i casi di velocità di impatto, viene colta in
maniera non del tutto esatta con una sovrastima della velocità di rimbalzo.
Nonostante i ripetuti tentativi di simulazione con parametri di rottura

86
diversi, abbiamo concluso che il best fitting richiedeva una calibrazione
più accurata dei parametri di Mullins. L’analisi delle mappe di
deformazione prodotte numericamente danno altresì conferma che il
criterio di rottura di Mohr ipotizzato è corretto. Sulla parte più periferica
del provino infatti, non essendoci l’azione contenitiva di un vincolo che
promuove una certa compressione idrostatica, il materiale abbassa la
propria resistenza e va in crisi a causa di deformazioni da taglio. La rottura
che sopraggiunge in questi casi è una tipica rottura fragile e il modello
implementato ne ha colto, in maniera sufficientemente accurata, i tratti
tipici.

87
Riferimenti bibliografici:
 Characterization of hyperelastic rubber-like materials by biaxial
and uniaxial stretching tests based on optical methods (Sasso,
Amodio);
 Costitutive modelling of large strain time-dependent behaviour of
elastomers (Bergstrom, Boyce);
 Large elastic deformation of isotropic materials (Rivlin);
 Modellazione costitutiva e analisi numerico sperimentale del
comportamento a impatto di schiume di pvc (Sasso, Lattanzi);
 Modelling impact response of agglomerated cork (Sousa);
 Static and dynamic characterization of agglomerated cork and
related sandwich structures (Sasso, Sarasini);
 Testing Elastomers for Hyperelastic Material Models in finite
element Analysis (Axel);
 The viscoelastic property of cork (Mano);
 Non-linear elastic deformations (Ogden);
 Non-linear finite element for continua and structure (Belytschko);
 Polymer engineering science and viscoelasticity (Brinson)
 Scienza delle costruzioni Vol.1 (Baldacci);
 Manuale di Abaqus 2021

88