2007.DellaPorta

UNIVERSITÀ DEGLI STUDI DI NAPOLI “FEDERICO II”

FACOLTA’ DI INGEGNERIA

DIPARTIMENTO DI INGEGNERIA STRUTTURALE

TESI DI LAUREA

Modellazione Numerica d’Elementi Strutturali Sottoposti a Carichi da Esplosione

Relatori Candidato Ch. mo Prof. Ing. Gaetano Manfredi Giuseppe Della Porta Ch. mo Prof. Ing. Andrea Prota Matr. 437/99 Correlatore Ing. Domenico Asprone

ANNO ACCADEMICO 2005 – 2006

A quanti hanno saputo sostenermi con amore, spirito di sacrificio e sopportazione, a voi tutti va la mia gratitudine e

il mio affetto

“È la somma che fa il totale” Principe Antonio de Curtis

In arte Totò

IINNDDIICCEE

4

Indice

Introduzione ______________________________________ 7

Capitolo 1 _______________________________________ 10

Caratterizzazione dell’esplosione ____________________ 10

1.1 – Propagazione dell’onda d’urto (Shock Waves)__________________________ 12 1.1.1 – Esplosione al suolo_____________________________________________________14 1.1.2 – Esplosione in Aria _____________________________________________________16 1.1.3 –L’esplosione in centri abitati ______________________________________________23

1.2 – Inneschi__________________________________________________________ 26 1.3 – Dispositivi Esplosivi ________________________________________________ 28 1.4 – Materiali Esplosivi _________________________________________________ 30

Esplosivi Semplici ___________________________________________________________31 Esplosivi Composti __________________________________________________________33

Capitolo 2 _______________________________________ 36

Descrizione di un test d’esplosione su una piastra in c.a. _ 36

2.1 – Descrizione della Prova Sperimentale _________________________________ 37 2.2 – Evidenze Sperimentali______________________________________________ 41 2.3 – Finalità della Prova Sperimentale ____________________________________ 43

Capitolo 3 _______________________________________ 44

Modellazione delle azioni prodotte dall’esplosione. _____ 44

3.1 – Curva di pressione adottata _________________________________________ 46 Coefficiente d’equivalenza_____________________________________________________47 Distanza Ridotta Z ___________________________________________________________48 Picco di Pressione Δp(Z) ______________________________________________________50 Coefficiente Sperimentale α. ___________________________________________________50 Durata della Fase Positiva t+ ___________________________________________________51 Curva di Pressione P(Z, t) _____________________________________________________52

Capitolo 4 _______________________________________ 60

Metodi di Discretizzazione__________________________ 60

4.1 – Finite Element Method (FEM) _______________________________________ 60 4.2 – Classificazione dei Metodi Particellari ________________________________ 65 4.3 – Applicazione dei Metodi Particellari __________________________________ 68 4.4 – Smoothed Particle Hydrodynamics (SPH) _____________________________ 69

4.4.1 – Metodi d’interpolazione SPH_____________________________________________72 4.4.2 – Teoria d’Approssimazione SPH___________________________________________78 4.4.3 – Approssimazione di Derivate SPH_________________________________________80

Formula gradiente I ________________________________________________________80 Formula gradiente II _______________________________________________________81

IINNDDIICCEE

5

Formula gradiente III_______________________________________________________82 Equazione di continuità in SPH_______________________________________________85 Equazione del Moto in SPH _________________________________________________86 Equazione dell’energia in SPH _______________________________________________89

4.4.4 – I limiti del Metodo SPH _________________________________________________91 Instabilità di tensione_______________________________________________________91 Condizioni di vincolo ______________________________________________________95

Capitolo 5 _______________________________________ 97

Modellazione Numerica con LS-Dyna Version 970. _____ 97 5.1 - Software LS-Dyna Version 970 _______________________________________ 97 5.2 – Modellazione della Prova Sperimentale_______________________________ 101

5.2.1 – Modello FEM ________________________________________________________102 5.2.2 – Modello SPH ________________________________________________________111

Capitolo 6 ______________________________________ 117

Legami Costitutivi del Calcestruzzo _________________ 117

6.1 – Pseudo – Tensor(Mat_016) _________________________________________ 117 6.2 – Concrete Damage (Mat_072) _______________________________________ 125

6.2.1 – Parametri che definiscono le superfici di rottura _____________________________128 6.4 – Dynamic Increase Factor (DIF) _____________________________________ 137

• Cowper and Symonds. __________________________________________________139 • Legame Tensione-Deformazione Parametrico. _______________________________139 • Curva di Strain-Rate. ___________________________________________________140 6.4.1 – Definizione della Curva di Strain-Rate ____________________________________140 6.4.2 – Origini della formulazione CEB-FIP ______________________________________141

6.3 – Equation of State Form 8 (Tabulated Compaction) _____________________ 149

Capitolo 7 ______________________________________ 161

Risultati delle Simulazioni Numeriche _______________ 161

7.1 – Prova Sperimentale I e II __________________________________________ 162 • Modello FEM, suddivisione della piastra in 11520 Element_Solid, Prova 1: ________163 • Modello FEM, suddivisione della piastra in 11520 Element_Solid, Prova 2: ________164 • Modello SPH, suddivisione della piastra in 11520 Element_SPH, Prova 1: _________166 • Modello SPH, suddivisione della piastra in 11520 Element_SPH, Prova 2: _________168

7.2 – Prova Sperimentale III ____________________________________________ 169 • Modello FEM, suddivisione della piastra in 11520 Element_Solid, Prova 3: ________171 • Modello SPH, suddivisione della piastra in 11520 Element_SPH, Prova 3: _________176

7.3 – Prova Sperimentale IV ____________________________________________ 181

CONCLUSIONI _________________________________ 186

Indice delle Figure _______________________________ 191

Indice dei Grafici ________________________________ 192

Indice delle Tabelle_______________________________ 193

IINNDDIICCEE

6

Bibliografia _____________________________________ 195

Ringraziamenti __________________________________ 198

IINNTTRROODDUUZZIIOONNEE

7

Introduzione

Nel campo dell’Ingegneria Civile si assiste ad un sempre crescente

interesse nei confronti delle azioni dinamiche prodotte da esplosioni

applicate alle strutture.

Infatti, i recenti attacchi terroristici portati al cuore del mondo

occidentale negli ultimi anni, hanno fatto sì che si alzassero i livelli di

guardia dei paesi maggiormente colpiti ed interessati.

La crescente preoccupazione dei paesi occidentali, ha quindi focalizzato

l’attenzione sull’argomento della protezione di quelle strutture che, per

le proprie sensibilità strategiche e sociali, rientrano tra i probabili

obiettivi terroristici. Tra queste si annoverano le ambasciate, le sedi di

governo, le basi militari e logistiche, ma anche scuole, ospedali, stadi,

centri commerciali, luoghi di grande affollamento tra cui anche

metropolitane, aeroporti e mezzi di trasporto collettivi.

Lo stesso ambito scientifico si è adoperato per caratterizzare il

comportamento di tali strutture sottoposte a carichi dinamici d’elevata

intensità e breve durata, quali sono quelli prodotti dalle esplosioni.

Il problema si presenta al quanto complesso, e si può articolare in tre

diversi aspetti fondamentali, che sono così elencati:

• Valutazione dei carichi dinamici prodotti da un’esplosione ed

applicati su un elemento strutturale;

• Caratterizzazione del comportamento meccanico dei materiali per

uso strutturale, sottoposti a carichi dinamici d’elevata intensità;

• Individuazione ed implementazione, mediante l’uso del

calcolatore, di un metodo d’analisi numerico per la descrizione e

soluzione del problema fisico-meccanico così definito.


8

Per quanto concerne la valutazione delle azioni prodotte da

un’esplosione, queste possono essere definite da curve di pressione

d’elevate intensità ma di durata ridottissima. Di quest’onda di pressione

è possibile considerare una fase positiva di sovrapressione e una negativa

di depressione.

Sperimentalmente si è potuto osservare come il calcestruzzo e l’acciaio,

sottoposti a carichi dinamici d’elevata intensità, presentino dei valori di

resistenza meccanica maggiori rispetto a quelle calcolate in condizioni di

carico quasi statiche. Ad oggi ancora non appare chiaro il motivo fisico

in virtù del quale si manifesta quest’incremento di resistenza.

Ciononostante, per la modellazione di un elemento strutturale realizzato

in c.a. non si può prescindere dal considerare tal effetto, che va sotto il

nome di Strain-Rate-Sensitivity.

Le metodologie d’analisi maggiormente impiegate, per la descrizione del

fenomeno fisico, prevedono una complessa discretizzazione del

continuo. A tal fine, è possibile distinguere diversi gruppi di metodi, tra

cui:

• Lagrangiani, utilizzati soprattutto per la discretizzazione

d’elementi solidi, e per la risoluzione di problemi meccanici, in

cui si osservano deformazioni non eccessivamente elevate;

• Euleriani, utilizzati soprattutto in fluidodinamica, e a differenza

dei precedenti offrono la possibilità di trattare problemi che

prevedono grandi deformazioni;

• ALE (Arbitrarian Lagrangian Eulerian), consentono di trattare

deformazioni di qualsiasi entità, la loro precisione è funzione di

quanto sia corretta la previsione delle deformazioni attese;


9

• SPH (Smoothed Particle Hydrodynamics), appartengono ai metodi

di tipo Meshfree, discretizza il continuo mediante un numero finito

d’elementi particellari.

L’oggetto del lavoro di tesi risiede nell’affrontare ognuna delle

problematiche in precedenza descritte. Ragion per cui, si è scelto di

modellare numericamente una prova sperimentale condotta su una

piastra in c.a. sottoposta a carichi da esplosione. Nei capitoli successivi

si descriveranno le metodologie adottate per il calcolo delle azioni

esercitate sull’elemento strutturale prodotte dall’esplosione, e quelle

utilizzate per la discretizzazione e l’analisi numerica della prova

simulata.

Inoltre, si farà un esplicito riferimento al codice di calcolo impiegato per

la suddetta modellazione, e sulle scelte riguardanti il legame costitutivo

implementato al fine di modellare il calcestruzzo armato, in modo tale da

poter considerare il già menzionato effetto di Strain-Rate-Sensitivity.

A conclusione della tesi si riporteranno i risultati ottenuti mediante

l’utilizzo di due diversi metodi d’analisi, tra cui il tradizionale metodo

degli elementi finiti (FEM), che rientra tra i metodi Lagrangiani prima

citati, e il più innovativo Smoothed Particle Hydrodynamics (SPH).

Dal confronto operato tra i risultati conseguiti mediante i due metodi

d’analisi, e tra gli stessi e quelli ottenuti in seguito alle prove

sperimentali se ne trarranno le successive conclusioni.

CCAAPPIITTOOLLOO 11 –– CCAARRAATTTTEERRIIZZZZAAZZIIOONNEE DDEELLLL’’EESSPPLLOOSSIIOONNEE

10

Capitolo 1 Caratterizzazione dell’esplosione

L’esplosione è definita come un improvviso e violento rilascio d’energia

meccanica, chimica o nucleare, normalmente con produzione di gas ad

altissima temperatura e pressione[1]. Questo rilascio d’energia si

attribuisce ad una trasformazione dello stato di una massa gassosa, solida

o liquida, impiegata per la costruzione della carica esplosiva. Per Energia

Meccanica s’intende la somma dell’energia Cinetica e di quella

Potenziale, mentre l’Energia Chimica è causata della formazione o

rottura di legami chimici di qualsiasi tipo. Con Energia Nucleare

s’intendono tutti quei fenomeni in cui si ha la produzione di energia in

seguito a trasformazioni nei nuclei atomici, quali possono essere fissioni

o fusioni nucleari. I volumi di gas prodotti dall’esplosione sono circa

10000-30000 volte quelli iniziali, la temperatura che ne consegue

raggiunge i 3000 gradi centigradi, mentre la pressione si aggira intorno

alle 20000 atm con picchi di 150000 atm[2]. Il gas d’esplosione, ad

elevata temperatura e pressione, si propaga nell’ambiente circostante

sotto forma di un onda d’urto, che in assenza d’ostacoli si espande come

una superficie sferica centrata nel centro dell'esplosione. Incontrando

degli ostacoli esercita su di loro una forza tanto maggiore quanto la

superficie investita, e quanto più vicina al centro dell'esplosione [1].

L’effetto distruttivo di un’esplosione dipende da vari fattori tra cui la

velocità dell’esplosione, il calore e i gas da lei prodotti, e chiaramente

dalle pressioni realizzabili.

Di conseguenza gli effetti prodotti da esplosivi d’uso militare e civile

sono chiaramente differenti, avendo caratteristiche diverse in particolar

modo gli esplosivi militari tra cui quelli plastici, come il C-4, hanno un


11

enorme potere distruttivo legato soprattutto alla loro elevata velocità di

detonazione, all’effetto prodotto dall’onda di pressione sviluppata

dall’esplosione, e in misura minore agli effetti prodotti dalle schegge. Gli

esplosivi commerciali d’uso civile, impiegati per la demolizione di

strutture esistente o per l’apertura di scavi in roccia, sono solitamente

utilizzati con cariche intasate, in altre parole introdotte in fori realizzati

nella roccia o nel terreno.

Le esplosioni chimiche sono suddivise in deflagrazioni, nelle quali la

propagazione della reazione chimica d’esplosione è una forma di

combustione che procede nel materiale a velocità subsonica, e

detonazioni, nelle quali la reazione chimica d’esplosione non è una

combustione ma una decomposizione diretta della molecola d’esplosivo,

innescata direttamente dall'onda d'urto, la reazione d’esplosione procede

quindi alla velocità del suono in quella particolare sostanza attraverso

tutto il materiale, e la pressione e temperatura finale dei prodotti di

reazione sono quindi molto più elevati [1]. Per caratterizzare la potenza

di un esplosivo è possibile condurre una prova sperimentale, che prevede

l’esplosione di una carica all’interno di un blocco di piombo le cui

dimensioni sono standardizzate. Eseguendo in seguito una misurazione

del volume della cavità creatasi al suo interno, è possibile stabilire una

correlazione tra la potenza della gelatina esplosiva, ritenuto uno degli

esplosivi più potenti, con altre tipologie di cariche esplosive, associando

convenzionalmente alla gelatina una potenza pari a 100 [1]. In ambito

militare ma anche in campo civile si preferisce definire un carico

esplosivo equivalente, calcolando dei coefficienti d’equivalenza rispetto

al Trinitrotulene (TNT), per la quale si assume un coefficiente unitario

[1].


12

1.1 – Propagazione dell’onda d’urto (Shock Waves) I gas ad alta temperatura e pressione, prodotti dall’esplosione, sono

confinati dal mezzo circostante, e si propagano al suo interno sotto forma

di un’onda di pressione. L’onda d’urto produce effetti devastanti sia a

breve sia a lunga distanza.

Gli effetti a lunga distanza sono funzione del mezzo attraverso il quale

avviene la propagazione.

Infatti, la massima pressione del gas d’esplosione, dipende dalle

caratteristiche fisico-meccaniche del mezzo che circonda la carica

esplosiva, che sia aria, acqua o suolo.

In aria si avrà uno spostamento d’aria, in acqua uno scoppio subacqueo,

se invece l’esplosione avviene nel terreno l’onda d’urto si propaga allo

stesso modo delle onde sismiche, arrecando quindi danni alle strutture

adiacenti cosi come alle persone a contatto con la superficie investita

dall’onda.

L’eventuale presenza di una superficie di separazione, tra il gas prodotto

dall’esplosione, e il mezzo attraverso il quale si propaga, produce al

tempo stesso una riflessione e una trasmissione dell’onda di pressione al

mezzo circostante (area o acqua).

L’onda di riflessione, si muove in verso opposto all’onda di detonazione.

Questa riflessione può essere di compressione o di rarefazione, a seconda

che la densità del mezzo di propagazione è maggiore o minore di quella

del gas d’esplosione [2].

Ad esempio, se la carica è applicata su una piastra infinitamente rigida di

calcestruzzo armato, siccome la densità della piastra è molto maggiore

rispetto a quella del gas prodotto dall’esplosione, l’onda di riflessione

sarà di pressione, ed interamente rivolta verso l’alto.


13

Allo stesso modo, l’intensità della pressione dell’onda trasmessa al

mezzo circostante, è funzione della densità del mezzo. In ogni caso,

all’aumentare della distanza dal centro d’esplosione, l’intensità di

pressione tende a diminuire.

L’onda esplosiva genera inizialmente una sovrapressione e in seguito in

tempi più lunghi una depressione nota come risucchio, dovuta all’aria

che torna violentemente verso il centro dell’esplosione. Gli effetti

prodotti dall’azione combinata di queste due onde sono molto distruttivi,

in quanto in alcuni casi la sola onda di pressione che impatta su un

oggetto, quale può essere ad esempio un pannello in muratura, provoca

su di lui inizialmente delle semplici lesioni, senza alterare la propria

staticità, sarà poi l’onda di risucchio a provocarne il successivo crollo.

Tutto ciò avviene in un periodo di tempo dell’ordine di grandezza di

pochi millisecondi, tale da non far apprezzare l’effetto prodotto dalle due

onde singolarmente.

Gli effetti prodotti a breve distanza sono invece da attribuire ad onde

d’urto pulsanti, le quali attraversano e si riflettono sulle superfici libere

degli oggetti incontrati, provocando su di loro un incremento delle

tensioni che ne determina la rottura. A questi si sommano gli effetti

prodotti dalle fiammate, dai corpi incandescenti scagliati ad alte velocità,

e dalla propagazione del calore che può essere causa d’incendi e

d’ustioni gravi.

In alcuni casi i danni prodotti a breve distanza possono essere aggravati

dalla presenza all’interno dell’ordigno di chiodi o altri corpi contundenti

(Bomba Sporca), quali schegge di vetro o acciaio di varie dimensioni,

che l’esplosione scaglia ad alte velocità, tra i 1000 e 1500 m/s, possono

provocare lesioni a persone e cose d’entità rilevanti [1].


14

1.1.1 – Esplosione al suolo Ad esplosione avvenuta, dalle dimensioni radiali del cratere formatosi al

suolo, e note le caratteristiche geotecniche dello stesso, sarà possibile

risalire all’entità del carico esplosivo W. Noti, infatti, la tipologia del terreno, e con lui le sue caratteristiche

fisiche e meccaniche, attraverso delle relazioni empiriche, è possibile

calcolare il raggio del cratere R, attraverso la seguente espressione:

Dove con R s’indica il raggio del cratere, k è un coefficiente numerico

adimensionale, che assume valori diversi in funzione del tipo di terreno,

W è il peso della carica esplosiva, e infine n è l’indice di potenza

dell’esplosivo stesso [3].

Questa formulazione empirica, ricavata quindi dalle evidenze

sperimentali, è in accordo con la seguente teoria.

Si prenda in considerazione una carica esplosiva, applicata per semplice

contatto sulla superficie del suolo. Indichiamo rispettivamente con W ed

E il suo peso e la sua Energia Specifica di esplosione, ovvero l’energia

rilasciata dall’esplosione stessa.

L’energia specifica necessaria alla rottura di un volume V di terreno sarà:

Er = Vσrεr

Dove, σr e εr sono rispettivamente la tensione e la deformazione di

rottura del materiale attraversato dall’esplosione. L’energia specifica di

rottura del materiale sarà:

Es = σrεr

nkWR =


15

Dall’uguaglianza di queste due espressioni, ricaviamo il volume di

terreno asportato dall’esplosione, definito come segue:

Sapendo che l’energia necessaria per provocare la rottura di un volume V

di terreno, è uguale all’energia trasferita dalla carica al suolo (K=WβE).

Nell’ipotesi in cui la densità di massa del suolo è inferiore rispetto a

quella del gas prodotto dall’esplosione, approssimativamente si può

considerare un cratere di forma semisferica. Per cui, sapendo che il

volume della sfera si calcola come:

Sostituendo nella relazione precedente, ricaviamo l’espressione del

raggio:

Avendo posto:

espresso quindi in funzione delle caratteristiche meccaniche del suolo e

dell’energia d’esplosione della carica[3].

ss

r

EEW

EEV β

==

32 3RV π

=

3332

3 WkWEER

s

== βπ

βπ

32

3

sEEk =


16

1.1.2 – Esplosione in Aria Si vogliono descrivere gli effetti prodotti da un’esplosione che avviene

in aria libera, considerata puntuale ed istantanea, alla quale corrisponde

il rilascio di una certa quantità d’energia E.

Ciò che interessa definire, è l’incremento di pressione del mezzo

circostante la carica esplosiva, in un dato istante di tempo t, ad una certa

distanza radiale R(t) dal centro d’esplosione, supponendo che sia di tipo

sferico.

Rispetto alle condizioni imperturbate del mezzo, alle quali corrisponde

una certa pressione p1, in seguito al passaggio del fronte d’onda, si avrà

una nuova pressione p2, maggiore di quella iniziale.

L’ipotesi di base è che le trasformazioni di natura termodinamiche

associate all’esplosione, siano adiabatiche, in altre parole comporteranno

uno scambio di calore tra l’ambiente interno, in cui si ha l’esplosione, e

quello esterno, in cui si ha la propagazione della stessa. Inoltre, l’aria

attraverso la quale si propaga l’onda d’urto, è considerata come un gas

politropico, per il qual è vera la seguente espressione:

pvn = cost.

Dove p è la pressione, v è il volume specifico per unità di massa, e n è un

coefficiente numerico adimensionale, che assume valori diversi a

seconda che la trasformazione del gas sia isoterma o adiabatica, in

questo caso essendo adiabatica vale 4,21 [3].

Per definire la velocità, la pressione, e la densità di massa del mezzo

attraversato dal fronte d’onda ad un generico istante di tempo t, sarà

quindi necessario descrivere il moto del gas.


17

Tale moto è univocamente determinato, una volta noti la densità di

massa indisturbata dell’aria ρ1, l’energia sprigionata dall’esplosione E,

ed una volta fissati l’istante di tempo t e la distanza radiale r dal centro

d’esplosione.

Dalla successiva relazione si ricava la posizione radiale del fronte d’onda

ad un generico istante di tempo t:

Dove k è un parametro adimensionale, che assume valore pari a 1, 2, o 3,

a seconda che l’espansione del gas è monodimensionale, bidimensionale

o tridimensionale [3].

Da cui, la posizione radiale dell’onda d’urto all’istante t sarà:

Dove β è una costante adimensionale, definita in funzione del rapporto

politropico γ tra il calore specifico del gas a pressione costante, e a

volume costante dell’aria[3].

Attraverso la precedente relazione, si è quindi in grado di individuare ad

ogni istante, in seguito all’esplosione, la posizione del fronte d’onda.

Noto R(t), si può calcolare la velocità con la quale si propaga l’onda

d’urto nel mezzo circostante, come:

( ) kEtr += 21

21ρλ

( )kEttR

+

⎟⎟⎠

⎞⎜⎜⎝

⎛=

21

1

2

ρβ

dtdRtU =)(


18

Nel caso in questione si avrà:

Bisogna quindi risolvere un sistema di tre equazioni nelle incognite v2,

ρ2, p2, che indica rispettivamente la velocità, la densità di massa, e la

pressione dell’aria in seguito al passaggio del fronte d’onda, ad una

distanza radiale R(t), e all’istante di tempo t [3].

Le equazioni da scrivere per la risoluzione del problema, derivano dal

principio di conservazione della massa, della quantità di moto, e

dell’energia.

Le soluzioni del sistema d’equazioni saranno:

La sola densità di massa si esprime indipendentemente dalla velocità di

propagazione dell’onda d’urto U [3].

Nota l’espressione della velocità di propagazione del fronte d’onda U,

sostituendola in quella della pressione p2, si ottiene la pressione sul

fronte d’onda:

Questa sarà quindi espressa in funzione dell’energia specifica

d’esplosione E, del tipo di propagazione espressa mediante k, del

rapporto tra il calore specifico a pressione e a volume costante γ, e della

posizione radiale R [3].

5,0

1

24

22)( ⎟⎟

⎠

⎞⎜⎜⎝

⎛+

=+

k

k

RE

ktU

ρβ

12

2 +=

γUv1

112 −

+=

γγρρ

1

2

2 12 ργ +

=Up

( ) ( ) kk

RkEp 1

128

24

2 ++= +

γβ


19

La pressione del fronte d’onda è quindi una funzione inversamente

proporzionale della distanza R dal centro d’esplosione.

In particolar modo al variare della costante k, che come già detto assume

valori diversi a seconda che la propagazione del fronte d’urto sia

monodimensionale, bidimensionale, o tridimensionale, questa

proporzionalità inversa sarà lineare, quadratica o cubica.

Nel caso di propagazione tridimensionale, l’attenuazione della pressione

sul fronte d’onda, sarà ancora più rapida all’aumentare della distanza

radiale dal centro d’esplosione.

Il picco di pressione, in altre parole la massima pressione conseguita in

un punto generico dello spazio, a distanza R dal centro d’esplosione, si

definisce come somma di 4 contributi:

Dove le 4 costanti, A1, A2, A3, e A4, sono determinate attraverso

un’interpolazione lineare dei risultati numerici di prove sperimentali, che

dipendono dal peso della carica esplosiva W, mediante la cosiddetta

Distanza Ridotta, a sua volta definita come:

Il concetto di distanza ridotta, s’introduce per individuare la distanza dal

centro d’esplosione, alla quale due cariche esplosive, con la stessa

energia specifica d’esplosione, ma con peso differente, W1 e W2,

generano lo stesso effetto in termini di pressione. Analogamente, per due

cariche esplosive con lo stesso peso W, ma con energie specifiche

433

221 A

RA

RA

RAp +++=Δ

3 WRZ =


20

d’esplosione diverse, E1 ed E2, s’individua la distanza alla quale le due

cariche generano la stessa pressione [3].

Nelle applicazioni pratiche, per la valutazione della pressione in un

generico punto dello spazio, individuato dalla distanza ridotta Z dal

centro d’esplosione, e in un generico istante di tempo t, sono impiegate

formulazioni semiempiriche come la seguente:

Dove p+y è la massima pressione della fase positiva, p0 è la pressione in

condizioni indisturbate del mezzo di propagazione, pressione

atmosferica nel caso dell’aria, t+ è la durata della fase positiva di

sovrapressione, α è un coefficiente di natura sperimentale definito in

funzione della distanza ridotta Z, e Δp(Z) è il picco di pressione alla

distanza ridotta Z [3].

La durata della fase positiva di sovrapressione, la si può valutare come

segue (M. A. Sadowsky):

Dove, il tempo è espresso in secondi, B è una costante che si assume

nella pratica pari a 1,3 [3].

+−

++ ⎟

⎠⎞

⎜⎝⎛ −Δ=−= t

t

y ettZppptZp

α1)(),( 0

RWBt 6310 −+ ⋅=


21

px

0p

+y

p

Fase negativa

Fase positiva

t(s)

2p(Kg/cm )

+tta

Figura 1. 1 – Andamento qualitativo di una Curva di Pressione prodotta da un’esplosione.

Per il picco di pressione Δp(Z), è possibile utilizzare diverse

formulazioni d’origini sperimentali, proposte da vari autori, e riferite al

solo TNT.

In ogni caso, tali formulazioni sono estendibili a qualsiasi tipologia

d’esplosivo, mediante l’introduzione di un coefficiente d’equivalenza al

TNT. Questo coefficiente, è calcolato come il rapporto tra il calore

specifico dell’esplosivo in questione, Q, e quello del TNT, QTNT [3]:

Per ognuna di queste espressioni, è possibile definire un campo di

validità in funzione della distanza ridotta Z. Nei rispettivi campi

TNTQQq =


22

d’applicazione, tra le varie formulazioni sussistono delle differenze.

Nella pratica, si preferisce riferirsi ad un valore medio di quello calcolato

con le singole formulazioni.

Di seguito si riportano le espressioni per il calcolo del picco di pressione,

e i rispettivi campi di validità [3]:

• Brode H.L.

• Petrowsky G.I.

• Sadowsky M.A.

• Henrich J.

232

23

/101,0019,085,5455,1975,0)(

/1017,6)(

cmKgpZZZ

Zp

cmKgpZ

Zp

≤Δ≤−++=Δ

≥Δ+=Δ

117,10)( 3 ≤−=Δ ZZ

Zp

1515,655,276,0)( 32 ≤≤++=Δ ZZZZ

Zp

1288,305,4662,0)(

13,01324,23262,01938,6)(

3,005,000625,03572,05397,5072,14)(

32

32

432

≥++=Δ

≤≤+−=Δ

≤≤+−+=Δ

ZZZZ

Zp

ZZZZ

Zp

ZZZZZ

Zp


23

1.1.3 –L’esplosione in centri abitati Nel momento in cui l’esplosione avviene in un centro abitato, la sua

propagazione in aria può essere descritta come se avvenisse attraverso un

insieme di canali [3]. In questo caso l’onda di sovrapressione incontrerà

strade, piazze ed edifici, che ne determineranno delle riflessioni. Le

strade sono quindi considerate come dei canali lateralmente confinati

dalla presenza degli edifici.

Si supponga di prendere in considerazione un’esplosione di tipo

puntuale, che avviene ad una distanza R0 dalla sezione d’imbocco del

nostro canale. L’onda d’urto prodotta dall’esplosione, inizialmente

sferica, raggiunta la sezione iniziale del canale, si trasforma da sferica in

piana. Questa trasformazione si deve alle riflessioni multiple dell’onda

contro le pareti laterali del canale, e alle sovrapposizioni di queste

ultime. In questo modo gli effetti dell’esplosione si risentono anche a

distanze maggiori rispetto al caso d’assenza d’ostacoli.

E’ interessante determinare la pressione prodotta nell’atmosfera ad una

certa distanza r dal centro d’esplosione. In generale, si può affermare che

questa sia pari a:

Dove r è la distanza dal centro d’esplosione, ed A una costante [3].

Se s’indica con r = R1 la distanza dal centro d’esplosione della sezione

iniziale del canale, dove la pressione vale p0, la costante A si calcola

come:

A = R1p0

rArp =)(


24

La pressione in una generica sezione x del canale sarà calcolata come:

Dove, p0 è la pressione calcolata sul fronte d’onda sferico alla distanza

R0 dalla carica esplosiva [3].

Dalla precedente espressione si evince come all’aumentare della distanza

dal centro d’esplosione si abbia una riduzione della pressione prodotta,

essendo questa ultima inversamente proporzionale alla distanza x.

In questa formulazione non si tiene conto del graduale passaggio della

superficie d’onda da sferica a piana, dell’influenza delle onde di

riflessione sulle pareti che limitano lateralmente il canale, e degli effetti

dissipatevi dovuti alla rugosità delle pareti del condotto.

E’ quindi necessario introdurre dei coefficienti correttivi che terranno

conto di questi fenomeni. S’introduce quindi un coefficiente α,

compreso tra 0 e 1, che moltiplicato per R0 riduce la distanza del centro

d’esplosione dalla sezione d’imbocco [3]. Il coefficiente amplificativo

Ω, che tiene conto della sovrapposizione delle onde di riflessione, tale

coefficiente si assume sempre inferiore all’unità. Infine un termine

esponenziale, in funzione del quale all’aumentare della distanza x lungo

il canale, si ha una riduzione della pressione in virtù degli effetti

dissipativi. Sperimentalmente, è possibile dimostrare che tali effetti

dissipativi, siano funzione della rugosità delle pareti laterali del condotto,

espresse mediante un coefficiente di attrito ξ, a sua volta funzione della

dimensione trasversale del canale D, e dell’altezza relativa delle asperità

h [3].

xRRpxp+

=1

10)(


25

L’espressione corretta della pressione prodotta ad una generica ascissa x

del canale sarà [3]:

Da questa formulazione, si evince che per ridurre la pressione è

sufficiente aumentare la dimensione trasversale D del canale stesso.

Questo si può fare realizzando delle camere d’espansione, che sono degli

allargamenti della sezione trasversale corrente. A valle di questa camera

d’espansione, nel momento in cui si ha un successivo restringimento del

condotto, la pressione è calcolata come [3]:

Dove S1 e S2, sono rispettivamente la dimensione trasversale del

condotto prima e dopo l’allargamento, cosi come p1 e p2, sono le relative

pressioni, all’ingresso del canale d’espansione, e nel canale d’uscita,

quest’ultima p2 sarà incognita.

Quando lungo il canale principale, l’onda di sovrapressione incontra

delle curve o delle diramazioni, la pressione che si avrà a valle del

cambio di geometria del condotto, sarà valutata come un’aliquota della

pressione iniziale p1:

p2 = kp1

Dove k è un coefficiente di proporzionalità, che varia tra 0 e 1, e assume

valori diversi in funzione del tipo di diramazione o cambio di geometria

del canale [3].

( )Dx

exR

Rpxpξ

αα 4,0

1

10)(−

+Ω=

8,0

1

2

2

1⎟⎟⎠

⎞⎜⎜⎝

⎛=

SS

pp


26

1.2 – Inneschi Gli Inneschi possono essere realizzati con spolette, corde di detonazione,

micce d’accensione o detonatori, in grado di provocare l’esplosione nei

tempi e nei modi desiderati; molti di questi elementi sono

sufficientemente piccoli da poter essere facilmente nascosti in un veicolo

e quindi passare inosservati ad un controllo [2].

In alcuni casi l’innesco di un esplosivo può indurre ad una reazione a

catena, dovuto all’innesco successivo d’altri esplosivi nelle vicinanze.

Ciò non sempre avviene, in quanto questo fenomeno e vincolato da una

serie di fattori tra cui la distanza e la sensibilità degli esplosivi, la

violenza dell’esplosione e il mezzo attraverso il quasi si propaga.

La scelta del sistema d’innesco dipende non solo dal tipo d’esplosione e

di danno che si vuole provocare ma anche dal tipo d’esplosivo.

Fatta eccezione per la polvere nera che può facilmente essere innescata

mediante l’accensione di una miccia, tutti gli altri esplosivi hanno

bisogno di un detonatore[2].

La detonazione è un fenomeno chimico-fisico costituito da una

esplosione che si propaga ad una velocità supersonica costante, che

genera un'onda d'urto il cui campo di velocità a valle può essere ancora

supersonico (detonazione forte) oppure subsonico (detonazione debole)

[2]. La velocità di detonazione è una funzione della densità del solido

impiegato per la costruzione della carica esplosiva. Superato un certo

limite di densità, la velocità di detonazione tende a diminuire a causa

delle difficoltà con le quali si sviluppano le reazioni chimiche, che sono

artefici della detonazione stessa. I Detonatori utilizzati per innescare gli

esplosivi normali sono solitamente dei tubi sottili d’alluminio o di rame,

contenenti vari tipi d’esplosivi primari e secondari presenti in piccole


27

quantità, che li rendono molto sensibili alle azioni esterne tra cui

percussioni, shock e calore, per tale motivo vanno maneggiati con

cautela e mai trasportati insieme all’esplosivo. Gli esplosivi primari sono

molto sensibili agli urti, agli sfregamenti e al calore, sono usati nei

detonatori per innescare l'esplosivo secondario. Gli esplosivi secondari,

tranne qualche eccezione, non sono sensibili alle sollecitazioni

meccaniche e termiche, e pertanto si possono definire stabili.

I detonatori possono essere innescati a loro volta tramite uno shock non

elettrico, i detonatori non elettrici si ottengono mediante l’accensione di

una miccia vincolata alla sua estremità o da un pezzo di corda di

detonazione, altresì da un detonatore elettrico costituito da un filamento

imbevuto in una miscela incendiaria, che è resa incandescente al

passaggio della corrente elettrica prodotta da una particolare batteria

definita esploditore[2].

Mentre i detonatori non elettrici sono istantanei, quelli elettrici possono

essere dotati di un dispositivo a tempo basato su un segnale elettronico

che agisce a distanza, prodotto ad esempio da un semplice telefono

cellulare, o basato su un periodo di innesco prefissato.

Le spolette a tempo sono per lo più utilizzate per la detonazione

d’ordigni militari (Bombe a Mano) [2].

Una gran varietà di detonatori è facilmente reperibile in commercio e

con varie dimensioni, tra cui le corde di detonazione dette anche micce

ordinarie a lenta combustione, utilizzate per l’innesco non elettrico dei

detonatori stessi o d’altri esplosivi a debita distanza e con sufficiente

ritardo di tempo. Queste sono realizzate con una polvere nera finissima

che non emette odori identificabili, avvolta con un filo in una pellicola o

in uno strato di plastica. Le stesse possono anche essere realizzate

sostituendo alla polvere nera un esplosivo secondario ad alta velocità di


28

detonazione, in questo caso si parla di “miccia detonante“, usata

soprattutto per la detonazione contemporanea di più cariche disposte a

distanza l’una dalle altre.

1.3 – Dispositivi Esplosivi E’ possibile individuare diverse tipologie di dispositivi esplosivi che si

contraddistinguono per innesco, cariche esplosive e contenitore.

Secondo i loro impieghi, tali ordigni si distinguono in tre diverse

categorie, Militari, Commerciali o Civili ed “Improvvisati“.

Gli Esplosivi Militari, come il trinitrotulene (TNT) e i vari esplosivi

plastici tra cui il C-4, possono essere usati singolarmente o mescolati tra

loro o con altre sostanze, tra cui la polvere d’alluminio, per migliorarne

le prestazioni. Lo stesso esplosivo plastico molto diffuso negli Stati Uniti

si ottiene proprio mescolando esplosivi secondari con sostanze plastiche

quali ad esempio i polimeri sintetici. Questi esplosivi sono caratterizzati

da una densità simile a quella dell’acqua, sono tolleranti alle condizioni

d’umidità e resistenti alle temperature estreme, si presentano sotto forme

e colori diversi fortemente dipendenti dai plastificanti e i coloranti

utilizzati durante il processo di produzione. Data la loro tolleranza agli

ambienti umidi possono essere facilmente nascosti in liquidi acquosi,

facendo attenzione ad evitare i solventi organici nei quali potrebbero

dissolversi.

L’innesco avviene singolarmente o tra più carichi esplosivi attraverso

l’uso di uno o più “blasting cap“ (cappuccio esplosivo) [2].

E’ stato provato che molti attacchi terroristici siano stati realizzati

proprio avvalendosi d’ordigni militari, con la possibilità di nasconderli in

solventi organici, quale il gasolio, a patto di isolare completamente


29

l’esplosivo dal solvente stesso; il beneficio che se ne trae è quello di

confondere l’odore dell’esplosivo rendendolo quindi meno intercettabile

dai controlli di sicurezza tenuti dalle unità cinofile.

Gli Esplosivi Commerciali d’uso civile si distinguono per dimensione,

colore e consistenza, tra solidi e gelatinosi, tra questi maggiormente

impiegati sono gli esplosivi da mina a base di Nitrato di Potassio o

d’Ammonio, e la gelatina esplosiva formata per oltre il 90% da

Nitroglicerina e la restante parte da Cotone Collodio. Questi ultimi si

annoverano tra i più pericolosi, in quanto i rischi sono connessi

soprattutto al trasporto e alla lavorazione, essendo, infatti, la

Nitroglicerina, un materiale esplosivo molto suscettibile agli urti [2].

Tali esplosivi hanno una densità molta variabile, usualmente simile a

quella dell’acqua, ciò li rende tolleranti al contatto con l’acqua stessa,

sebbene alcuni esplosivi commerciali non possono essere detonati in

completa immersione in acqua o in altri solventi. La detonazione avviene

con un singolo “blasting cap“ dato che l’innesco di un solo contenitore di

materiale esplosivo può provocare la detonazione di altri a lui adiacente

(Reazione a Catena) [2].

Gli esplosivi commerciali non hanno una lunga vita utile, soprattutto se

esposti per lunghi periodi ad elevate temperature, le quali possono

provocare la fuoriuscita d’oli esplosivi o d’altri elementi. Molti di questi

esplosivi, ed in particolar modo la Dinamite, emanano un forte odore,

l’esposizione a tali vapori può causare rapidamente un’intensa emicrania

della durata di alcuni minuti.


30

Gli Esplosivi “Improvvisati, sono realizzati clandestinamente e

artigianalmente, il loro potere distruttivo è funzione dei materiali

impiegati e delle abilità tecniche degli stessi produttori.

Tali ordigni presentano notevoli difficoltà, può, infatti, verificarsi una

detonazione prematura o incompleta, oltre che la fuoriuscita di un forte

vapore acido che può corrodere il contenitore metallico dell’esplosivo o

decomporsi quando esposto al calore. I contenitori degli esplosivi

maggiormente utilizzati sono tubi d’acciaio o di plastica con entrambe le

estremità tappate, solitamente di grandi dimensioni e per tanto facili

obiettivi delle tecniche d’ispezione, e per tale motivo sono usualmente

nascosti negli scompartimenti di un veicolo [2].

1.4 – Materiali Esplosivi

I materiali esplosivi, utilizzati per la fabbricazione delle bombe sono di

vario tipo. Questi possono essere distinti tra cariche concentrate e

allungate a seconda che siano ammassate in modo globulare o disposte in

tubi esplosivi. Possono essere interne o esterne a seconda che siano

disposte in cavità realizzate nel corpo da far esplodere o semplicemente

appoggiate su di essi. Si contraddistinguono in funzione delle loro

caratteristiche chimiche, fisiche, del colore e dell’odore. Variano

soprattutto da un luogo all’altro in funzione della loro reperibilità.

Alcune pubblicazioni scientifiche attestano che i migliori materiali

esplosivi rintracciabili negli Stati Uniti sono la Polvere Nera, utilizzata

per la produzione di bombe a tubo, la Dinamite, una miscela di Nitrato di

Ammonio e Oli Combustibili (ANFO), l’esplosivo al plastico (C-4) e il

TNT, questi ultimi due utilizzati soprattutto per operazioni militari. I più

usati esplosivi artificiali sono esplosivi chimici, che normalmente

comprendono una rapida e violenta reazione d’ossidazione che produce


31

una notevole quantità di gas ad alta temperatura. Sono esplosivi molto

versatili, compatti, disponibili in quantità e in ogni tipo.

Le caratteristiche chimiche e fisiche dei materiali esplosivi incidono

chiaramente sui loro vari impieghi, tra cui quelli militari e purtroppo

anche quelli terroristici. Tali caratteristiche influenzano la scelta operata

dai terroristi tra i vari materiali per la costruzione del proprio ordigno, al

fine di sottrarsi ai normali controlli di sicurezza all’interno degli

aeroporti o altri luoghi controllati dai metal detector, ed altri dispositivi

di sorveglianza.

Molti materiali esplosivi appaiono sotto forma di polvere bianca o nera,

solitamente con una struttura di tipo cristallina, ma quando caricate nel

contenitore queste sono generalmente bagnate, al fine di evitare

esplosioni premature, il che può chiaramente appesantirle e conferirgli

una configurazione pastosa.

E’ possibile distinguere diversi esplosivi utilizzati per la fabbricazione di

bombe, tra cui Esplosivi Semplici, come il TNT (trinitrotulene), la

NITROCELLULOSA e la NITROGLICERINA, ed Esplosivi Composti

quali la DINAMNITE, la POLVERE NERA, gli esplosivi plastici (C-4;

SEMTEX) e l’ANFO [2].

Esplosivi Semplici

• Il TNT è il più comune tra gli esplosivi utilizzati per la

produzione d’armi militari, si presenta sotto forma di polvere

cristallina o a fiocchi, di colore variabile tra il giallo e il marrone

acceso, se esposto all’ossigeno e ai raggi ultravioletti può subire

uno scolorimento e una riduzione di stabilità all’impatto, data la

sua composizione può essere sciolto per fargli assumere altre


32

forme, è assolutamente inodore e può essere utilizzato come

componente per molti altri esplosivi.

• NITROCELLULOSA è un materiale altamente infiammabile

costituito da cellulosa e acido nitrico. È più stabile e brucia più

velocemente della polvere nera emettendo gas molto caldi.

Esplode quando innescata, ed è solitamente impiegata per la

produzione di propellente o più semplicemente per il lancio di

proiettili in armi leggere; se asciutta è molto sensibile agli

impatti, all’attrito e alle scintille.

• NITROGLICERINA è un esplosivo molto potente tra i maggiori

elementi utilizzati per la produzione d’ordigni. Allo stato puro si

presenta come un liquido pesante di colore chiaro con la stessa

consistenza di un olio; può essere assorbita dal corpo mediante

inalazione o a contatto con la pelle inducendo stordimento e forte

emicrania, oltre ad emettere fumi molto tossici durante la sua

combustione.

• RDX, conosciuto anche come ciclonite, esogeno o T4. É una

nitroammina, ed è un materiale esplosivo ampiamente usato dai

militari. Si presenta come un solido cristallino di colore bianco.

E’ usato solitamente in miscele con altri esplosivi e plastificanti,

oppure desensibilizzanti. È un esplosivo stabile e si può

conservare per molto tempo se immagazzinato bene. E’

considerato il più potente tra tutti gli esplosivi militari. E’

prodotto facendo reagire l'acido nitrico concentrato

sull'esammina, nitrato d’ammonio, acido acetico, ed anidride


33

acetica. Comincia a decomporsi a circa 170°C, brucia piuttosto

che esplodere, per farlo esplodere necessita di un detonatore, ed è

molto sensibile quando cristallizzato a temperature inferiore ai -

4°C [4].

Esplosivi Composti

• DINAMITE è un esplosivo commerciale a differenza dei primi

tre, che sono invece d’uso militare, realizzata con materiale

poroso imbevuto con Nitroglicerina ed avvolto in un foglio di

cera, in un cartone o in una pellicola di plastica di colore

variabile secondo il produttore, solitamente confezionata come

un tubo di diverso diametro. La Nitroglicerina impiegata per la

sua realizzazione varia in percentuali comprese tra il 5% e il

90%, conferendogli un forte odore che può provocare immediati

stordimenti; è solitamente innescata con una corda di

detonazione o con un blasting cap.

• POLVERE NERA è un esplosivo chimico tra i più vecchi

conosciuti. Si realizza mediante una miscela di polveri tra cui

Nitrato di Potassio o Nitrato di Sodio, Carbone di Legna e Zolfo.

Il suo colore varia tra il nero e il marrone ed appare sottoforma di

polvere sottile o granulare. E’ solitamente arricchita con della

grafite che ha il compito di ridurre l’attrito tra i singoli granuli

che potrebbe anche provocarne l’innesco accidentale, dato che

questo materiale risulta molto sensibile all’attrito oltre che al

calore, all’impatto e alle scintille, ciò fa di questo materiale

esplosivo uno dei più pericolosi da maneggiare. Inoltre, risulta


34

essere sensibile all’elettricità, pertanto deve essere lavorato con

attrezzi di legno o di plastica a bassa conducibilità elettrica. Deve

essere conservata asciutta data la sua sensibilità all’acqua e

all’umidità che ne possono compromettere il corretto innesco.

• C-4, noto anche come esplosivo al plastico per la sua

composizione e modalità d’applicazione, è anch’esso un

esplosivo usato per scopi militari, in particolar modo come carico

per le demolizioni. E’ costituito da un materiale esplosivo (RDX)

arricchito con degli additivi che lo rendono simile ad uno stucco

facile da modellare. La peculiarità di quest’esplosivo è

chiaramente la sua flessibilità che lo rende facile da modellare e

da applicare su superfici irregolari, oltre che a un’elevata velocità

di detonazione.

• SEMTEX è un esplosivo plastico al pari del C-4 prodotto

originariamente nell’Europa dell’Est, è costituito da due

esplosivi il RDX e il PETN, arricchiti anch’essi con additivi che

lo rendono flessibile e malleabile, è assolutamente inodore e

caratterizzato da un colore che varia tra il giallo e il nero, preso

singolarmente è un materiale relativamente innocuo che può

facilmente essere manipolato. Per la sua detonazione si usa

solitamente un blasting cap o una corda di detonazione. La sua

vita utile si aggira intorno ai 10 anni oltre i quali iniziano a

perdere la sua flessibilità diventando friabile o rigido.

• ANFO è composto di una miscela di Nitrato d’Ammonio e Oli

Combustibili. Il Nitrato è solitamente di colore bianco, ha una


35

forma simile a quella di un fertilizzante ed è spesso imballato in

contenitori impermeabili che servono ad evitare che possa

perdere il suo potere e sensitività all’aumentare dell’umidità.

Alcune variazioni si possono ottenere mescolando il Nitrato

d’Ammonio con Enitrometano o con lo Zolfo. Ha una densità

pari a circa l’85% di quella dell’acqua, non è sensibile alle

temperature estreme e ciò nonostante deve essere tenuto asciutto

e non può quindi essere nascosto in acqua [4].

Si riporta di seguito, una tabella illustrativa dei parametri d’alcuni

esplosivi classici sopra elencati [3]:

Esplosivo

Volume specifico

d’Esplosione V [l/Kg]

Calore specifico

d’Esplosione Q [Kcal/Kg]

Temperatura d’Esplosione

T [°C]

Velocità di Detonazione

D [m/s]

Nitroglicerina 717 1470 4110 8000

Esogeno 908 1500 3850 8300

TNT 728 1000 2950 6800

Fulminato di Mercurio 304 368 4810 5400

Dinamite (62%) 634 1200 4040 6600

Tabella 1. 1 – Parametri che caratterizzano fisicamente alcuni materiali esplosivi

CCAAPPIITTOOLLOO 22 -- DDEESSCCRRIIZZIIOONNEE DDII UUNN TTEESSTT DD’’EESSPPLLOOSSIIOONNEE SSUU UUNNAA PPIIAASSTTRRAA IINN CC..AA..

36

Capitolo 2 Descrizione di un test d’esplosione su una piastra in c.a.

Oggetto della tesi è stata la modellazione numerica di una piastra in c.a.

sottoposta all’azione di un carico dinamico prodotto da un’esplosione.

A tale fine, è stato necessario avvalersi dei risultati di una prova

sperimentale, condotta presso il Dipartimento d’Ingegneria Civile,

Architettonica ed Ambientale, dell’Università del Missouri-Rolla (USA).

Realizzata per opera degli Ingegneri Pedro F. Silva, Binggeng. Lu &

Antonio Nanni [5].

Nel caso specifico si tratta di una condizione d’esplosione (blast), di cui

si vogliono descrivere gli effetti prodotti in termini di pressione,

attraverso l’uso di una trattazione numerica di seguito riportata.

Operando il confronto tra i risultati noti a priori del test sperimentale, e

quelli della modellazione numerica, si ha la possibilità di verificare

l’accuratezza della simulazione, che sarà tanto più attendibile, quanto

maggiore sarà la fedeltà nella riproduzione della prova, e nel

conseguimento di risultati simili a quelli ottenuti empiricamente.

La simulazione numerica, è stata realizzata utilizzando il software di

calcolo LS-DYNA Version 970 [11]. Il programma, è stato impiegato con

il fine di implementare un metodo di calcolo numerico, per la

discretizzazione prima e la risoluzione analitica poi, dell’elemento

strutturale.

Tra i possibili metodi a disposizione, si è scelto di implementare un

Metodo Particellare di tipo Meshfree, in particolare lo Smoothed

Particle Hydrodynamics (SPH) [7]. Utilizzando lo stesso software, si è

potuto ripetere la stessa simulazione, implementando un altro metodo di

calcolo, il più tradizionale Metodo degli Elementi Finiti (FEM) [6].


37

In seguito, si è operato un confronto tra i risultati ottenuti, con i due

diversi metodi impiegati, al fine di capire quale tra questi meglio si

presta per la modellazione di un elemento strutturale sottoposto a carico

da esplosione.

2.1 – Descrizione della Prova Sperimentale Il test si è svolto sottoponendo all’azione di un carico dinamico, una

piastra in calcestruzzo armato. La sollecitazione dinamica è stata

prodotta dalla detonazione di una carica esplosiva di peso W, sospesa in

aria alla distanza R dal centro dell’elemento. Il materiale esplosivo

utilizzato è l’RDX (Esogeno), considerato tra i più potenti esplosivi d’uso

militare.

Di seguito si descrivono:

• Caratteristiche Geometriche della piastra in c.a.

Le dimensioni in pianta sono 48x48 in (1,2x1,2 m), mentre lo

spessore è di 3,5 in (0,0875 m). Le armature metalliche

bidirezionali, sono disposte sul fondo della piastra, con copriferro

di 0.5 in (0,0125 m). Le barre hanno sezioni trasversali di 0.11 in2

(φ9), e interasse di 6 in (0,15 m). Si ha quindi una percentuale

d’armatura in entrambi le direzioni di 0,528%.

• Condizioni di vincolo;

I campioni di prova sono stati appoggiati su due travi parallele

d’acciaio, con sezione a doppia T, le cui flangie avevano una

larghezza di 6 in (0,15 m). Si possono osservare su entrambe le

travi d’acciaio due file di spinotti metallici, a cui la piastra è stata


38

vincolata, e che formano in questo modo un incastro su entrambi i

bordi esterni.

Sezione A-A ; B -B

1,2

0,15

0,09

AA

y

x

B

φ9/15 cm

9φ9

B

1,2

0,15

Figura 2. 1 – Caratteristiche Geometriche della piastra, e dell’armatura metallica.

Figura 2. 2 – Setup di Prova.


39

Come si può notare, la piastra oltre ad essere geometricamente

simmetrica, rispetto ad una coppia di assi cartesiani con origine nel

centro, lo è anche nei confronti delle condizioni di vincolo e di carico,

essendo la carica esplosiva applicata nel centro della piastra, ed essendo

questa ultima semplicemente appoggiata su due lati.

Di seguito si riportano, le proprietà dei materiali usati per la prova

sperimentale [5]:

Materiali Strain-rate [sec-1]Modulo di

Elasticità di Young E [Kg/cm2]

Rapporto di Poisson ν

Acciaio 100 2,90E+07 0,3

Calcestruzzo 100 3,78E+06 0,2

Tabella 2. 1 – Proprietà elastiche dei materiali impiegati per la costruzione della piastra, allo Strain-Rate di 100 sec-1.

Tensione di Rottura a

Compressione [Kg/cm2]


Trazione [Kg/cm2]

Deformazione Ultima εu [‰]

Calcestruzzo

717,13 195,8 4,76

Tabella 2. 2 – Proprietà dinamiche del calcestruzzo allo Strain-Rate di 100 sec-1.

Strain-rate [sec-1]Tensione di

Snervamento [Kg/cm2]

1,00E-04 3958,28 Acciaio

174 5582,37

Tabella 2. 3 – Tensione di Snervamento dell’acciaio, per due diversi valori di Strain-Rate.


40

Mentre per il calcestruzzo sono fornite le sole proprietà dinamiche

(Strain-Rate di 100 sec-1), per l’acciaio, si fa distinzione tra la tensione

di snervamento calcolata per condizioni di carico quasi statiche (Strain-

Rate 1,00E-04 sec-1), e per le condizioni di carico dinamiche (Strain-

Rate 174 sec-1). Da questi dati si nota un incremento della tensione di

snervamento del materiale, all’aumentare della velocità di deformazione

ε [sec-1], mentre il modulo di elasticità si mantiene costante.

Note le caratteristiche meccaniche del calcestruzzo per condizioni di

carico dinamiche, sono state da queste derivate le stesse caratteristiche in

condizioni quasi statiche. Utilizzando la formulazione proposta dal

bollettino CEB-FIP Model Code 1990 [20].

Ottenendo i seguenti risultati:

Materiali

Modulo d’Elasticità di

Young E [Kg/cm2]

Modulo d’Elasticità

Tangenziale G [Kg/cm2]

Calcestruzzo 179843,71 74934,88


Compressione [Kg/cm2]

Tensione di Rottura a Trazione [Kg/cm2]

Calcestruzzo

212,65 28,67

Tabella 2. 4 – Proprietà del calcestruzzo, per velocità di deformazione quasi statica.


41

2.2 – Evidenze Sperimentali Si riportano ora i risultati delle prove eseguite, in termini di danno

prodotto [5]:

Test Peso della carica esplosiva RDX.

W [Kg]

Distanza della carica dal

centro della piastra. R [m]

Livello di Danno prodotto

1 0,10 0,91 I (No Cracking)

2 0,50 0,91 II (Minor Cracking)

3 0,90 0,30 III (Major Cracking)

4 1,71 0,20 IV (Severe Damage)

Tabella 2. 5 – Condizioni di carico ed Evidenze Sperimentali ad esse associate.

Dalle evidenze sperimentali si evince, che all’aumentare dell’entità del

carico, e contemporaneamente al ridursi della distanza tra la piastra e il

centro d’esplosione, si osservano danni crescenti in termini di

fessurazione, fino a giungere alla condizione di rottura dell’elemento.

Nelle immagini successive, sono raffigurati i risultati delle prove 3 e 4:

Figura 2. 3 – Immagine della piastra in seguito alla terza Prova.


42

Figura 2. 4 – Immagini della piastra in seguito alla quarta Prova.

Si osserva che:

• Al seguito della prima prova non si riscontrano danni percepibili;

• Nel secondo caso si osservano piccole lesioni d’entità in ogni

modo non apprezzabili;

• Nel terzo caso la piastra presenta fessure rilevanti. Il quadro

fessurativo che si presenta in questo caso è caratterizzato da

lesioni radiali, che dal centro della piastra, ovviamente

maggiormente sollecitata perché nelle immediate vicinanze della

carica esplosiva, si estendono verso i bordi della stessa, con

ampiezze decrescenti linearmente. La lesione maggiore nel centro

della piastra misura un’ampiezza di circa 3 mm. La deformazione

dell’elemento strutturale è data dalla somma di due aliquote, una

deformazione elastica e una plastica. La deformazione plastica si

manifesta mediante uno spostamento residuo, che si può misurare

al termine della prova. La piastra presenta nel centro un

abbassamento residuo di 21 mm, e di 11 mm in prossimità della

mezzeria del bordo esterno.

• Nella quarta ed ultima prova, aumentando il carico esplosivo e

riducendo ulteriormente la distanza della carica esplosiva dal


43

centro della piastra, si provoca la rottura della stessa, così come

mostrato nella precedente immagine.

2.3 – Finalità della Prova Sperimentale Il test di prova è stato condotto con l’intenzione di verificare la validità

del metodo DBD (Displacement Based Design method), usato nel

valutare la dimensione del carico esplosivo W, e la distanza di tolleranza

R, tali da produrre un livello di danno stabilito a priori [5].

Oggetto della ricerca era l’individuazione di uno strumento utile nel

determinare, entro un grado di precisione ragionevole, il livello di

vulnerabilità delle strutture esistenti. Inoltre, nel modellare un’esplosione

è di fondamentale importanza definire un livello di danno tollerabile. Gli

effetti prodotti da un’esplosione possono variare dal danno minore fino

al completo collasso della struttura, con le conseguenti e considerevoli

perdite di vite umane. Questo programma di ricerca ha mostrato risultati

promettenti nell’impiego del metodo per predire il blast load (carico

esplosivo), in termini di distanza e peso della carica esplosiva di

sicurezza.

Le evidenze sperimentali hanno mostrato che i livelli di danno realizzati

erano simili ai valori predetti.

Dato che le finalità della prova sperimentale e della tesi differiscono tra

loro, si preferisce omettere la descrizione del metodo DBD, non essendo

questo di nostro interesse. Sarà quindi sufficiente avvalersi dei dati dei

test eseguiti e dei risultati da loro conseguiti, da utilizzare per la

modellazione degli stessi.

CCAAPPIITTOOLLOO 33 –– MMOODDEELLLLAAZZIIOONNEE DDEELLLLEE AAZZIIOONNII PPRROODDOOTTTTEE DDAALLLL’’EESSPPLLOOSSIIOONNEE

44

Capitolo 3 Modellazione delle azioni prodotte dall’esplosione.

Al fine di valutare le azioni prodotte sulla superficie della piastra

dall’esplosione, sarà necessario definire per ogni prova, la curva di

pressione a cui la piastra sarà sottoposta.

La curva di pressione definisce al variare del tempo, in un dato punto

dello spazio alla distanza ridotta Z dal centro della carica, l’andamento

della pressione generata dall’esplosione. Di questa curva di carico, è

possibile distinguere una fase positiva ed una negativa. La fase positiva è

dovuta all’incremento di pressione rispetto alla normale pressione

atmosferica, mentre quella negativa si deve ad una depressione, in virtù

del ritorno verso il centro di esplosione dell’aria che è stata spostata [3].

Nella modellazione effettuata dagli autori della prova sperimentale,

l’onda di pressione è stata semplificata con una curva di carico

triangolare. I parametri principali richiesti per definire la blast load sono

il Picco di Pressione ΔPs, e la durata dell'impulso di blast td, o durata

della fase positiva. Assumendo nullo il tempo d’arrivo ta del fronte

d’onda sulla superficie piana, l’andamento qualitativo della curva di

pressione sarà come mostrato nella successiva immagine [5]:

td t

P

ΔPs

Figura 3. 1 – Curva di Pressione semplificata, usata dagli autori per la modellazione della prova.


45

Semplici espressioni possono essere usate per riferire ΔPs e td al peso

della carica esplosiva e alla distanza, indicate rispettivamente con W e R.

Il picco di pressione ΔPs, è valutato come una funzione inversamente

proporzionale della distanza ridotta Z, secondo la formulazione proposta

da H.L. Brode [5].

Al variare del campo di validità, ΔPs sarà definito come:

Dove, la distanza ridotta Z si esprime in funzione di R e W così come

segue:

La durata della fase positiva, o impulso dell’esplosione td, si definisce

sempre in funzione di R ed W, a partire dalle seguenti funzioni

logaritmiche [5]:

La curva di pressione così calcolata nel solo centro della piastra, è poi

estesa all’intera superficie di carico.

232

23

/101,0019,085,5455,1975,0)(

/1017,6)(

cmKgpZZZ

Zp

cmKgpZ

Zp

≤Δ≤−++=Δ

≥Δ+=Δ

3 WRZ =

0,1log95,175,2log310310 ≤⎟

⎠

⎞⎜⎝

⎛⋅+−≅⎟⎠

⎞⎜⎝

⎛ ZWR

Wt d

0,1log27,075,2log310310 ≥⎟

⎠

⎞⎜⎝

⎛⋅+−≅⎟⎠

⎞⎜⎝

⎛ ZWR

Wt d


46

Nella nostra modellazione, sia il picco di pressione sia la durata della

fase positiva sono stati valutati in maniera diversa. Per la definizione di

una curva di pressione, che si avvicini quanto più verosimilmente alla

descrizione del reale fenomeno fisico provocato dall’esplosione, si è

preferito seguire un diverso metodo analitico, la cui trattazione teorica è

stata esplicitata nel Capitolo 1.

3.1 – Curva di pressione adottata Essendo la piastra simmetrica, sia geometricamente, che per quanto

concerne le condizioni di vincolo e di carico, sarà sufficiente definire la

curva di pressione per un solo quarto della piastra.

Per conseguire una maggiore accuratezza nella modellazione della prova,

la piastra viene discretizzata in un numero finito d’elementi, d’uguali

caratteristiche fisiche, geometriche e meccaniche. Si suddividerà quindi,

la superficie di carico della piastra in elementi di dimensione 10x10 cm,

sia lungo l’asse x sia y, del sistema di riferimento cartesiano introdotto.

6

5

4

3

2

1

FEDCBA

x

yABCDEF

F E D C B A A B C D E F

1

2

3

4

5

6

1

2

3

4

5

6

6

5

4

3

2

1

Figura 3. 2 – Discretizzazione della superficie di carico.


47

Limitando il calcolo al solo quadrante positivo del sistema di riferimento

scelto, di origine coincidente con il baricentro della piastra, si definisce

una curva di pressione per ogni superficie elementare. Data la

discretizzazione adottata, e le dimensioni dei singoli elementi, saranno

36 le curve di pressione da definire, per un totale di 144 curve di carico

per ognuna delle 4 prove sperimentali. Questa semplificazione, è resa

possibile dalla simmetria d’applicazione del carico esplosivo, essendo

questo centrato rispetto al baricentro della piastra stessa. Inoltre, la

propagazione dell’onda di sovrapressione avviene in maniera sferica dal

centro d’esplosione, essendo questo sollevato rispetto alla piastra e non

ostacolato se non dalla piastra stessa, che n’è investita. La pressione sui

singoli elementi, e in un dato istante di tempo, sarà quindi una funzione

della sola distanza effettiva dal centro d’esplosione. La curva di

pressione quindi si ripete identicamente sulle superfici di 4 elementi

simmetrici rispetto al centro della piastra.

Si descrive in seguito le operazioni successive, che si è scelto di adottare

per definire la curva di pressione prodotta dall’esplosione, facendo

riferimento all’i-esimo elemento di dimensioni finite, in cui è stata

suddivisa l’area di carico.

Coefficiente d’equivalenza

Nel considerare il peso della carica esplosiva è necessario introdurre un

coefficiente d’equivalenza che consente di esprimere il peso della carica

esplosiva, usata nella prova sperimentale, in funzione del TNT

(Trinitrotulene). L’esigenza di equiparare il peso di una qualunque carica

esplosiva al TNT, nasce dalla necessità di impiegare formule empiriche

nel calcolo della curva di pressione, ricavate attraverso prove


48

sperimentali in cui il TNT è stato usato come materiale esplosivo di

riferimento.

Il coefficiente d’equivalenza q è dato dal rapporto tra i Calori Specifici

d’Esplosione Qw, dei due materiali [3].

Dal coefficiente d’equivalenza si definisce il peso equivalente in TNT

della carica esplosiva usata nel test sperimentale:

Nel caso specifico si avrà:

Caratteristiche dell'esplosivo

Esplosivo Qw [Kcal/Kg]

RDX 1500 TNT 1000

Coefficiente d’equivalenza q 1,5

Tabella 3. 1 – Calore specifico degli esplosivi e Coefficiente d’Equivalenza.

Ciò significa che, gli effetti prodotti da 1 Kg di RDX sono equivalenti a

quelli prodotti da 1,5 Kg di TNT.

Distanza Ridotta Z

Noti i pesi equivalenti della carica esplosiva WTNT e le distanze effettive

tra il centro d’esplosione e quello della superficie dell’i-esimo elemento

Ri, è possibile calcolare la Distanza Ridotta Zi d’ogni elemento come [3].

TNT

RDX

QQ

q =

RDXTNT WqW ⋅=


49

Dove la distanza effettiva dal centro d’esplosione Ri è così definita:

Così come si può semplicemente dedurre dalla seguente immagine, in

cui si rappresenta la pianta e la sezione di un quarto di piastra, e l’i-

esimo elemento a lei appartenente:

W

bi

a

Ri

i

y

x

A B C D E F

1

2

3

4

5

6

Pianta della piastra

Sezione Trasversale

Figura 3. 3 – Distanza effettiva dell’i-esimo elemento dal centro d’esplosione.

3TNT

ii W

RZ =

22ii baR +=


50

Picco di Pressione Δp(Z)

In funzione delle Distanze Ridotte, si definisce il picco di pressione

Δp(Z), per ogni prova sperimentale e per singolo elemento in cui la

piastra è stata in precedenza discretizzata.

Si è scelto di definire Δp(Z) come la media dei valori forniti dalle varie

formulazioni in precedenza esposte (Brode H.L. Henrich J.; Petrowsky

G.I.; Sadowsky M.A.) [3], al variare della distanza ridotta, e quindi del

relativo campo di validità, a differenza di quanto fatto nella definizione

della curva di pressione semplificata, in cui si è scelto di adottare una

sola delle precedenti formulazioni.

Coefficiente Sperimentale α.

Per definire la Curva di Pressione, è necessario esprimere in funzione

della distanza ridotta Z, un coefficiente di natura sperimentale α.

Tale coefficiente si può trovare tabellato [3], in funzione però di soli

alcuni valori di Z. Mediante un’operazione d’estrapolazione, dai valori

noti del coefficiente, è stato possibile ricavare una funzione di Z, che ci

consentiva di trovare i diversi valori di α al variare della distanza ridotta.


51

Estrapolazioney = 2,651x-0,7694

R2 = 0,9539

0

0,5

1

1,5

2

2,5

3

3,5

0 20 40 60

Z

α αPotenza (a)

Grafico 3. 1 - Estrapolazione numerica del coefficiente sperimentale α.

Durata della Fase Positiva t+

È l’intervallo di tempo entro cui si annulla la sovrapressione prodotta

dall’esplosione nell’aria circostante. L’espressione usata per la sua

definizione è la seguente [3]:

Allo stesso modo di quanto visto nella definizione della curva di

pressione semplificata, anche in questo caso la durata della fase positiva

si esprime in funzione di W e di R, ma attraverso una più semplice

espressione semiempirica.

iTNT RWBt 6310 −+ ⋅=


52

Curva di Pressione P(Z, t)

La curva di pressione non è altro che l’andamento della pressione

prodotta dall’esplosione in un dato istante di tempo t, e ad una certa

distanza ridotta Z. Analiticamente l’andamento della pressione può

essere calcolato come segue [3]:

Fissata la distanza ridotta Z, da cui dipendono il picco di pressione

Δp(Z), e il coefficiente sperimentale α, attraverso quest’espressione si

ottiene la curva di pressione al variare del tempo. Essendo Z a sua volta

funzione della distanza effettiva dal centro d’esplosione R, al variare di

questa otteniamo diversi valori di Z, e quindi diverse curve di pressione.

Ragion per cui sarà necessaria definire per ogni elemento della piastra

una curva di pressione diversa.

Quindi, per semplicità d’esposizione si riportano di seguito i risultati

conseguiti nella definizione di una singola curva di pressione, in

particolare quella calcolata nel centro della piastra. Essendo, infatti, il

centro della superficie superiore della piastra, il punto più vicino al

centro d’esplosione, sarà anche quello maggiormente sollecitato.

PROVA DATI

1 2 3 4

a (m) 0,91 0,91 0,30 0,20

bi (m) 0,000 0,000 0,000 0,000

Ri (m) 0,91 0,91 0,30 0,20 WRDX (Kg) 0,10 0,50 0,90 1,71

WTNT (Kg) 0,15 0,75 1,35 2,57 Zi 1,71 1,00 0,27 0,15

Tabella 3. 2 - Distanza Effettiva Ri, Carico Equivalente WTNT, Distanza Ridotta Zi.

+−

+ ⎟⎠⎞

⎜⎝⎛ −Δ= t

t

ettZptZp

α1)(),(


53

Δp(Z)[Kg/cm2]

PROVA Autore

1 2 3 4

Campo di validità

2,33 7,67 336,00 2149,19 Δp>10 Kg/cm2 Brode H.L.

2,21 8,23 315,82 1950,47 0,1<Δp<10 Kg/cm2

10,03 19,22 110,31 255,00 0,05<Z<0,3

3,93 7,98 125,01 710,81 0,3<Z<1 Henrich J.

2,42 7,97 221,81 1248,47 Z>1 Petrowsky

G.I. 1,13 9,65 534,00 3429,69 Z<1

Sadowsky M.A. 2,61 9,77 362,41 2208,72 1<Z<15

Tabella 3. 3 - Picco di Pressione Δp(Z), calcolato secondo le varie formulazioni, e per le singole prove sperimentali.

I valori segnati in grassetto sono quelli compatibili con il campo di

validità, per tanto impiegati nel calcolo del valore medio del Picco di

pressione, per ogni singola prova.

Prova 1 2 3 4

α 1,75 2,65 7,26 11,41

Δp(Z)[Kg/cm2] 2,41 8,66 326,77 1944,62

Tabella 3. 4 - Picco di Pressione medio, Coefficiente α.

PROVA Durata fase positiva 1 2 3 4

t+ (s) 0,00090 0,00118 0,00075 0,00068

t+ (ms) 0,904 1,182 0,749 0,680 B 1,3

Tabella 3. 5 - Durata della fase positiva t+, e costante numerica B.


54

Curva di pressione calcolata nel centro della piastra. Prova 1 Prova 2 Prova 3 Prova 4

t (ms) P(Z,t)[Kg/cm2] P(Z,t)[Kg/cm2] P(Z,t)[Kg/cm2] P(Z,t)[Kg/cm2]

0,00 2,41 8,66 326,77 1944,62 0,05 2,07 7,41 187,77 778,80 0,10 1,77 6,33 107,34 309,94 0,15 1,51 5,40 61,00 122,43 0,20 1,28 4,59 34,42 47,93 0,25 1,08 3,90 19,26 18,56 0,30 0,90 3,30 10,67 7,09 0,35 0,75 2,78 5,84 2,66 0,40 0,62 2,34 3,14 0,98 0,45 0,51 1,95 1,66 0,35 0,50 0,41 1,63 0,85 0,12 0,55 0,33 1,35 0,42 0,04 0,60 0,25 1,11 0,19 0,01 0,65 0,19 0,91 0,08 0,00 0,70 0,14 0,73 0,02 0,00 0,75 0,10 0,59 0,00 0,00 0,80 0,06 0,47 -0,01 0,00 0,85 0,03 0,36 -0,01 0,00 0,90 0,00 0,27 -0,01 0,00 0,95 -0,02 0,20 -0,01 0,00 1,00 -0,04 0,14 -0,01 0,00 1,05 -0,05 0,09 0,00 0,00 1,10 -0,06 0,05 0,00 0,00 1,15 -0,07 0,02 0,00 0,00 1,20 -0,08 -0,01 0,00 0,00 1,25 -0,08 -0,03 0,00 0,00 1,30 -0,09 -0,05 0,00 0,00 1,35 -0,09 -0,06 0,00 0,00 1,40 -0,09 -0,07 0,00 0,00 1,45 -0,09 -0,08 0,00 0,00 1,50 -0,09 -0,08 0,00 0,00 1,60 -0,08 -0,08 0,00 0,00 1,70 -0,08 -0,08 0,00 0,00 1,80 -0,07 -0,08 0,00 0,00 1,90 -0,07 -0,07 0,00 0,00 2,05 -0,06 -0,06 0,00 0,00 2,20 -0,05 -0,05 0,00 0,00 2,35 -0,04 -0,04 0,00 0,00 2,50 -0,03 -0,04 0,00 0,00 2,65 -0,03 -0,03 0,00 0,00 2,80 -0,02 -0,02 0,00 0,00 3,00 -0,02 -0,02 0,00 0,00 3,20 -0,01 -0,01 0,00 0,00 3,40 -0,01 -0,01 0,00 0,00 3,60 -0,01 -0,01 0,00 0,00 3,80 0,00 0,00 0,00 0,00

Tabella 3. 6 - Curve di Pressione P(Z, t), calcolate per le quattro prove sperimentali.


55

Nel definire l’andamento della pressione al variare del tempo, si trascura

per semplicità di calcolo, la determinazione del tempo d’arrivo ta,

essendo irrilevante rispetto all’intera durata della curva di pressione. Il

tempo d’arrivo si definisce come il tempo impiegato dal fronte d’onda ad

impattare contro la superficie superiore della piastra. In realtà le curve di

pressione dovrebbero essere traslate di una quantità pari a ta lungo l’asse

delle ascisse.

Il metodo sopra esposto sarà impiegato nella definizione delle curve di

pressione, per ogni singolo elemento della piastra, al variare della

distanza effettiva tra il centro d’esplosione e il baricentro dell’elemento,

e per ognuna delle quattro prove sperimentali di cui si conoscono i dati

iniziali.

Di seguito si riportano le curve di pressione calcolate nel centro della

piastra, e per ogni prova sperimentale.

56

Curva delle Pressioni Prova_1

-1

0

1

2

3

0,00 0,15 0,30 0,45 0,60 0,75 0,90 1,05 1,20 1,35 1,50 1,65 1,80 1,95 2,10 2,25 2,40 2,55 2,70 2,85 3,00 3,15 3,30 3,45 3,60 3,75 3,90

t [ms]

P [Kg/cm2] Prova_1

+−

+ ⎟⎠⎞

⎜⎝⎛ −Δ= t

t

ettZptZP

α1)(),(

57


-1

0

1

2

3

4

5

6

7

8

9

10

0,00 0,15 0,30 0,45 0,60 0,75 0,90 1,05 1,20 1,35 1,50 1,65 1,80 1,95 2,10 2,25 2,40 2,55 2,70 2,85 3,00 3,15 3,30 3,45 3,60 3,75 3,90

t [ms]

P [Kg/cm2]

+−

+ ⎟⎠⎞

⎜⎝⎛ −Δ= t

t

ettZptZP

α1)(),(

58


-10

10

30

50

70

90

110

130

150

170

190

210

230

250

270

290

310

330

0,00 0,15 0,30 0,45 0,60 0,75 0,90 1,05 1,20 1,35 1,50 1,65 1,80 1,95 2,10 2,25 2,40 2,55 2,70 2,85 3,00 3,15 3,30 3,45 3,60 3,75 3,90

t [ms]

P [Kg/cm2]

+−

+ ⎟⎠⎞

⎜⎝⎛ −Δ= t

t

ettZptZP

α1)(),(

59


-100

50

200

350

500

650

800

950

1100

1250

1400

1550

1700

1850

0,00 0,15 0,30 0,45 0,60 0,75 0,90 1,05 1,20 1,35 1,50 1,65 1,80 1,95 2,10 2,25 2,40 2,55 2,70 2,85 3,00 3,15 3,30 3,45 3,60 3,75 3,90

t [ms]

P [Kg/cm2]

+−

+ ⎟⎠⎞

⎜⎝⎛ −Δ= t

t

ettZptZP

α1)(),(

CCAAPPIITTOOLLOO 44 –– MMEETTOODDII DDII DDIISSCCRREETTIIZZZZAAZZIIOONNEE

60

Capitolo 4 Metodi di Discretizzazione

I metodi di discretizzazione che si sono scelti di impiegare per la

simulazione numerica, sono il Finite Element Method (FEM) e lo

Smoothed Particle Hydrodynamics (SPH). Al termine della

modellazione, avendo utilizzato entrambi i metodi per le stesse prove

sperimentali, si avrà la possibilità di operare uno studio parametrico tra i

due. Si potrà capire quale tra questi meglio si presta per la modellazione

numerica di un elemento strutturale sottoposto a carichi da esplosione.

Questo capitolo è dedicato ad una breve descrizione del più noto ed

utilizzato FEM, a prescindere dal quale non sarebbe possibile introdurre i

Meshfree Particle Methods. Infatti, questi ultimi sono stati sviluppati con

l’intento di estendere i primi, ampliandone il raggio d’applicazione, e

rendendoli più versatili.

In seguito si dedicheranno maggiori attenzioni ai Meshfree Particle

Methods, facendo riferimento in particolare allo Smoothed Particle

Hydrodynamics (SPH), impiegato nella modellazione.

4.1 – Finite Element Method (FEM) Sin dalla sua invenzione, risalente alla metà del secolo scorso, il metodo

degli elementi finiti (FEM), è stato largamente impiegato per le

computazioni numeriche nell’ambito ingegneristico [6].

Il FEM rientra tra i metodi approssimati per la risoluzione delle

equazioni fondamentali della statica e della dinamica.

La sua principale caratteristica risiede nella possibilità di suddividere un

continuo in un numero discreto d’elementi, che non si sovrappongono


61

nel dominio di nostro interesse. Tale operazione è così detta di

Discretizzazione del continuo in Elementi Finiti.

I singoli elementi sono vincolati tra loro da una mappa topologica,

definita Mesh (Maglia), sopra la quale sono costruite le funzioni di

interpolazione [6].

Ad ogni elemento in cui il continuo è stato suddiviso, vanno associate

delle funzioni spostamento mediante il vettore spostamento così definito:

}{ tzyxss ,,,=

Questi vettori sono definiti attraverso la funzione di spostamento dei

singoli nodi [6]:

( )( )( )( )

[ ]{ }Δ=⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧= N

tzyxwtzyxvtzyxu

tzyxs,,,,,,,,,

,,,

Dove [N] è la matrice delle funzioni di forma, questa lega gli

spostamenti generici all’interno dell’i-esimo elemento agli spostamenti

dei nodi {Δ}, che lo individuano nello spazio. Le funzioni di forma sono

in numero pari al numero di nodi della Mesh, e per definizione valgono 1

se calcolate nel nodo a cui si riferiscono, e 0 negli altri nodi

dell’elemento i-esimo. Le funzioni di forma hanno un ruolo

fondamentale nel metodo d’analisi, in quanto rappresentano il polinomio

che approssima il valore esatto degli spostamenti, e determinano quindi

la convergenza del metodo [6].

La Mesh assicura la compatibilità delle interpolazioni, sebbene non

sempre le condizioni di compatibilità numerica coincidono con quelle

fisiche del continuo stesso.


62

In oltre, quando uno stesso corpo presenta caratteristiche locali diverse, o

geometrie complesse, è necessario suddividerlo in un numero molto

elevato di elementi finiti, realizzando quindi una maglia sempre più fitta,

ed aumentando l’onere di computazione.

Il FEM è quindi un metodo limitato, che necessità di modifiche e

adeguamenti tra i suoi limiti si annoverano [6]:

• La compatibilità numerica è sola un’approssimazione della reale

compatibilità fisica del continuo. Questa approssimazione non

sempre risulta essere precisa ed accurata. Basti pensare ad

un’eccessiva distorsione della Mesh, che può provare il termine

della computazione o in ogni modo un errore rilevante nella

stessa;

• Campi d’interpolazione d’elevati ordini di grandezza sono difficili

da riprodurre per domini con geometrie arbitrarie e

pluridimensionali. Ciò comporta una restrizione nell’applicare i

metodi di computazione basati sugli elementi finiti, per la

risoluzione di problemi pratici, come ad esempio per la

simulazione di Plates e Shell (Piatti e Gusci), con alti gradienti

d’elasticità e plasticità;

• Quando s’impiegano contemporaneamente Mesh differenti per la

simulazione numerica di uno stesso fenomeno fisico, è possibile

ottenere modelli differenti di disintegrazione per uno stesso

materiale. Tale fenomeno è noto come “suscettibilità” del FEM

nei confronti della rottura graduale del materiale;


63

• In alcune applicazioni, la Mesh stessa può essere fonte di disturbo

nella computazione numerica. Un esempio ben noto è la

simulazione del problema di localizzazione degli sforzi;

• Questi Metodi tradizionali, presentano notevoli difficoltà

d’applicazione per problemi che presentano deformabilità dei

vincoli, superfici indipendenti, o interfacce in movimento;

• Le funzioni di spostamento non sempre assicurano la continuità

degli spostamenti richiesta tra gli elementi adiacenti, può essere

quindi violata la condizione di congruenza.

Per questi ed altri motivi l’operazione di discretizzazione con il metodo

degli elementi finiti, non sempre risulta essere vantaggiosa nelle

computazioni numeriche.

Si rende quindi necessario un adeguamento del FEM, o almeno una

ricerca di metodi alternativi da utilizzare in particolari applicazioni.

Una prima soluzione è stata individuata nello sviluppo del così detto

“Arbitrary Lagrangian Eulerian” (ALE).

L’obiettivo di tale formulazione risiede nel rendere indipendente le

caratteristiche del materiale dalla Mesh, così che la sua distorsione può

essere minimizzata.

Purtroppo, in alcune formulazioni con l’ALE, quali ad esempio

modellazioni di fluidi o di grandi deformazioni del continuo, la

distorsione della maglia è ancora presente, e ciò provoca errori nella

computazione numerica.

Inoltre, gli effetti di trasporto convettivi in ALE, spesso conducono ad

oscillazioni spurie che hanno bisogno di essere stabilizzate da una

diffusione artificiale o una stabilizzazione di Petrov-Galerkin.


64

Per superare le difficoltà e i limiti sopra elencati, sarà quindi necessario

trovare un nuovo metodo di discretizzazione del continuo.

E’ in questa ottica che nascono i Meshfree Particle Methods (Metodi di

Discretizzazione Particellari), con l’obiettivo di superare i limiti dei

metodi FEM.

I Meshfree Particle Methods, in quanto metodi numerici, presentano

alcune caratteristiche in comune con altri metodi numerici, quali il già

citato Metodo degli Elementi Finiti, o il Metodo delle Differenze Finite.

Tuttavia si differenziano gli uni dagli altri in virtù dei processi

d’implementazione.

Attraverso questi metodi, è computazionalmente efficace e fisicamente

accurato, discretizzare un continuo attraverso un set di punti nodali, o

particelle distribuite nello spazio in maniera del tutto arbitraria, senza in

pratica costrizioni di maglia, che provveda alla loro connessione.

I vantaggi di tali metodi possono essere così elencati [7].

• E’ possibile trattare facilmente grandi deformazioni, essendo le

connessioni tra i nodi generate come una parte della

computazione stessa, che può subire variazioni nel tempo;

• Il metodo interpreta facilmente il danno degli elementi del

modello, tale capacità si dimostra molto utile per la modellazione

della rottura graduale;

• I nodi possono essere aggiunti piuttosto facilmente, cosi com’è

possibile realizzare modifiche delle caratteristiche geometriche

del continuo modellato;

• E’ possibile usare modelli pluri-dimensionali del continuo, al fine

di modellare grandi deformazioni di strutture a guscio sottili, come

i nanotubi;


65

• Mediante la discretizzazione Meshfree si possono rappresentare

oggetti geometrici tridimensionali.

4.2 – Classificazione dei Metodi Particellari I metodi particellari possono essere classificati in conformità a due

differenti criteri, principi fisici o formulazioni computazionali.

Sulla base dei principi fisici, questi possono essere suddivisi in due

differenti classi: quelli basati su modelli probabilistici, e su modelli

deterministici [7].

In base alle formulazioni computazionali si possono suddividere in due

diverse tipologie, a seconda che servono come approssimazione della

Strong Form o della Weak Form, delle equazioni differenziali parziali

(PDE) [7].

Per approssimare la Strong Form dell’equazione differenziale parziale,

questa viene solitamente discretizzata da una tecnica di collocazione

specifica, quale ad esempio lo Smoothed Particle Hydrodynamics (SPH),

il Vortex Method, il metodo generalizzato delle differenze finite, ed altri

ancora.

Sia l’SPH sia il Vortex Method, furono sviluppati inizialmente come

metodi probabilistici, mentre attualmente per alcune applicazioni

vengono impiegati come metodi deterministici.

Ciononostante, la maggior parte dei metodi particellari appartenenti a

questa categoria, sono basati su principi probabilistici, o usati come

strumenti di simulazioni probabilistiche.

Gli Strong Form Methods sono semplici da implementare, efficienti dal

punto di vista computazionale, e non necessitano di alcuna integrazione

nella definizione del sistema discreto di equazioni da risolvere. Tuttavia,


66

spesso sono instabili e poco precisi, specialmente quando si discretizza

l’elemento in un insieme di nodi distribuiti in maniera irregolare.

In questa categoria è possibile distinguere tre metodi principali: le

Dinamiche Molecolari, la simulazione diretta di Monte Carlo (DSMC), e

il Lattice Gas Automaton (LGA), da cui deriva il Lattice Boltzmann

Equation Method (LBE) [7].

Questo ultimo, è un esempio di come non tutti i metodi particellari

possono essere considerati metodi Meshfree. In quanto il LBE richiede

per la sua applicazione la definizione di una griglia, la cui funzione è

assimilabile a quella di una maglia.

Alla seconda classe di metodi particellari appartengono i Galerkin

Methods. Questi sono usati per l’approssimazione della Weak Form della

PDE; tra cui il Metodo degli Elementi Diffusi (DEM), Element Free

Galerkin Methods (EFGM), Reproducing Kernel Particle Methods

(RKPM), Metodo di partizione dell’unità (Partition of Unity Method),

Meshless local Petrov-Galerkin Method (MLPG), ed altri ancora [7].

I Weak Form Methods, sono stabili e molto precisi, tuttavia richiedono

una conoscenza sia complessiva sia locale, della Mesh realizzata, per la

risoluzione del modello, per tanto non sono considerati dei veri e propri

Meshfree Methods [7].

A questa classificazione esistono delle eccezioni, in quanto alcuni

metodi particellari possono essere usati in ambo le posizioni, quella di

forma forte, così come per la discretizzazione della forma debole. Un

esempio è il Particle-in-Cell (PIC). Per tale metodo, infatti, è possibile

distinguere due versioni, una per la Strong Form, e l’altra per la Weak

Form.

Per semplicità d’esposizione si riporta la seguente tabella grafica:

67

Latti

ce G

as A

utom

aton

(LG

A)

La S

imul

azio

ne D

iretta

di M

onte

Car

lo (D

SMC

)Le

Din

amic

he M

olec

olar

iIl

Met

odo

Gen

eral

izza

to d

elle

Diff

eren

ze F

inite

Vor

tex

Met

hod

Smoo

thed

Par

ticle

Hyd

rody

nam

ics (

SPH

)

Mes

hles

s loc

al P

etro

v-G

aler

kin

Met

hod

(MLP

G)

Parti

tion

of U

nity

Met

hod

Rep

rodu

cing

Ker

nel P

artic

le M

etho

ds (R

KPM

)El

emen

t Fre

e G

aler

kin

Met

hod

(EFG

M)

Met

odo

degl

i Ele

men

ti D

iffus

i (D

EM)

Gal

erki

n M

etho

ds

PDE

Wea

k Fo

rmPD

ESt

rong

For

mM

etod

i Pr

obab

ilist

ici

Met

odi

Det

erm

inis

tici

Crit

eri d

i cla

ssifi

cazi

one

Form

ulaz

ioni

Com

puta

zion

ali

Prin

cipi

Fis

ici

Mes

hfre

e Pa

rticl

e M

etho

ds

Tab

ella

4. 1

- C

lass

ifica

zion

e de

i Met

odi P

artic

ella

ri


68

4.3 – Applicazione dei Metodi Particellari In ambiti quali l’astrofisica, la biofisica, la biochimica e la biomedica, è

facile imbattersi in corpi che possono essere modellati non come

continui, bensì come insiemi di particelle.

In questi casi, i metodi particellari rappresentano una scelta naturale per

la simulazione numerica.

Esempi rilevanti sono le simulazioni di fenomeni quali la formazione di

sistemi stellari, movimenti di milioni d’atomi in uno stato di non

equilibrio, ed interazioni dinamiche di varie molecole.

In questo momento, i metodi particellari non sono quindi impiegati solo

come metodi di discretizzazione, per la soluzione di problemi del

continuo tridimensionale, come l’SPH, il Vortex Method, e i Galerkin

Meshfree Methods, ma anche come modelli per la descrizione dei

comportamenti fisici del continuo [7].

Gli ultimi esempi sono l’impiego del Lattice Boltzmann Equation

Method (LBE), per la risoluzione di problemi di fluido-meccanica, e il

Molecular Dynamics, per risolvere problemi inerenti alla meccanica

della frattura [7].

I metodi Meshfree Galerkin, sono largamente impiegati per la

risoluzione di problemi nell’ambito della meccanica dei corpi solidi. Tra

questi la simulazione della rottura graduale, da molti ritenuta una delle

più ostiche computazioni nell’ambito della meccanica della frattura. In

tal senso, i metodi Meshfree offrono notevoli vantaggi rispetto ai

tradizionali metodi degli elementi finiti. Per tale applicazione, infatti, si

rende necessaria un’operazione detta di Remeshing, in altre parole di

ricostruzione della maglia in seguito alla rottura. Questa operazione può

essere invece evitata utilizzando un metodo Meshfree, quale il Galerkin

o anche il Metodo di partizione dell’unità [7].


69

Tra gli altri, il Meshless Local Petrov-Galerkin (MLPG), sviluppato nel

1998 dagli studiosi Atluri e Zhu [7], risulta essere quello maggiormente

impiegato per l’analisi di travi e strutture placcate, non essendo, infatti,

necessaria una conoscenza complessiva della Mesh creata per il modello,

ma è sufficiente una conoscenza locale del modello stesso.

Un'altra applicazione in cui si preferisce usare un metodo di

discretizzazione di tipo Meshfree, rispetto ai tradizionali metodi degli

elementi finiti, è la modellazione della localizzazione degli sforzi. È

possibile dimostrare, infatti, che i metodi Meshfree hanno la capacità di

sostenere grandi distorsioni della maglia. Questo perché non c’è un

coinvolgimento diretto della maglia nella simulazione.

Invece, utilizzando il metodo FEM, avremmo un incremento della banda

di taglio lungo l’i-esimo elemento vincolato, che compromette la

simulazione del reale comportamento fisico dell’elemento stesso.

La maglia in pratica non è più in grado di seguire, per grandi

deformazioni, l’andamento reale del continuo, sotto l’azione dei carichi

esterni.

4.4 – Smoothed Particle Hydrodynamics (SPH) Tra i diversi metodi particellari di tipo Meshfree, ai fini della

modellazione numerica, si è scelto di focalizzare l’attenzione sul metodo

SPH.

Lo Smoothed Particle Hydrodynamics, è uno dei principali metodi

particellari impiegati nella meccanica computazionale.

Tale metodo fu formulato nell’1977 simultaneamente da Lucy, Gingold

& Monaghan, al fine di descrivere alcuni fenomeni Astrofisici, tra cui la

formazione, l’evoluzione, e il collasso delle stelle, delle galassie, fino

alla “Modellazione dell’intero Universo” [7].


70

Il movimento collettivo di questo insieme di particelle, è, infatti,

assimilabile a quello di un liquido o di un flusso di gas, può quindi

essere modellato utilizzando le equazioni che governano l’Idrodinamica

Classica di Newton.

L’SPH, fu ampiamente adottato come un efficiente tecnica

computazionale, per la risoluzione di problemi di meccanica applicata,

tra cui la fluido-dinamica, e l’ idrodinamica classica.

Al contrario di un convenzionale metodo di discretizzazione, che

suddivide un sistema continuo in un numero discreto di elementi finiti,

l’SPH unifica un numero discreto di particelle, accomunate dalle stesse

caratteristiche, quali massa e volume, generando così un campo continuo

locale, con il fine di rappresentare il comportamento collettivo di un

sistema discreto [7].

Basti pensare alle prime applicazioni in Astrofisica e nella Meccanica

Quantistica, dove il sistema fisico reale a già per se discreto. In tal senso,

l’SPH può essere considerato come un Metodo Particellare, basato su

una tecnica di omogeneizzazione.

La principale assunzione su cui si fonda il metodo SPH è che il valore di

una funzione f ad un certo punto i può essere approssimato moltiplicando

la stessa funzione per una funzione d’omogeneizzazione W (Kernel

function), di ampiezza h, ed integrando sul dominio computazionale nel

modo seguente [9]:

dove fi = f(xi), mentre Wi (x) = W (x−xi, h) è la funzione

d’omogeneizzazione. Questa integrazione avviene su un dominio

costituito da un insieme di particelle, il cui volume elementare è definito

( )dxxWxff ii ∫= )(


71

dal rapporto tra la massa e la propria densità. E’ quindi necessario

definire una Kernel Function W(x, h), individuando le coordinate del

nucleo W, rispetto ad un generico sistema di riferimento cartesiano, e il

periodo d’omogeneizzazione h (Smoothing Lenght). Di solito, nel

tradizionale metodo SPH, si presume che la funzione

d’omogeneizzazione soddisfi la condizione di normalizzazione nel

proprio dominio di pertinenza:

In ogni modo, è facile dimostrare che questa condizione di

normalizzazione per la funzione d’omogeneizzazione è valida solo per i

punti interni al dominio, ma non è soddisfatta per i punti di confine.

L’equazione precedente viene discretizzata come segue:

Dove Wij = W(xj - xi, h), mj e ρj sono rispettivamente la massa e la

densità di massa della j-esima particella, il rapporto tra questi due è il

volume elementare associato alla particella stessa, N è il numero totale di

particelle in cui si è discretizzato il continuo.

Tale funzione sarà poi impiegata per approssimare la strong form di un

equazione differenziale parziale (PDE), attraverso una complessa

integrazione.

Con il metodo tradizionale SPH tale approssimazione sarà:

1),( =−∫ dxhxxW i

ijj

N

j j

ji Wf

mf ∑

=⎟⎟⎠

⎞⎜⎜⎝

⎛=

1 ρ


72

considerando la proprietà della funzione nucleo:

L’equazione può essere riscritta come:

Le approssimazioni particellari delle derivate parziali possono quindi

essere scritte come:

Il valore numerico di una funzione f e le sue derivate parziali, possono

essere quindi ottenuti con una funzione kernel su un insieme di particelle

piuttosto che su una maglia. Questa è l'essenza del metodo tradizionale

SPH.

4.4.1 – Metodi d’interpolazione SPH Essendo l’SPH un metodo Meshfree, l’interpolazione è costruita su un

insieme di punti nodali, indipendenti tra loro e non vincolati da nessuna

maglia o griglia.

Questa operazione d’interpolazione è basata sul seguente concetto [7]:

( ) ( ) ( ) '''

xRdxAxxxA

nΩ−= ∫ δ Ω∈∀x

dxWxff iixxi ∫= )(

0=∫ dxW iix

[ ] dxWxfxff ixixi ∫ −= )()(

[ ] ijiij

N

j j

ji Wxfxf

mf ∇−⎟

⎟⎠

⎞⎜⎜⎝

⎛= ∑

=

)()(1 ρ


73

dove δ(x) è la funzione Dirac Delta, definita a sua volte come il limite

della seguente funzione δε(x) [7].

con:

La funzione Dirac Delta, non può però essere impiegata né per un

processo d’interpolazione, né tanto meno per un processo di

collocazione, questo perché è una funzione di tipo generalizzata, in altre

parole affetta da alcune patologie. Inoltre, manca d’alcune proprietà

fondamentali, non essendo né continua né derivabile nel proprio dominio

di definizione.

Per superare tali inconvenienti, s’introduce uno Smooth Kernel W(x, h)

(Nucleo omogeneo), al fine di utilizzare la funzione Dirac Delta per

processi d’interpolazione e collocazione.

L’SPH è quindi un operatore di media (Average Operator), definito nel

seguente modo [7]:

Ed approssimato come:

)()( lim0

xx εε

δδ→

=

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧

><<−

−<=

→ 2/2/2/

2/

;0;/1

;0)( lim

0 εεε

εεδ

εε

xx

xx

( ) ( ) ( ) ''' ,

xRdxAhxxWxA

dΩ−>=< ∫

( ) ( ) II

N

iI VxAhxxW Δ−≈ ∑

= 1,


74

Dove la notazione < … > indica l’operatore di media, mentre h, come

abbiamo già detto, è la Smoothing Lenght, in altre parole l’ampiezza o

periodo d’omogeneizzazione, infine ΔVI è il volume elementare dell’i-

esima particella.

Assumendo che il continuo ha una densità di massa distribuita

uniformemente, possiamo dividerlo in un numero discreto d’elementi

volumetrici, cui associa una massa mI. Ogni elemento volumetrico con

massa può essere attribuito ad una particolare particella.

Per cui avremo [7]:

Dove con ρI indichiamo la densità di massa dell’i-esima particella.

E’ possibile notare, che l’operatore di media <A(x)>, è una funzione

della variabile spaziale x, infatti, <A(x)> può essere visto come una

rappresentazione non locale di A(x).

La funzione Kernel deve soddisfare le seguenti proprietà [7]:

La seconda proprietà è la cosiddetta condizione di normalizzazione,

insieme alla terza sono proprietà comuni alla funzione Dirac Delta.

La terza proprietà assicura la convergenza del metodo, l'ultima deriva

dal requisito che la funzione deve essere derivabile più di una volta.

( )( ) 1;),(

0;),(

1),(

0),(

0

3

≥∈

→→

=Ω

≥

∫

kRChxW

hxhxW

dhxW

hxW

dk

xR

δ

I

II

mVρ

=Δ


75

Il vantaggio di usare un nucleo analitico, risiede nella possibilità di

valutare una Kernel Function in un qualsiasi punto dello spazio, senza

conoscere la distribuzione locale delle particelle.

La rappresentazione a nucleo non solo consente di compiere una

discretizzazione uniforme di un equazione differenziale parziale, ma

fornisce anche uno schema interpolante per un insieme di particelle in

movimento.

Con l’SPH è quindi possibile trattare un corpo solido al pari di un fluido

viscoso, cioè come un insieme di particelle, di assegnata massa e

volume, indipendenti le une dalle altre.

Sfruttando tali caratteristiche, questo metodo può essere impiegato per la

risoluzione di problemi di natura meccanica.

Ad esempio, Libersky applicando questo metodo ai solidi meccanici, è

riuscito a simulare in 3D alcuni fenomeni, tra cui l’esplosione di una

bomba, problemi di frammentazione, impatto e penetrazione di solidi.

La simulazione d’impatti e penetrazioni con l’ausilio dell’SPH, sono

state eseguite anche da Johnson, il quale impiega tale metodo per la

modellazione di problemi inelastici, di danno e di grandi deformazioni.

E’ possibile distinguere diverse funzioni Kernel, che sono

comunemente impiegate nelle computazioni. Ad esempio, la Kernel

Function può essere una funzione positiva di tipo Gaussiana, così

definita [7]:

( )31

,exp1),( 2

2

2/2

≤≤

⎥⎦

⎤⎢⎣

⎡−=

n

hx

hhxW n

π


76

Dove n è la dimensione dello spazio in cui la funzione è definita, x è la

coordinata rispetto ad un generico sistema di riferimento cartesiano

dell’i-esima particella, ed h è la Smoothing Lenght.

Altri esempi di Kernel Function sono la funzione di Spline cubica e

quartica, di seguito definite [7]:

Cubic Spline Function

Dove q = x/h; mentre C è il fattore di normalizzazione, che assume

valori diversi in funzione della dimensione dello spazio in cui è definita

la funzione:

Quartic Spline Function

Dove la costante C è determinata dalla condizione di normalizzazione,

che la funzione deve soddisfare.

( )

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

≤≤−

≤≤+−

=

casialtriin

qq

qqq

hChxW n

;0

21;241

10;4/32/31

),( 3

32

⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧

===

=3;/12;7/101;3/2

ddd

Cπ

π

⎭⎬⎫

⎩⎨⎧

≥<<−+−

=1;0

10;3861),(

432

qqqqq

hChxW n


77

Vediamo come si sceglie l’ampiezza d’omogeneizzazione h (Smoothing

Length).

Fissata l’ampiezza d’omogeneizzazione h, questa definisce una regione

che limita la maggior parte della funzione Kernel. La precisione

dell’interpolazione dipende dal numero discreto di particelle contenute

entro questa regione. La scelta della Smoothing Length è di

fondamentale importanza nell’implementazione al computer del metodo

SPH. Tale scelta, è solitamente operata in funzione della precisione

d’interpolazione che si vuole conseguire. Per una data precisione del

calcolo, infatti, è fissato il numero di particelle contenute nella regione

definita da h. In realtà, non è possibile determinare univocamente

l’ampiezza d’omogeneizzazione h, in quanto in una simulazione

dinamica, la densità di particelle contenute in detta regione varia al

variare del tempo. Per tanto, affinché si mantenga costante la precisione

numerica e l’efficienza della computazione, la Smoothing Lenght

dovrebbe variare anch’esso nel tempo.

Per fare questo, è necessario scegliere una funzione Kernel di tipo

Gaussiana, e assumere che ogni particella abbia approssimativamente la

stessa massa Δm [7].

Si deve quindi assumere che:

Una volta scelta la massa da attribuire ad ogni singola particella Δm, in

virtù dell’interpolazione SPH, potremo scrivere quanto segue:

( )t

hxx

I J

ij cosexp 2

2

==⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧ −

−∑ ∑Λ∈ Λ∈

γ

( )( ) ( ) ( )

ddI J

ijdd h

mmh

xxh 2/

22

2

2

2/ exp1πγ

πρ Δ

=Δ⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧ −

−>=< ∑ ∑Λ∈ Λ∈


78

Dove, con i pedici i e j s’indicano le generiche particelle appartenenti al

domino Λ, d è la dimensione del dominio, e xi e xj le posizioni delle

generiche particelle rispetto all’origine del sistema cartesiano.

Da quest’ultima relazione possiamo quindi trarre la seguente

conclusione [7]:

Se tutte le particelle contenute nel dominio Λ hanno quindi la stessa

quantità di massa Δm, e la stessa ampiezza d’omogeneizzazione h,

questa ultima dovrà quindi essere proporzionale al reciproco

dell’inverso della densità media delle particelle.

Questo perché, la Smoothing Lenght è proporzionale al volume

elementare ΔVi dell’i-esima particella, e quindi alla distanza di

separazione tra le particelle stesse.

4.4.2 – Teoria d’Approssimazione SPH Essendo l’SPH un metodo di computazione numerico, in quanto tale è

fondato su un insieme di principi e regole di approssimazione. Queste,

sono le cosiddette regole di Monaghan, e rappresentano le linee guida

per discretizzare un equazione differenziale [7].

Si possono distinguere quattro diverse regole d’approssimazione,

proposte dallo stesso Monaghan. La prima di queste, è quella secondo

cui è sempre meglio scegliere come funzione Kernel una funzione di tipo

Gaussiana. Tra le varie funzioni Kernel, quella Gaussiana è la migliore

nell’analizzare la coerenza fisica di un modello realizzato in SPH.

dh /1/1 ><∝ ρ


79

La seconda regola d’approssimazione, è che l’insieme del prodotto di

due funzioni, può essere approssimato dal prodotto degli insiemi

individuali:

Dove, A e B sono due funzioni, mentre <…> è l’operatore di media

(Average Operator).

La terza regola d’approssimazione è la seguente:

ed è valida per ogni campo scalare A.

La quarta regola di approssimazione, non è altri che una proprietà della

terza regola sopra citata, ovvero:

Quando la funzione Kernel è limitata, il secondo termine del secondo

membro della precedente equazione si annulla. Infatti, sfruttando la terza

regola d’approssimazione (3), possiamo scrivere:

essendo:

In virtù della seconda proprietà della funzione Kernel, sostituendo si

avrà:

><⋅>>=<⋅< BABA)2(

><∇>=∇< AA)3(

>∇<−><∇>=∇< WAAA)4(

><∇>≈∇< WW

jj

j VhxxWhxW Δ−= ∑Λ∈

),(),(

0)1( =∇>≈<∇ W


80

4.4.3 – Approssimazione di Derivate SPH L’SPH viene impiegato per approssimare la Strong Form di un

equazione differenzia parziale. Al contrario del metodo delle differenze

finite, che per discretizzare un equazione differenziale parziale, impiega

una griglia, l’SPH, essendo un metodo particellare di tipo Meshfree,

discretizza la PDE sulla base delle sopra citate regole di

approssimazione.

Si possono distinguere tre diverse formule gradienti, per

l’approssimazione di un equazione differenziale parziale [7].

Formula gradiente I

Una formula d’approssimazione diretta, impiegando un metodo

d’interpolazione SPH, può essere la seguente:

Un’approssimazione così diretta è solitamente poco accurata, e spesso

distrugge la proprietà di conservazione del sistema continuo associato.

Quando l’approssimazione è combinata con un termine addizionale, che

contiene un’espressione nulla del tipo <W(x,h)> = 0, questo può

produrre un migliore risultato.

Dalla regola d’approssimazione (3) di Monaghan, si ricava la formula

approssimata del gradiente per un campo scalare A(xi) come:

( ) )(,)( IVAhxxWxA III

I Δ−∇=∇ ∑Λ∈

( ) jJIj

IJI VWAAxA Δ∇−=∇ ∑Λ∈

,)(


81

La stessa approssimazione è valida anche per un campo vettoriale A:

Formula gradiente II

Nella formulazione in SPH, la densità è una grandezza di fondamentale

importanza. La densità di massa, infatti, stabilisce una correlazione che

sussiste tra la massa e il volume di una singola particella.

Per conseguire una precisione maggiore nella formulazione del gradiente

d’approssimazione della PDE, è preferibile scrivere il gradiente del

campo scalare e del campo vettoriale come segue:

Nella prima espressione A è un operatore scalare, nella seconda un

operatore vettoriale.

Quindi, il gradiente dell’interpolazione in SPH, può essere approssimato

nel seguente modo:

( ) jJIj

IJI VWAAxA Δ∇⋅−=⋅∇ ∑Λ∈

,)(

I

II V

mΔ

=ρ

[ ] ρρρ /)( ∇−∇=∇ AAA

[ ] ρρρ /)( ∇⋅−⋅∇=⋅∇ AAA

)()( IImWAAA jJ

IJIIJII ∑Λ∈

∇−=>∇<ρ


82

Rispetto alla prima formulazione non si fa nient’altro che esprimere il

volume dell’i-esima particella in funzione della sua densità di massa, e

della massa stessa.

Una simile formulazione può essere impiegata per approssimare il

gradiente del campo vettoriale A, facendo distinzione tra i prodotti

scalari, tensoriali, e vettoriali tra l’operatore gradiente e il campo

vettoriale:

Formula gradiente III

Nonostante le precedenti formule d’approssimazione del gradiente siano

sufficientemente accettabili, queste sono in ogni modo affette da una

deficienza comune. Queste formulazioni, infatti, non sono simmetriche

rispetto agli indici I e J dell’operatore A. Per ottenere

un’approssimazione simmetrica del gradiente, è necessario introdurre

una nuova regola d’approssimazione. Introduciamo quindi la successiva

identità:

jJ

IJIIJII mWAAscalare ∑Λ∈

∇⋅−=>⋅∇<→ ρ

jJ

IJIJIII mAWAtensoriale ∑Λ∈

⊗∇−=>⊗∇<→ ρ

jJ

IJIIJII mWAAvettoriale ∑Λ∈

∇×=>×∇<→ ρ

⎟⎟⎠

⎞⎜⎜⎝

⎛∇+⎟⎟

⎠

⎞⎜⎜⎝

⎛∇=

∇−−− 121

11σσσσ ρρρρρAAA


83

Dove σ è un numero intero. In base a questa identità è possibile scrivere

la seguente regola di approssimazione SPH:

A seconda che il parametro σ assumerà valore pari a 1 o 2, avremo le

seguenti formule d’approssimazione del gradiente, sia per il campo

scalare sia per quello vettoriale.

Quando σ = 1, per il campo scalare A(x) avremo:

Mentre, per il prodotto scalare tra l’operatore gradiente e il campo

vettoriale A, si avrà:

Quando invece σ = 2, per il campo scalare A(x) avremo:

Mentre, per il prodotto scalare tra l’operatore gradiente e il campo

vettoriale A, si avrà:

jIJIj JI

I

JI

J

I

mWAAA∇⎟⎟

⎠

⎞⎜⎜⎝

⎛+=⎟⎟

⎠

⎞⎜⎜⎝

⎛ ∇ ∑Λ∈

−− σσσσ ρρρρρ 22

)()( IIIaVWAAA JJ

IJIJII Δ∇+=>∇< ∑Λ∈

JJ

IJIJII VWAAA Δ∇⋅+=>⋅∇< ∑Λ∈

)(

)()( 22 IIIbmWAAAJ

JIJI

I

I

j

jI ∑

Λ∈

∇+=>∇

<ρρρ

JJ

IJII

I

j

jI mWAAA ∑

Λ∈

∇⋅+=>⋅∇

< )( 22 ρρρ


84

L’SPH non è solo uno schema di interpolazione, ma fornisce anche un

insieme di regole di approssimazione, che consentono di costruire un

sistema dinamico discreto del continuo.

Prima di capire come usare queste regole d’approssimazione, è

necessario descrivere le leggi di conservazione nella meccanica del

continuo.

Queste sono [7]:

• Equazione di Continuità

Dove ρ è la densità di massa, e v il campo vettoriale delle velocità.

• Equazione del Moto

Dove σ è la tensione di Cauchy.

• Equazione dell’Energia

Dove E è l’energia specifica interna, somma dell’energia Potenziale e di

quella Cinetica.

Esistono diversi modi per scrivere un’equazione differenziale in SPH,

utilizzando le regole d’approssimazione, da impiegare in seguito nella

pratica computazionale, tra cui l’equazione di Continuità, del Moto e

0=⋅∇+ vdtd ρρ

σρ

⋅∇−=1

dtdv

vdtdE

∇−= :1 σρ


85

dell’Energia, che serviranno per descrivere l’evoluzione nel tempo del

contino discretizzato in SPH [7].

Equazione di continuità in SPH

L’equazione di continuità SPH impiegata nella computazione, è ricavata

da differenti regole d’approssimazione.

Dall’equazione di continuità sopra indicata, introducendo l’operatore

medio (<….>) si ottiene quanto segue:

Utilizzando nell’ordine le regole d’approssimazione 2, 3 e 4, in

precedenza elencate, il secondo membro dell’equazione di continuità può

essere scritto come:

Sapendo che:

Sostituendole nell’equazione di continuità, scritta con le regole di

approssimazione, si avrà:

>⋅∇<−>=< vdtd ρρ

4

3

2

regolaWvv

regolav

regolav

>∇<⋅+><∇⋅><−≈

><∇⋅><−≈

>⋅∇><<−≈

ρρ

ρ

ρ

IJI

N

J j

j

IJIj

N

j j

j

Wm

W

Wvm

v

∇>=∇<

∇>=<∇⋅

∑

∑

=

=

1

1

ρ

ρ


86

Quella appena scritta altri non è che l’equazione differenziale di

continuità discretizzata come un equazione algebrica.

Si è visto come, utilizzando le regole d’approssimazione dell’SPH, sia

possibile discretizzare un equazione differenziale, in una più semplice

equazione algebrica.

Equazione del Moto in SPH

Per descrivere l’evoluzione di un sistema utilizzando l’SPH, è necessario

riscrivere le equazioni del moto in una forma discretizzata.

In funzione del tipo di regola d’approssimazione del gradiente, impiegata

per la discretizzazione dell’equazione differenziale del moto, si possono

distinguere tre diverse equazioni del moto in SPH, usate nella pratica.

Infatti, dalla seguente equazione del moto:

Dove, ρ è la densità di massa dell’i-esima particella, v è il campo

vettoriale delle velocità, e σ la tensione di Cauchy, impiegando la IIIa

regola d’approssimazione del gradiente (σ = 1), il secondo membro di

quest’equazione può essere cosi discretizzato:

( ) IJIJI

N

j j

jI

I Wvvm

dtd

∇⋅−= ∑= 1 ρ

ρρ

)0(>⋅∇<−>=< σρdtdv

( ) JIJI

N

JJI VW Δ∇⋅+−>≈⋅∇<− ∑

= 1σσσ


87

Sapendo che la densità di massa si esprime come:

Sostituendo si ottiene l’espressione discretizzata SPH dell’equazione del

moto:

Se invece di usare la regola d’approssimazione del gradiente IIIa,

s’impiega la IIIb, per la quale σ = 2, dividendo per la densità di massa,

l’equazione del moto può anche essere scritta come:

Impiegando la IIIb riferita al campo vettoriale σ, il secondo membro di

quest’equazione può essere approssimato come:

Moltiplicando ambo i membri dell’equazione per la massa mi, osi ottiene

la seguente espressione dell’equazione del moto:

Una terza espressione dell’equazione del moto, si può ottenere

utilizzando la regola I d’approssimazione del gradiente. Infatti, sapendo

I

II V

mΔ

=ρ

( ) )1(1

IJI

N

JJIJI

II WVV

dtdvm ∇⋅+ΔΔ−= ∑

=

σσ

>⋅∇

<−>=<ρ

σdtdv

JIJI

N

J J

J

I

I mW∇⋅⎟⎟⎠

⎞⎜⎜⎝

⎛−−>≈

⋅∇<− ∑

= 122 ρ

σρσ

ρσ

)2(1

22 IJI

N

J J

J

I

IJI

II Wmm

dtdvm ∇⋅⎟⎟

⎠

⎞⎜⎜⎝

⎛−−= ∑

= ρσ

ρσ


88

che il secondo membro dell’equazione (0), in virtù della regola

d’approssimazione 4, si può scrivere come:

Essendo la funzione Kernel limitata, per quanto già dimostrato, il

secondo termine al secondo membro della precedente equazione, sarà

nullo.

Di conseguenza:

Il cui secondo membro può essere discretizzato secondo la regola I

d’approssimazione del gradiente, valida per il campo vettoriale σ:

Sostituendo l’approssimazione, e ricordando l’espressione della densità

di massa ρ, l’equazione del moto diventa:

In queste tre diverse espressioni dell’equazione del moto, il primo

membro, che è chiaramente sempre lo stesso, fisicamente esprime la

variazione nel tempo della quantità di moto. Infatti, la quantità di moto

dell’i-esima particella, è data dal prodotto della massa mI, che è costante

e quindi indipendente dal tempo, e della velocità vI, variabile in funzione

del tempo.

)( >∇><<−><∇⋅−>=⋅∇<− Wσσσ

><−∇⋅>=⋅∇<− σσ

IJIIJ

N

JJ WV ∇⋅−Δ>≈<∇⋅− ∑

=

)(1

σσσ

( ) )3(1

IJIIJ

N

JJI

II WVV

dtdvm ∇⋅−ΔΔ−= ∑

=

σσ


89

Equazione dell’energia in SPH

Il secondo membro dell’equazione dell’energia, utilizzando la seconda

regola d’approssimazione di Monaghan, può essere riscritto come:

Questo termine, può essere ulteriormente discretizzato, impiegando la II

regola d’approssimazione del gradiente per il campo vettoriale σ,

l’equazione dell’energia può quindi essere espressa come:

Le equazioni discretizzate con il metodo SPH, non sempre preservano le

caratteristiche fisiche del corrispondente sistema continuo, come

l’energia e la quantità di moto. Questo problema incide notevolmente

sulla precisione e la qualità delle soluzioni numeriche.

Dal punto di vista computazionale, un buono schema d’interpolazione,

dovrebbe essere almeno capace di rappresentare correttamente il moto di

un corpo rigido. Purtroppo, il metodo d’interpolazione SPH non è in

grado di rappresentare propriamente il moto di un corpo rigido. Ciò

induce ad una serie d’ostacoli nella computazione numerica. Sebbene, un

campo d’interpolazione non riesca a rappresentare il movimento

completo del corpo rigido, essendo questo composto di una traslazione e

una rotazione rigida, la sola traslazione può essere preservata in una

computazione discreta.

E’ possibile dimostrare che le equazioni di continuità e dell’energia,

prima approssimate con il metodo SPH, sono degli invarianti Galileiani.

>∇<><><

−>≈∇<−>≈∇<− vvv ρσρ

ρσρ

σρ

:1)(:1:122

)1(:1

2 IJIIJ

N

JJ

I

II Wvmdt

dE∇⊗⎟⎟

⎠

⎞⎜⎜⎝

⎛= ∑

=ρσ


90

Infatti, se si scrivono le equazioni di continuità e dell’energia rispetto ad

un sistema di riferimento cartesiano assoluto, le stesse equazioni

permangono inalterate in un sistema di riferimento cartesiano in moto

relativo rispetto a quello assoluto.

La posizione e la velocità della generica particella, in cui si è

discretizzato il continuo, rispetto al riferimento relativo, possono essere

così espresse:

Dove x’ e v’, sono rispettivamente la posizione e la velocità della

particella rispetto al sistema di riferimento relativo, mentre x e v sono

definite rispetto a quello assoluto, a e t sono invece l’accelerazione e

l’istante di tempo.

Ovviamente, la velocità si ottiene in funzione della posizione dividendo

questa ultima per il tempo.

L’equazione di continuità e dell’energia in SPH, scritte nel nuovo

riferimento diventano quindi:

Dato che, la differenza tra le velocità in seguito al cambio di riferimento

resta in ogni modo invariato:

le equazioni stesse saranno invarianti rispetto al riferimento cartesiano.

avvatxx

+=+=

''

( ) IJI

N

j j

jI

I Wvvm

dtd

JI∇⋅−= ∑

=

''

1 ρρρ

( ) IJIJ

N

JJ

II WVvmdt

dEI

I

∇⋅−= ∑=

''

12ρ

σ

( ) ( )JIJI vvvv −=− ''


91

4.4.4 – I limiti del Metodo SPH I metodi SPH presentano notevoli vantaggi nelle computazioni

numeriche, tra cui la facilità di trattare grandi deformazioni della Mesh, o

la semplicità d’implementazione al calcolatore.

Ciononostante, lo stesso metodo presenta una serie di limiti ed

inconvenienti di seguito elencati [8].

• Instabilità di tensione;

• Insufficienza di consistenza d’interpolazione;

• Difficoltà nel far rispettare le condizioni di vincolo.

Attraverso gli anni, sono state sviluppate varie tecniche per il

miglioramento del metodo, e l’eliminazione delle varie patologie che si

manifestano nelle computazioni numeriche, di seguito definite.

Instabilità di tensione

Si definisce in questo modo la situazione in cui, in una regione del

dominio dove sussiste uno stato tensionale, il moto delle particelle

diventa instabile. Ciò significa che, una piccola perturbazione della

posizione delle particelle, può causare un esponenziale aumento della

velocità delle stesse, che può indurle ad un moto oscillatorio.

Diversi autori hanno condotto studi con il fine di individuare le cause di

questo fenomeno, e in seguito per porvi rimedio.

Per quanto attiene alle origini di tale fenomeno si fa riferimento ai

risultati conseguiti dagli studi condotti da Sweegle [8].


92

Ipotizzando che le masse delle particelle adiacenti siano le stesse, e che

sia soddisfatta la condizione secondo la quale:

i criteri di stabilità di Sweegle vengono di seguito elaborati.

Considerando una distribuzione di particelle in cui la distanza tra due

particelle successive è pari alla Smoothing Lenght h, assumendo l’origine

del sistema di riferimento cartesiano coincidente con la generica

particella XI, e facendo riferimento ad un semplice caso

monodimensionale, la particella I+1 è posizionata alla distanza XI + h,

mentre, il punto x in cui la derivata della funzione Kernel W’(x) si

annulla, sarà compreso tra la particella XI e I+1, come di seguito

mostrato.

-W'

W'

stabile instabile

stabileinstabile

Trazione

Compressione

h

σ

XXI I+1X

Figura 4. 1 – Instabilità d trazione delle particelle.

0=I

II

dxdW


93

Assumendo per convenzione che la tensione negativa sarà di

compressione, e quella positiva di trazione, sarà possibile distinguere nel

piano x-σ le seguenti situazioni [8].

Trazione (σ>0)

• Nel tratto in cui la derivata della Funzione Kernel -W’(x) è

crescente [-W’(x)>0], all’aumentare della variabile indipendente x

le particelle si allontanano, e la tensione di trazione σ allo stesso

tempo aumenta. Viceversa, al diminuire della variabile x le

particelle si avvicinano, e la tensione σ si riduce. Il sistema in

questo caso si definisce stabile.

• Nel tratto in cui la derivata della Funzione Kernel -W’(x) è

decrescente [-W’(x)<0], all’aumentare della variabile indipendente

x le particelle si allontanano, e la tensione di trazione σ in questo

caso diminuisce. Viceversa, al diminuire della variabile x le

particelle si avvicinano, e la tensione σ aumenta. Il sistema in

questo caso si definisce instabile.

Compressione (σ<0)

• Nel tratto in cui la derivata della Funzione Kernel W’(x) è

crescente [W’(x)>0], all’aumentare della variabile indipendente x

le particelle si allontanano, e la tensione di compressione σ

diminuisce. Viceversa, al diminuire della variabile x le particelle si

avvicinano, e la tensione σ aumenta. Il sistema in questo caso si

definisce instabile.

• Nel tratto in cui la derivata della Funzione Kernel W’(x) è

decrescente [W’(x)<0], all’aumentare della variabile indipendente


94

x le particelle si allontanano, e la tensione di compressione σ allo

stesso tempo aumenta. Viceversa, al diminuire della variabile x le

particelle si avvicinano, e la tensione σ si riduce. Il sistema in

questo caso si definisce stabile.

Tra le possibili soluzioni a questo problema Morris propose di introdurre

una speciale Kernel Function, ma tale soluzione produce risultati

soddisfacenti soltanto in alcuni casi particolari [8].

Dyka propose il cosiddetto Stress Point Method, l’idea principale su cui

si fonda, risiede nell’introdurre dei punti addizionali, oltre a quelli

introdotti per la rappresentazione delle singole particelle, nel momento in

cui si vuole valutare lo stato tensionale o altri stati variabili. Mentre le

variabili cinematiche, quali lo spostamento, la velocità, e l’accelerazione

sono ancora riferite alle particelle SPH, a questi punti addizionali sono

invece associati gli stati tensionali, o altri stati variabili nel tempo.

Tali punti sono definiti Stress Point, e non sono niente altro che delle

particelle fittizie. Di seguito si riporta una rappresentazione grafica:

Figura 4. 2 - Stress Point Method in uno spazio bi-dimensionale.


95

Condizioni di vincolo

I metodi particellari in generale hanno delle difficoltà nel rafforzare le

condizioni di vincolo essenziali. Per l’SPH, sono stati sviluppati dei

metodi correttivi per risolvere questa patologia. Notevolmente impiegato

nelle computazioni pratiche, è il così detto “Ghost Particle”, proposto

da Randles e Libersky [8]. In virtù di questo metodo, si supponga che l’i-

esima particella i sia vincolata, mentre le altre particelle all’interno dello

stesso dominio N(i), possano essere suddivise in tre sottoinsiemi:

• I(i): punti interni al dominio, adiacenti la particella i;

• B(i): punti di confine, adiacenti la particella i;

• G(i): punti esterni al dominio, adiacenti la particella i.

L’insieme N(i), può quindi essere espresso come l’unione dei tre

sottoinsiemi descritti:

Dove G(i) sono appunto le particelle fantasma introdotte (Ghost

Particle) [8].

In quest’approccio, si può ricavare l’equazione di vincolo per un

generico campo scalare fi, infatti, assumendo:

)()()()( iGiBiIiN ∪∪=

ij

bcj

ij

VV

iGjff

iBjff

Δ=Δ

∈∀=

∈∀=

)(

)(


96

la formula di correzione per il generico campo scalare f, è definita come

segue:

Dove:

• fbc è il prescritto valore per x = xi;

• ΔVj è il volume elementare della j-esima particella;

• Wij è la funzione Kernel.

Il principale vantaggio nell’utilizzare questa formula correttiva, risiede

nel fatto che questa dipende esclusivamente dalle particelle interne, e

non da quelle esterne (Ghost Particle).

( )

⎟⎟⎠

⎞⎜⎜⎝

⎛Δ−

Δ−+=

∑

∑

∈

∈

)(

)(

1iBj

ijj

iIjijjbcj

bci

WV

WVffff

CCAAPPIITTOOLLOO 55 –– MMOODDEELLLLAAZZIIOONNEE NNUUMMEERRIICCAA CCOONN LLSS--DDYYNNAA VVEERRSSIIOONN 997700

97

Capitolo 5 Modellazione Numerica con LS-Dyna Version 970.

Modellare numericamente una prova sperimentale significa riprodurre,

grazie all’ausilio di un codice di calcolo, lo stesso elemento oggetto della

sperimentazione. Le evidenze sperimentali serviranno per poter

verificare l’attendibilità del modello creato al calcolatore, oppure i

risultati dei test potranno essere impiegati per calibrare il modello stesso.

Tutto ciò è finalizzato alla costruzione di un modello che consente di

eseguire sperimentazioni numeriche, senza ricorrere necessariamente alle

prove sperimentali da realizzare in laboratorio, con i costi e le

problematiche che esse comportano.

Individuate quali sono le principali finalità della modellazione numerica,

bisogna poi in seguito scegliere tra i codici di calcolo strutturali a

disposizione, quello da impiegare per la simulazione stessa.

Ai fini della modellazione, si è scelto di impiegare il programma ad

elementi finiti d’uso generale LS-Dyna Version 970 [11], che consente

la simulazione di problemi complessi, ed è largamente accettato, come

software d’analisi, per le più avanzate applicazioni dell'ingegneria.

5.1 - Software LS-Dyna Version 970 LS-Dyna è un programma di calcolo ad elementi finiti, che permette di

analizzare dei comportamenti fisici altamente non lineari con dei metodi

numerici detti espliciti o impliciti. Frequenti applicazioni si hanno nel

caso di fenomeni in cui grandi deformazioni avvengono in un breve

intervallo di tempo, come ad esempio in un Crash test, o in

un’esplosione.


98

LS-Dyna conta numerosi utilizzatori nell’industria automobilistica,

aeronautica, aerospaziale e nella formatura di metalli. Il calcolo

accoppiato fluido-struttura permette di simulare l’interazione tra una

struttura e un fluido come nell’aquaplanning o nella caduta di un

recipiente contenente un fluido.

Le origini del LS-Dyna risalgono alla metà degli anni settanta, quando fu

sviluppata la sua prima versione, nota come DYNA3D, per opera del

Lawrence Livermore National Laboratory.

Il DYNA3D nasce come software d’analisi, sulla base della teoria degli

elementi finiti, per la risoluzione di problemi strutturali e di meccanica

del continuo. Grazie alla sua rapidità nell’integrazione delle equazioni

del moto, è specialmente impiegato per la risoluzione di problemi

dinamici.

Figura 5. 1 - Modellazione di una Turbina con LS-Dyna.

Essa possiede una libreria di materiali (Material Model), da poter

impiegare nelle modellazioni, tra i quali sono inclusi materiali isotropici-

elastici, ortotropici-elastici ed elasto-plastici. Consente la simulazione

d’elementi solidi, shell, beam, spring, e damper. Inoltre si possono


99

generare contatti tra superfici, in modo tale da considerare le interazioni

tra corpi diversi. Col passare degli anni dal 1976 fino al 1988, il

DYNA3D è stato di volta in volta arricchito con nuove caratteristiche e

possibilità, riguardanti materiali ed elementi da modellare. Nelle versioni

che si sono succedute sono stati apportati miglioramenti ai tempi

d’integrazione, per la risoluzione delle equazioni del moto, riducendo

notevolmente i tempi d’analisi, e consentendo discretizzazione sempre

più fitte degli elementi.

Verso la fine del 1988, dal DYNA3D nasce il più recente codice di

calcolo LS-Dyna. Questa evoluzione del DYNA3D è dettata dalla

crescente richiesta d’applicazioni automobilistiche. Gli sforzi per lo

sviluppo del codice di calcolo furono intensificati, proprio per far fronte

alla maggiore richiesta di risoluzione di problemi automobilistici. Agli

inizi del 1989 la scrittura del codice passò dalla Lawrence Livermore

National Laboratory (LLNL), alla Livermore Software Technology

Corporation, California. Per opera della LSTC sono state migliorate

diverse capacità del software. Da questo momento, per i successivi 18

anni, passando attraverso le diverse versioni del programma, dalla

Version 940 fino alla recentissima Version 971, innumerevoli sono state

le migliorie adottate. Tra queste si riportano di seguito solo alcune delle

più innovative ed interessanti, che caratterizzano la versione più recente,

impiegata per la modellazione numerica. Le sue capacità attuali

consentono la soluzione di Nonlinear dynamics, Rigid body dynamics,

Quasi-static simulations, Normal modes, Linear static, Thermal analysis,

Fluid analysis, Eulerian capabilities, ALE (Arbitrary Lagrangian-

Eulerian), Fluid-structure interactions, FEM-rigid multi-body dynamics

coupling (MADYMO, CAL3D), Underwater shock, Failure analysis,

Crack propagation, Real-time acoustics, Design optimization, Implicit


100

springback, Multi-physics coupling, Structural-thermal coupling,

Adaptive Remeshing, Smooth particle hydrodynamics, Element-free

meshless method [11].

LS-Dyna utilizza per la soluzione delle equazioni del moto sia il metodo

esplicito sia quello implicito. La disponibilità in un unico prodotto

software di metodi sia espliciti sia impliciti è importante, perchè un

singolo metodo di soluzione non è convenientemente applicabile in

qualsiasi situazione. Con un solo prodotto, LS-Dyna, si è in grado di

portare a soluzione le problematiche più diverse, tra cui:

• Solid Mechanics.

• Nonlinear elements for large deformations, Reduced and fully

integrated, linear elements for eigenvalues, superelements, and

linear structural analyses.

• Dynamics: Explicit methods for short duration transient problems,

implicit methods for static and long duration problems,

instantaneously switch between methods, Fluid Mechanic.

• Flow regim: Incompressible flow, Compressible flow, Acoustics.

Il metodo esplicito è in pratica l'unico applicabile in diversi casi, come

ad esempio le simulazioni d’impatto (crash e drop test), il metal forming

(stampaggi a caldo e a freddo, laminazioni, imbottiture), le simulazioni

d’eventi prodotti da esplosioni, sia in aria sia in acqua, oppure gli effetti

prodotti dai terremoti sulle strutture. Può essere impiegato per le

simulazioni su materiali non metallici quali compositi o schiume

polimeriche. Le schiume sono i principali componenti, ad esempio, dei

sedili d’auto e pertanto costituiscono un elemento di primaria importanza


101

per la sicurezza degli occupanti. In ultimo, ma non meno importante, con

il metodo esplicito si è in grado di simulare tutte le fasi d’esercizio di un

airbag, dall'esplosione della carica, al dispiegamento durante il

gonfiaggio, all'impatto del corpo umano sul cuscino gonfiato, allo

sgonfiaggio controllato.

Tra le potenzialità del software in precedenza elencate, quelle che

maggiormente ci hanno indotto a sceglierlo per la simulazione numerica,

sono senza dubbio la possibilità di implementare due metodi di

discretizzazione, dettagliatamente descritti nel capitolo precedente, il

Finite Element Method (FEM), e lo Smoothed Particle Hydrodynamics

(SPH), e in secondo luogo la capacità del codice di calcolo di

considerare grandi deformazioni che si verificano in un intervallo di

tempo di pochi millisecondi, così come accade nella simulazione di

un’esplosione.

5.2 – Modellazione della Prova Sperimentale Con l’ausilio del codice di calcolo LS-Dyna Version 970, si sono

modellate le quattro prove sperimentali in precedenza descritte.

Si è scelto di impiegare le seguenti unità di misura [12].

Lunghezze Tempo Pressione Massa Temperatura Densità [m] [s] [Pa] [Kg] [K] Kg/m3]

Tabella 5. 1 – Unità di Misure adottate dal LS-Dyna Version 970.

Tale modellazione è stata condotta utilizzando due metodi di

discretizzazione differenti, il FEM e l’SPH. Si descrivono in seguito

quali sono le caratteristiche dello stesso modello realizzato con i due

metodi differenti.


102

5.2.1 – Modello FEM L’elemento strutturale piastra viene in questo caso realizzato mediante

un Box_Solid [12]. Ovviamente la scelta del numero d’elementi che

costituiscono il Box, e quindi la piastra, è effettuata in funzione del grado

d’accuratezza che si vuole conseguire nel calcolo, associato al numero

d’informazioni che da lui si vogliono ottenere, e tenendo presente che un

numero eccessivo d’elementi solidi e di nodi che costituiscono la Mesh,

incrementano i tempi d’analisi del processore. Questo perchè il numero

d’equazioni che il processore dovrà risolvere aumenta all’aumentare del

numero di nodi. La piastra è quindi suddivisa in 11520 elementi solidi

(Solid Element) [12]. Ovviamente, trattandosi di una discretizzazione in

elementi finiti, ad ognuno di questi saranno associate le stesse

dimensioni geometriche e caratteristiche dei materiali. Le dimensioni del

singolo elemento sono (2,5x2,5x1,8)cm3, in questo modo si avrà una

maglia di 48x48 elementi distribuiti su 5 strati sovrapposti. Ogni

elemento è contraddistinto da 8 nodi, ognuno dei quali a sua volta

individuato, rispetto all’origine degli assi del sistema di riferimento

coincidente con il baricentro geometrico della piastra, dalle sue

coordinate cartesiane (x, y, z).

Ad esempio per l’elemento individuato dal numero identificativo eid

(element identification) = 1 si avrà: *ELEMENT_SOLID $# eid pid n1 n2 n3 n4 n5 n6 n7 n8 1 1 1 2 51 50 2402 2403 2452 2451

Tabella 5. 2 – Element Solid.

Dove con ni s’indica il numero identificativo del nodo che appartiene

all’elemento, coincidente con uno dei suoi vertici, mentre pid è il

numero identificativo della Part associata all’elemento.


103

La Part è un comando del codice di calcolo che consente di assegnare a

più elementi le stesse proprietà, quali ad esempio lo stesso Material, la

stessa Equation of State (EOS) e la stessa Section, richiamando

semplicemente i rispettivi numeri identificativi (mid, eosid, secid) [12]. *PART $# title Part $# pid secid mid eosid hgid grav adpopt tmid 1 1 1 1

Tabella 5. 3 – Part.

Nella Section sono definite alcune proprietà, quali la formulazione

dell’elemento, le regole d’integrazione o gli spessori nodali, secondo la

tipologia d’Element a cui la Section è associata. Nel caso della

Section_Solid si attribuisce all’elemento l’Element Formulation Optinon

(elform) [12], in questo caso specifico si sceglie di considerare costante

la tensione nel singolo elemento.

*SECTION_SOLID_TITLE Section_Solid $# secid elform aet 1 1

Tabella 5. 4 – Section Solid.

Costruito il modello geometrico della piastra, e la mappa topologica

degli elementi finiti in cui questa è stata discretizzata, in seguito si

definiscono le condizioni di vincolo e di carico.

Per quanto riguarda le condizioni di vincolo, si crea un Set di Nodi [12],

in altre parole un insieme di Nodi, a cui andrà poi in seguito associato la

condizione di vincolo. In questo caso si è scelto di vincolare lungo l’asse

y del sistema di riferimento cartesiano, i due bordi esterni della piastra. Il

Set di nodi a loro corrispondente è il seguente:


104

*SET_NODE_LIST_TITLE Vincoli $# sid da1 da2 da3 da4 37 $# nid1 nid2 nid3 nid4 nid5 nid6 nid7 nid8 1 50 99 148 197 246 295 344 393 442 491 540 589 638 687 736 785 834 883 932 981 1030 1079 1128 1177 1226 1275 1324 1373 1422 1471 1520 1569 1618 1667 1716 1765 1814 1863 1912 1961 2010 2059 2108 2157 2206 2255 2304 2353 2402 2451 2500 2549 2598 2647 2696 2745 2794 2843 2892 2941 2990 3039 3088 3137 3186 3235 3284 3333 3382 3431 3480 3529 3578 3627 3676 3725 3774 3823 3872 3921 3970 4019 4068 4117 4166 4215 4264 4313 4362 4411 4460 4509 4558 4607 4656 4705 4754 4803 4852 4901 4950 4999 5048 5097 5146 5195 5244 5293 5342 5391 5440 5489 5538 5587 5636 5685 5734 5783 5832 5881 5930 5979 6028 6077 6126 6175 6224 6273 6322 6371 6420 6469 6518 6567 6616 6665 6714 6763 6812 6861 6910 6959 7008 7057 7106 7155 7204 7253 7302 7351 7400 7449 7498 7547 7596 7645 7694 7743 7792 7841 7890 7939 7988 8037 8086 8135 8184 8233 8282 8331 8380 8429 8478 8527 8576 8625 8674 8723 8772 8821 8870 8919 8968 9017 9066 9115 9164 9213 9262 9311 9360 9409 9458 9507 9556 9605 9654 9703 9752 9801 9850 9899 9948 9997 10046 10095 10144 10193 10242 10291 10340 10389 10438 10487 10536 10585 10634 10683 10732 10781 10830 10879 10928 10977 11026 11075 11124 11173 11222 11271 11320 11369 11418 11467 11516 11565 11614 11663 11712 11761 11810 11859 11908 11957 12006 12055 12104 12153 12202 12251 12300 12349 12398 12447 12496 12545 12594 12643 12692 12741 12790 12839 12888 12937 12986 13035 13084 13133 13182 13231 13280 13329 13378 13427 13476 13525 13574 13623 13672 13721 13770 13819 13868 13917 13966 14015 14064 14113 14162 14211 14260 14309 14358 49 98 147 196 245 294 343 392 441 490 539 588 637 686 735 784 833 882 931 980 1029 1078 1127 1176 1225 1274 1323 1372 1421 1470 1519 1568 1617 1666 1715 1764 1813 1862 1911 1960 2009 2058 2107 2156 2205 2254 2303 2352 2401 2450 2499 2548 2597 2646 2695 2744 2793 2842 2891 2940 2989 3038 3087 3136 3185 3234 3283 3332 3381 3430 3479 3528 3577 3626 3675 3724 3773 3822 3871 3920 3969 4018 4067 4116 4165 4214 4263 4312 4361 4410 4459 4508 4557 4606 4655 4704 4753 4802 4851 4900 4949 4998 5047 5096 5145 5194 5243 5292 5341 5390 5439 5488 5537 5586 5635 5684 5733 5782 5831 5880 5929 5978 6027 6076 6125 6174 6223 6272 6321 6370 6419 6468 6517 6566 6615 6664 6713 6762 6811 6860 6909 6958 7007 7056 7105 7154 7203 7252 7301 7350 7399 7448 7497 7546 7595 7644 7693 7742 7791 7840 7889 7938 7987 8036 8085 8134 8183 8232 8281 8330 8379 8428 8477 8526 8575 8624 8673 8722 8771 8820 8869 8918 8967 9016 9065 9114 9163 9212 9261 9310 9359 9408 9457 9506 9555 9604 9653 9702 9751 9800 9849 9898 9947 9996 10045 10094 10143 10192 10241 10290 10339 10388 10437 10486 10535 10584 10633 10682


105

10731 10780 10829 10878 10927 10976 11025 11074 11123 11172 11221 11270 11319 11368 11417 11466 11515 11564 11613 11662 11711 11760 11809 11858 11907 11956 12005 12054 12103 12152 12201 12250 12299 12348 12397 12446 12495 12544 12593 12642 12691 12740 12789 12838 12887 12936 12985 13034 13083 13132 13181 13230 13279 13328 13377 13426 13475 13524 13573 13622 13671 13720 13769 13818 13867 13916 13965 14014 14063 14112 14161 14210 14259 14308 14357 14406

Tabella 5. 5 – Set_Node_List Vincoli.

Al Set di nodi così definito sono poi associate le condizioni di vincolo,

bloccando gli spostamenti di tali nodi nelle direzioni x y e z, realizzando

così un vincolo di tipo incastro.

*BOUNDARY_SPC_SET_ID $# cid heading 1Incasro $# nsid cid dofx dofy dofz dofrx dofry dofrz 37 0 1 1 1

Tabella 5. 6 – Boundary Condiction.

Come si può notare, nella Card che definisce le condizioni di vincolo

(Boundary Condiction) [12] una volta richiamato il nsid (Node Set

Identification), corrispondente al Set di nodi creato in precedenza, si

definiscono quali sono gli spostamenti e gli assi lungo i quali avvengono,

oppure le rotazioni e gli assi rispetto ai quali si realizzano, che devono

essere impediti.

Una volta calcolate le curve di pressione che agiscono all’estradosso

della piastra, così come descritto al Capitolo 3, queste dovranno essere

definite nel modello e in seguito assegnate agli elementi solidi.

Per definire una qualsivoglia curva, che indichi l’andamento della

pressione in funzione del tempo, o che sia in ogni modo l’andamento di

una variabile dipendente in funzione di un'altra indipendente (ese.

Temperatura e Tempo), si utilizza il comando Define_Curve[12].

Questo comando consente di assegnare in maniera tabulata, nel piano (x,


106

y), i valori da associare alle due variabili, indipendente lungo l’asse delle

ascisse e dipendente lungo l’asse delle ordinate. Inoltre, si possono anche

definire dei fattori di scala e d’Offset sia per i valori lungo l’asse delle

ascisse (sfa, offa), che per quelli definiti lungo l’asse delle ordinate (sfo,

offo). Ad ogni curva si associa un nome (Title), e un numero

identificativo lcid (Load Curve Identification) [12].

A titolo d’esempio si riporta una delle 36 curve di carico, definite per

modellare l’effetto prodotto sulla superficie della piastra, investita

dall’onda d’urto prodotta dall’esplosione: *DEFINE_CURVE_TITLE A-1 $# lcid sidr sfa sfo offa offo dattyp 1 0 1.000000 1.000000 $# a1 o1 0.000 1.6806463e+005 4.9999999e-005 1.4723016e+005 9.9999997e-005 1.2860335e+005 1.5000001e-004 1.1196810e+005 1.9999999e-004 97128.437500 2.5000001e-004 83906.765625 3.0000001e-004 72142.164063 3.4999999e-004 61688.882813 3.9999999e-004 52414.949219 4.4999999e-004 44200.937500 5.0000002e-004 36938.789063 5.5000000e-004 30530.804688 6.0000003e-004 24888.658203 6.5000000e-004 19932.556641 6.9999998e-004 15590.442383 7.5000001e-004 11797.274414 7.9999998e-004 8494.3789063 8.5000001e-004 5628.8549805 8.9999998e-004 3153.0300293 9.5000002e-004 1023.9760132 0.00100000 -796.94500732 0.00105000 -2344.4819336 0.00110000 -3649.8649902 0.00115000 -4741.1430664 0.00120000 -5643.4848633 0.00125000 -6379.4589844 0.00130000 -6969.2807617 0.00135000 -7431.0468750 0.00140000 -7780.9331055 0.00145000 -8033.3881836 0.00150000 -8201.3037109 0.00160000 -8328.1982422 0.00170000 -8239.0712891 0.00180000 -7994.9111328 0.00190000 -7643.4130859 0.00205000 -6993.4340820 0.00220000 -6273.6499023 0.00235000 -5543.2690430 0.00250000 -4839.4492188 0.00265000 -4184.0019531 0.00280000 -3588.2438965


107

0.00300000 -2893.4970703 0.00320000 -2310.4570313 0.00340000 -1829.9639893 0.00360000 -1439.5209961 0.00380000 -1125.7989502 0.00400000 0.00420000

Tabella 5. 7 – Curva di Pressione Tabulata.

Grafico 5. 1 – Andamento della Curva di Pressione, calcolata per l’area di carico

A-1 e relativa alla prima prova modellata.

Lo stesso procedimento è ripetuto per le 36 diverse curve di pressione

calcolate in precedenza, e per le quattro diverse condizioni di carico

realizzate nelle altrettante prove sperimentali eseguite.

Definite le curve di carico, bisognerà poi in seguito assegnare queste

condizioni di carico agli elementi che costituiscono la piastra stessa.

Per i modelli realizzati in FEM, si è scelto di applicare il carico in

maniera distribuita. In questo caso sono stati definiti dei Set_Segment

[12], vale a dire delle aree di carico, coincidenti con le stesse che

abbiamo in precedenza definito per il calcolo delle curve di pressioni

derivanti dall’esplosione. Per definire un Segment è sufficiente

individuare 4 nodi che lo delimitano, il singolo Segment sarà quindi una


108

semplice area quadrata. Ogni Set_Segment sarà invece costituito da 4

Segment diversi, tra loro simmetrici rispetto al centro della piastra. Si

ricorda, infatti, che la condizione di carico a cui la piastra è sottoposta è

perfettamente simmetrica rispetto al centro della stessa, essendo la carica

esplosiva posizionata ad una certa distanza dal suo centro. In definitiva

avremo 36 Set_Segment, per un totale di 144 Segment [12].

*SET_SEGMENT_TITLE A-1 $# sid da1 da2 da3 da4 1 $# n1 n2 n3 n4 a1 a2 a3 a4 14382 14186 14190 14386 12030 12034 12230 12226 14382 14378 14182 14186 12030 12226 12222 12026

Tabella 5. 8 - Set_Segment.

Ad ogni Set_Segment è poi associata la corrispondente curva di

pressione, sapendo che questa si ripete identicamente per i 4 Segment a

lei appartenenti.

A tale fine si utilizza il comando Load Set_Segment [12], che consente di

associare ad ogni ssid (Set_Segment Identification) la corrispondente

lcid (Load Curve Identification), moltiplicata per un fattore di scala sf, il

quale si assume unitario. Questo comando ci consente di applicare una

pressione distribuita su un’area quadrilatera definita da 4 nodi.

*LOAD_SEGMENT_SET $# ssid lcid sf at 1 1 1.000000 2 2 1.000000 3 3 1.000000 4 4 1.000000 5 5 1.000000 6 6 1.000000 7 7 1.000000 8 8 1.000000 9 9 1.000000 10 10 1.000000 11 11 1.000000 12 12 1.000000 13 13 1.000000


109

14 14 1.000000 15 15 1.000000 16 16 1.000000 17 17 1.000000 18 18 1.000000 19 19 1.000000 20 20 1.000000 21 21 1.000000 22 22 1.000000 23 23 1.000000 24 24 1.000000 25 25 1.000000 26 26 1.000000 27 27 1.000000 28 28 1.000000 29 29 1.000000 30 30 1.000000 31 31 1.000000 32 32 1.000000 33 33 1.000000 34 34 1.000000 35 35 1.000000 36 36 1.000000

Tabella 5. 9 – Load Set_Segment.

Figura 5. 2 - Sezione della Piastra sottoposta a Carico distribuito.

Per realizzare un modello quanto più attendibile è possibile, sarà

necessario attribuire a singoli elementi anche il proprio peso. Per fare

questo è necessario creare una Set_Part_List, in altre parole una lista di

Part alle quali associare un unico carico, il peso proprio. Si associa

quindi alla Set_Part_List l’unica Part di cui disponiamo [12].


110

*SET_PART_LIST_TITLE part list solid $# sid da1 da2 da3 da4 1 $# pid1 pid2 pid3 pid4 pid5 pid6 pid7 pid8 1

Tabella 5. 10 – Set_Part_List.

Si definisce quindi un ulteriore curva, in questo caso semplicemente

costate, essendo la forza peso sempre la stessa indipendentemente dal

tempo, si definisce un Load_Body_Z [12], in cui si moltiplica la curva di

carico in precedenza costruita per un fattore di scala sf, pari a 9,81 m/s2,

l’accelerazione di gravità, e si associa la nuova condizione di carico

all’intera Part.

*DEFINE_CURVE_TITLE Peso $# lcid sidr sfa sfo offa offo dattyp 37 0 1.000000 1.000000 $# a1 o1 0.000 1.00000000 0.00600000 1.00000000 *LOAD_BODY_PARTS $# psid 1 *LOAD_BODY_Z $# lcid sf lciddr xc yc zc 37 9.810000

Tabella 5. 11 – Carico da Peso Proprio.

Per completare il modello, non rimane altro che definire le caratteristiche

meccaniche dei materiali, calcestruzzo e acciaio, utilizzati per la

costruzione della piastra. A tal scopo si rimanda al Capitolo 6, in cui si

descriveranno in maniera dettagliata i Material Model impiegati, il modo

in cui si definisce l’Equation of State ad essi associati, ed infine le curve

di Strain-Rate, che descrivono il miglioramento del comportamento dei

materiali, per effetto dei carichi dinamici a cui sono soggetti. Queste

ultime Card che definiscono il modello sono le stesse sia per il Metodo

FEM sia per quello SPH.


111

Completato il modello in FEM per la prima condizione di carico, per

costruire gli altri tre modelli rimanenti, sarà sufficiente copiare lo stesso

modello modificando le sole curve di pressione riguardanti le altre tre

condizioni di carico.

5.2.2 – Modello SPH In questo caso per creare l’elemento piastra, a differenza di quanto fatto

per il modello FEM, non si costruisce una Mesh, bensì una mappa

topologica di soli nodi. In questo caso, infatti, la piastra non è suddivisa

in un numero discreto d’elementi finiti di forma prismatica, ma in un

numero discreto d’elementi nodali, definiti particelle. Ad ogni particella

è associata una massa ed una posizione nello spazio, individuata dalle tre

coordinate (x, y, z), sempre rispetto all’origine degli assi del riferimento

cartesiano considerato dal programma. Ad ogni nodo della mappa

topologica costruita in precedenza, si associa un elemento particellare

dotato di massa. La piastra è stata discretizzata in un numero totale di

particelle pari a 11520. Essendo il numero di particelle coincidenti con

quello degli Element Solid, in cui è stata discretizzata la piastra con il

metodo FEM, si potrà operare il confronto tra i due metodi. Questi

elementi particellari si distribuiscono su cinque superfici sovrapposte,

ognuna delle quali costituita da 48x48 elementi particellari. Nel piano (x,

y) le particelle distano tra loro 2,5 cm in entrambe le direzioni, mentre

lungo l’asse z i cinque strati sono posti ad una distanza di 2,25 cm.

La massa d’ogni singola particella è stata semplicemente calcolata

dividendo il peso totale della piastra per il numero di particelle che la

costituiscono.


112

Ptot = γs Vtot = 2500 Kg/m3 (0,09 x1,2 x1,2)m3 = 324 Kg

⇓ Mele = Ptot / Numele = 324 Kg / 6912 = 0,046875 Kg

*ELEMENT_SPH $# nid pid mass 1 1 0.04687500

Tabella 5. 12 – Element SPH.

Ad ogni nid (Node Identification) è attribuita una particella di massa

(mass) pari a 0,046875 Kg, ed una pid (Part Identification), a cui

compete la stessa definizione data per la modellazione in FEM [12].

La differenza principale che sussiste tra la modellazione in FEM e quella

in SPH, risiede nella definizione della Section da associare all’elemento

modellato. Mentre per il metodo FEM si è attribuito all’elemento

modellato una Section_Solid, in questo caso sarà ovviamente necessario

definire una Section_SPH [12]. Questo comando consente di associare

alcune proprietà alle singole particelle.

*SECTION_SPH_TITLE SPH_section $# secid cslh hmin hmax sphini death start 1 1.200000 0.200000 2.000000 0.0001.0000E+20

Tabella 5. 13 – Section SPH.

Il parametro numerico di maggiore interesse che deve essere in questo

caso definito, è la Smoothing Lenght (h), di cui già si è detto nel capitolo

precedente.

Il processore LS-Dyna utilizza una Smoothing Lenght variabile durante il

tempo d’analisi. Il Processore calcola un valore iniziale di Smoothing

Lenght (h0), come la massima tra le minime distanze tra ogni particella.


113

Ad ogni particella è associato la propria Smoothing Lenght, che varia nel

tempo in accordo con la seguente equazione:

Dove h(t) è la Smoothing Lenght in funzione del tempo, mentre div(v) è

la divergenza del flusso. La Smoothing Lenght aumenta, quando le

particelle si allontanano le une dalle altre e si riduce, quando aumenta la

densità delle particelle. Esso varia in modo tale da mantenere costante il

numero di particelle tra loro vicine.

La Smoothing Lenght varia tra un valore minimo e uno massimo:

Non è possibile associare un valore unitario né a quello massimo né a

quello minimo, Hmin e Hmax, in quanto così facendo si definirebbe una

Smoothing Lenght costante al variare del tempo e dello spazio, e

coincidente con quella iniziale h0 calcolata dal processore [11].

Nella Card della Section_SPH, è necessario definire oltre al valore

minimo e massimo della Smoothing Lenght (hmin, hmax), anche una

costante cslh, che moltiplica la Smoothing Lenght delle particelle. Tale

valore deve essere compreso tra 1,05 e 1,3. Il valore di Default è uguale

1,2 ed è ritenuto accettabile per diverse problematiche. Valori inferiori

all’unità sono ritenuti inammissibili, mentre quelli maggiori di 1,3

aumenterebbero i tempi di computazione del software.

I parametri sphini, death e start, che definiscono la Card, individuano

rispettivamente un valore iniziale della Smoothing Lenght, ed il tempo

d’arresto e d’inizio dell’approssimazione. Nel primo caso se si definisce

)()())(( vdivththdtd

=

0max0min )( hHthhH ⋅≤≤⋅


114

un valore di sphini diverso dallo 0,0 si sovrascrive il valore iniziale h0

calcolato dal processore all’inizio del calcolo [11].

Una successiva differenza rispetto al modello costruito in FEM, riguarda

la modalità d’applicazione dei carichi. Questi ultimi ovviamente non

variano per le stesse prove che sono state modellate, al variare dei due

metodi impiegati FEM e SPH.

In questo caso essendo la piastra discretizzata in un numero finito

d’elementi nodali (Element_SPH), è preferibile applicare il carico

nodalmente. Si passa quindi da una condizione di carico distribuita ad

una di carico concentrato. Ogni singolo elemento particellare che

costituisce l’estradosso della piastra, sarà quindi soggetto ad una forza

concentrata, d’intensità variabile da una particella all’altra, e diretta

lungo l’asse z verso il basso.

Per definire questa condizione di carico è quindi necessario costruire dei

Set di nodi, allo stesso modo di quanto visto per la definizione delle

condizioni di vincolo. Si definiranno quindi 36 Set_Node_List [12], ad

ognuno dei quali apparteranno 64 nodi differenti, riconducibili alle 4

aeree di carico, simmetriche rispetto al centro della piastra, così come

descritte in precedenza.

*SET_NODE_LIST_TITLE A-1 $# sid da1 da2 da3 da4 1 $# nid1 nid2 nid3 nid4 nid5 nid6 nid7 nid8 6885 6886 6887 6888 6837 6838 6839 6840 4773 4774 4775 4776 4725 4726 4727 4728 4681 4682 4683 4684 4633 4634 4635 4636 4777 4778 4779 4780 4729 4730 4731 4732 4677 4678 4679 4680 4629 4630 4631 4632 6789 6790 6791 6792 6741 6742 6743 6744 6889 6890 6891 6892 6841 6842 6843 6844 6793 6794 6795 6796 6745 6746 6747 6748

Tabella 5. 14 – Set_Node_List Carico.


115

Per applicare il carico in modo concentrato su ogni particella, si usa il

comando Load_Noad_Set [12]. In questo modo si attribuisce ad ogni Set

di nodi la corrispondente Curva di Carico (Load Curve). Per tener conto

che si sta passando da una pressione, e quindi una forza per unità di

superficie, ad una forza concentrata, la curva di carico dovrà essere

moltiplicata per un fattore di scala, pari all’area di carico che compete ad

ogni singola particella, che sarà pari a (2,5x2,5) cm2.

*LOAD_NODE_SET $# nsid dof lcid sf cid m1 m2 m3 1 3 1 -6.250E-4 2 3 2 -6.250E-4 3 3 3 -6.250E-4 4 3 4 -6.250E-4 5 3 5 -6.250E-4 6 3 6 -6.250E-4 7 3 7 -6.250E-4 8 3 8 -6.250E-4 9 3 9 -6.250E-4 10 3 10 -6.250E-4 11 3 11 -6.250E-4 12 3 12 -6.250E-4 13 3 13 -6.250E-4 14 3 14 -6.250E-4 15 3 15 -6.250E-4 16 3 16 -6.250E-4 17 3 17 -6.250E-4 18 3 18 -6.250E-4 19 3 19 -6.250E-4 20 3 20 -6.250E-4 21 3 21 -6.250E-4 22 3 22 -6.250E-4 23 3 23 -6.250E-4 24 3 24 -6.250E-4 25 3 25 -6.250E-4 26 3 26 -6.250E-4 27 3 27 -6.250E-4 28 3 28 -6.250E-4 29 3 29 -6.250E-4 30 3 30 -6.250E-4 31 3 31 -6.250E-4 32 3 32 -6.250E-4 33 3 33 -6.250E-4 34 3 34 -6.250E-4 35 3 35 -6.250E-4 36 3 36 -6.250E-4

Tabella 5. 15 – Load_Node_Set.

Dove il dof indica la direzione in cui il carico agisce, il fattore di scala

invece si pone negativo perché il carico è applicato verso il basso.


116

Le altre Card che completano il modello in SPH equivalgono a quelle

definite per il modello FEM.

Per rendere più realistica la nostra modellazione, è stato considerato

anche l’effetto di Damping, impostando lo stesso al valore tipico per un

calcestruzzo del 5%. Questo comando consente, infatti, di considerare lo

smorzamento della velocità con la quale si muovono gli elementi.

Tra le varie possibilità si è scelto il Damping_Frequency_Range [12],

che consente di applicare l’effetto di smorzamento alla Part

dell’elemento. *DAMPING_FREQUENCY_RANGE $# cdamp flow fhigh psid 0.050000 0.100000 100.00000 1

Tabella 5. 16 – Damping.

Una volta assegnata la costante di Damping cdamp, si definisce un

range di frequenze, flow e fhigh, entro cui si considera l’effetto di

smorzamento.

CCAAPPIITTOOLLOO 66 –– LLEEGGAAMMII CCOOSSTTIITTUUTTIIVVII DDEELL CCAALLCCEESSTTRRUUZZZZOO

117

Capitolo 6 Legami Costitutivi del Calcestruzzo

Un aspetto cruciale di questa modellazione risiede nella definizione delle

caratteristiche meccaniche del calcestruzzo. Il Software di calcolo LS-

Dyna Version 970 implementa più di duecento Legami Costitutivi,

ognuno dei quali impiegati per la modellazione di materiali differenti. Si

ha quindi la possibilità di scegliere tra un’ampia gamma di Material

Model, secondo gli obiettivi che si vogliono perseguire e dei modelli che

si vogliono creare.

In maniera molto sintetica è possibile classificare i Material Model in

funzione del tipo di materiale da modellare, delle caratteristiche peculiari

ad essi associate, e dei metodi di discretizzazione con essi compatibili. In

particolar modo è possibile distinguere modelli finalizzati alla

modellazione di materiali plastici, compositi, ceramici, fluidi, vetrosi,

idrodinamici, metallici, gommosi, schiumosi ed infine terreni e

calcestruzzi.

Inoltre, nel caso specifico della modellazione realizzata era necessario

scegliere un Material Model che consentisse di considerare il così detto

effetto di Strain-Rate, e il raggiungimento della condizione di rottura

(Failure) [12].

Tra tutti i Material Model quelli che meglio rispondono a dette necessità

sono lo Pseudo Tensor (Mat_016) e il Concrete Damage (Mat_072).

6.1 – Pseudo – Tensor(Mat_016) Nell’analisi di strutture sottoposte a carichi da esplosione e per grandi

valori di deformazioni, i codici di calcolo Lagrangiani basati sulla teoria


118

degli elementi finiti, necessitano di uno strumento efficace per la

modellazione del comportamento del materiale.

In particolar modo per la modellazione d’elementi in calcestruzzo

armato, l’implementazione di un’efficiente Material Model si è

dimostrata complessa e variabile. Diversi modelli sono propriamente

adatti per descrivere il comportamento costitutivo del calcestruzzo, che

va dalla risposta elastica del materiale fino alla sua rottura. Tra questi si

annovera lo Pseudo Tensor, che risulta essere uno dei migliori.

Questo modello può essere impiegato definendo due superfici in

funzione della pressione media p, una di rottura massima σmax, e l’altra di

tensione residua σr [13].

paapa

210max ++=σ

paapar

ff

210 +

+=σ

Δσ

P

Figura 6. 1 – Superficie di rottura massima e residua.

È necessario associare al modello del materiale un equazione di stato,

essa fornisce il valore corrente della pressione p come una funzione della

deformazione volumetrica. In funzione della pressione si definisce una


119

superficie di snervamento o di rottura (Yeld o Failure), che limita il

secondo invariante della parte deviatorica del tensore tensionale.

La superficie superiore indica la tensione massima di snervamento del

materiale σmax, al variare della pressione p, mentre quell’inferiore la

tensione di rottura σfailed. Come si può osservare le due superfici sono

definite in funzione di cinque parametri a0, a1, a2, a0f e a1f, ai quali è

quindi necessario associare un valore numerico nella Card del materiale.

Il codice di calcolo impiegato ai fini della modellazione, suggerisce un

metodo semplice per la definizione di detti parametri, l’Internal

Generation [14].

Quando si utilizza questo tipo di materiale per la modellazione del

calcestruzzo, si ha la possibilità di definire i parametri della Card del

materiale, in funzione della sola resistenza a compressione del

calcestruzzo non confinato f’c. Questa viene, infatti, assunta come la

tensione di rottura a compressione σf. Tutto questo è reso possibile

assumendo semplicemente un valore negativo per il parametro a0, ad

esempio a0 = - 1 [14].

In questo modo il parametro –a0 è considerato come un fattore di

conversione tra l’unità di misura usata dal LS-Dyna per la pressione, e la

stessa espressa in psi [14].

Inoltre, questo Material Model consente di considerare anche la presenza

di una percentuale di armatura di rinforzo, di cui se ne definiscono le

caratteristiche meccaniche. Per i suddetti motivi, il Mat_016 è

raccomandato per l’analisi d’elementi strutturali in cemento armato

soggetti a carichi impulsivi.

A titolo d’esempio, si riporta di seguito la Card del programma LS-Dyna

attraverso il quale si definisce le caratteristiche del materiale, impiegato


120

nella modellazione della piastra, qualora si volesse usare l’opzione

Internal Generation: *MAT_PSEUDO_TENSOR_TITLE C.A. $# mid ro g pr 1 2500.0000 7.3486E+9 0.200000 $# sigf a0 a1 a2 a0f a1f b1 per 9.7896E+6 -1.000000 0.000 0.000 0.000 0.000 0.000 0.527980 $# er prr sigy etan lcp lcr 1.9990E+11 0.300000 3.8820E+8 $# x1 x2 x3 x4 x5 x6 x7 x8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 $# x9 x10 x11 x12 x13 x14 x15 x16 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 $# ys1 ys2 ys3 ys4 ys5 ys6 ys7 ys8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 $# ys9 ys10 ys11 ys12 ys13 ys14 ys15 ys16 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

Figura 6. 2 – Card Material Model_016. Dove:

Calcestruzzo Acciaio ro Densità di massa per % di Armatura g Modulo di elasticità a taglio er Modulo di Young pr Rapporto di Poisson prr Rapporto di Poisson sigf Tensione di rottura a compressione sigy Tensione di snervamento a0 Coesione etan Tensione tangenziale

Tabella 6. 1 – parametri che definiscono la Card del Material Model.

Nel momento in cui al parametro a0 si associa un valore negativo, il

numero identificativo da associare all’equazione di stato (EOSid) nella

Part del modello, si pone uguale a zero. In questo modo sarà

automaticamente generata un equazione di stato Tabulated Compaction

(EOS 8) in funzione dei soli valori f’c, e del Rapporto di Poisson ν.

Questo tipo d’approccio è troppo semplificativo, e non consente di

controllare quanto fatto in automatico dal software. Per questo motivo è

preferibile assegnare tutti i parametri richiesti dalla Card del materiale,

ad ognuno dei quali si associa un significato fisico.

Lo Pseudo Tensor considera in maniera dissociata la risposta deviatorica

e quella sferica del calcestruzzo durante la fase di carico e ricarico.


121

Durante la fase di carico iniziale o di ricarico, la parte deviatorica si

mantiene elastica fino a raggiungere la superficie iniziale di snervamento

(Yeld Surface), raggiunta la quale può ancora aumentare fino a

raggiungere la superficie di rottura massima (Max Surface). Quando il

provino di calcestruzzo è poi scaricato, la parte deviatorica del tensore

delle tensioni raggiunge la superficie di rottura residua (Residual

Surface). Così come mostrato nella seguente figura [13].

Figura 6. 3 – Andamento della superficie di rottura corrente Δσ, in funzione della tensione sferica p.

La superficie di rottura corrente Δσ che limita la parte deviatorica dello

stato tensionale, è definita come una combinazione lineare di due

superfici, Δσm e Δσr, a loro volta definiti in funzione della tensione

sferica p [13].

rm σησησ Δ−+Δ=Δ )1(


122

Dove:

Con Δσm s’indica la superficie di rottura massima, mentre con Δσr quella

residua [13].

Nelle precedenti espressioni compare la tensione sferica p, che a sua

volta è definita in funzione delle tensioni normali:

Il parametro η indica la funzione di migrazione tra la superficie di

rottura residua e quella massima, ed è definito in funzione di λ,

quest’ultimo è una misura della deformazione plastica effettiva, ed

indica il danno accumulato dal materiale. La funzione η = η(λ) è di tipo

lineare, e va da 0 a 1, e da 1 a 0 [13].

La superficie di snervamento iniziale Δσy (Yeld Surface), è compresa tra

quella di rottura massima e residua, Δσm e Δσr, e si esprime come segue:

Dove ηy = η(0), in altre parole è il valore iniziale di η prima ancora che

abbia inizio la plasticizzazione del materiale, i risultati provenienti da

prove di compressione tri-assiale sul calcestruzzo, hanno mostrato che la

superficie di snervamento iniziale Δσy, dovrebbe essere circa il 45% di

quella massima. Dato che le due superfici diventano parallele per grandi

paapam

210 +

+=Δ σ

paapa

ffr

210 +

+=Δ σ

( )3

zyxpσσσ ++

−=

rymyy σησησ Δ−+Δ=Δ )1(


123

valori di p, esse non possono propriamente rappresentare il punto di

transizione tra il comportamento duttile e quello fragile.

Di seguito si descrive il significato delle superfici di rottura[13]:

Figura 6. 4 – Superfici di Rottura.

Compressive Meridian

La Compressive Meridian indica l’andamento della Δσm, i dati per la sua

costruzione si ottengono solitamente attraverso prove sperimentali

condotte su provini di calcestruzzo sottoposti a compressioni assiali non

confinati, e prove di compressione tri-assiali con vari livelli di

confinamento. I test usuali non forniscono in ogni modo risultati per

valori di pressione inferiori a f’c/3. Per il modello originale Pseudo

Tensor, sono necessari almeno due livelli di confinamento diversi da

zero affinché la Compressive Meridian sia definita da tre parametri a0, a1

e a2. Questi tre parametri che definisco la superficie di rottura massima


124

Δσm, solitamente sovrastimano gli sforzi, quando sono estrapolati per

valori inferiori a f’c/3.

Tensile Meridian

È noto che la Tensile Meridian della superficie di rottura per il

calcestruzzo è solitamente inferiore rispetto alla Compressive Meridian.

Dati sperimentali suggeriscono che il rapporto tra la Tensile e

Compressive Meridian, denotato da ψ, è compreso tra 0,5, quando le

pressioni sono negative, e 1, per livelli di confinamento elevati.

Tensile Cutoff

La Tensile Cutoff è inserita proprio per cercare di limitare i difetti del

modello per valori bassi di pressione, citati nella definizione della

Compressive Meridian. La Tensile Cutoff limita la tensione principale

massima allo sforzo di trazione ft. Ciò in ogni modo non risolve il

problema per valori di tensione sferica compresa tra 0 e f’c/3.

Pressure Cutoff

Il modello originale incorpora una Pressure Cutoff, che impedisce alla

pressione di andare sotto il valore di ft/3. La Pressure Cutoff, non incide

per le prove di trazione monoassiale, mentre limita la differenza di

tensione principale al valore di ft/2 per una prova di trazione biassiale, e

a ft/3 per una tri-assiale. Entrambi questi limiti contrastano con le

evidenze sperimentali, le quali mostrano che in entrambi i casi, la

differenza di tensione principale Δσ dovrebbe raggiungere

approssimativamente ft. Quantunque la Pressure Cutoff sia raggiunta, lo

stato delle tensioni corrente rimane in ogni modo inalterato.


125

Il comportamento elastico del materiale è descritto con lo Pseudo Tensor

mediante la definizione dell’equazione di stato Tabulated Compaction

(EOS8). Infatti, assegnando nell’EOS8 l’Unloading Bulk Modulus che è

una costante di elasticità definita in funzione del Modulo di Young E, e

del Rapporto di Poisson ν :

automaticamente in funzione di K il programma definisce un'altra

costante d’elasticità, che è il Modulo di elasticità a Taglio G.

Ciononostante, lo Pseudo Tensor contiene una serie di difetti, tra cui

l’incapacità di considerare la dilatazione a taglio, normalmente osservata

sui provini di calcestruzzo sottoposti a sollecitazioni di compressione.

Tali difetti sono stati corretti mediante l’introduzione del nuovo modello

Concrete Damage (Mat_072) [13].

6.2 – Concrete Damage (Mat_072) Una prima differenza rispetto al modello originale Pseudo Tensor risiede

nella definizione di una nuova Pressure Cutoff pc. Per valori negativi di

pressione, la funzione di migrazione η, serve non solo per scalare la

superficie di rottura corrente da quella massima a quella residua, ma

anche per aumentare la pc da –ft a 0. Ciò avviene controllando la

pressione calcolata dall’equazione di stato, e facendola coincidere con pc

se non sarà soddisfatta la condizione secondo la quale deve essere p ≥ pc

[13].

( )ν213 −=

EK


126

Dove pc è uguale a –ft, se non è raggiunta la superficie di rottura

massima, ed è uguale a –ηft, se invece è raggiunta. Questa modifica può

quindi sovrascrivere la pressione calcolata dall’equazione di stato[13].

L’altra innovazione riguarda le superfici di rottura. Nel Concrete

Damage, infatti, è inserita una terza superficie, indipendente dalle

precedenti due. Tale superficie rappresenta quell’iniziale di snervamento,

ed è definita in funzione di tre nuovi parametri a0y, a1y e a2y, così come

segue [13]:

Per quanto concerne le restanti due superfici, già viste nel modello

precedente, anche in questo caso sono apportate delle modifiche. Dato

che il calcestruzzo ha uno sforzo residuo in trazione nullo, si assume tale

anche il parametro indipendente dalla pressione p che definisce la

superficie di rottura residua, aof = 0. Sempre a proposito della superficie

di rottura residua, per considerare l’intersezione tra questa e quella di

rottura massima in un punto di transizione tra il comportamento duttile e

quello fragile, s’inserisce un nuovo parametro a2f, mentre nel modello

originale questo parametro era uguale per entrambe le superfici.

In conclusione la superficie residua adesso assume la seguente

espressione [13]:

paapa

yyyy

210 +

+=Δ σ

paap

ffr

21 +=Δ σ

⎩⎨⎧

−−

=t

tc f

fp

η


127

Nel nuovo modello, la superficie di rottura corrente dopo aver raggiunto

quella di snervamento iniziale, e prima ancora di raggiungere quella di

rottura massima, è definita come una combinazione lineare delle due:

Dove il parametro η è sempre compreso tra 0 e 1, ed è ancora definito

come una funzione della misura di deformazione plastica effettiva λ.

Una volta raggiunta la superficie massima, quella corrente di rottura è

ottenuta come interpolazione tra quella massima e quella residua:

Nel nuovo modello, la funzione di danno η = η(λ), è assegnata in

maniera tabulata dall’utente come una serie di valori (η, λ). Al variare di

λ si ottengono valori di η compresi comunque tra 0 e 1. Quando λ = 0

anche η = 0, quando λ = λm allora η = 1, e torna ad annullarsi, quando λ

assume valori maggiori di λm [13].

Finche λ non decresce, Δσ tende ad assumere i valori di Δσy, Δσm e Δσr.

Quando λ ≤ λm, allora Δσ si ottiene come un’interpolazione lineare tra

Δσy e Δσm, viceversa quando λ > λm, allora sarà data dall’interpolazione

tra Δσm e Δσr [13].

La superficie di rottura massima non subisce nessuna modifica rispetto a

quella vista nel modello originale.

Riassumendo, per il nuovo modello le tre superfici di rottura, quella

massima, residua e di snervamento, sono rispettivamente definite come

segue [13]:

( ) yym σσσησ Δ+Δ−Δ=Δ

( ) rrm σσσησ Δ+Δ−Δ=Δ


128

In prossimità del punto in cui si osserva la transizione tra il

comportamento fragile e quello duttile del materiale, avremo che sia la

Δσy sia la Δσr, saranno limitate dalla Δσm [13].

6.2.1 – Parametri che definiscono le superfici di rottura Vediamo come si determinano i parametri a0, a1 ed a2 che definiscono la

superficie di rottura massima Δσm. Questi tre parametri dovrebbero

essere determinati mediante prove di laboratorio su provini di

calcestruzzo sottoposti a sollecitazioni di compressione assiale, e prove

di compressione tri-assiale per diversi livelli di confinamento.

Purtuttavia è possibile ottenere una stima di questi tre valori anche

numericamente. Il primo dei tre, a0, rappresenta, infatti, il punto

d’intersezione della superficie di rottura massima con l’asse delle

ordinate, Δσ(p=0) = a0 [13].

Il secondo parametro, a1, si definisce come l’inverso del minimo della

funzione Δσ:

10

1adp

d

p

=⎥⎦

⎤⎢⎣

⎡Δ

=

σ

paapam

210 +

+=Δ σ

paap

ffr

21 +=Δ σ

paapa

yyyy

210 +

+=Δ σ


129

Infine, per grandi valori di tensione sferica p, il denominatore del

secondo termine dell’equazione che definisce Δσm, è governato dal

termine a2p, in quanto a1 diventa insignificante:

Per definire la superficie di snervamento iniziale Δσy, è invece

necessario determinare i valori dei parametri a0y, a1y e a2y. Da dati

sperimentali è possibile affermare che tale superficie sia il luogo dei

punti per Δσ = 0,45Δσm, per un percorso di compressione di tipo tri-

assiale. Ciò vuol affermare che preso un generico punto di coordinate (p,

Δσm), appartenente alla superficie di rottura massima, il corrispondente

punto sulla superficie iniziale di snervamento (p’, Δσy) sarà definito

come:

Da cui p può essere espresso in funzione di p’, mentre nota l’espressione

di Δσm, Δσy sarà calcolato come [13]:

Quindi, assumendo p’ = 0, avremo che Δσy = a0y, per ricavare gli altri

due parametri sarà sufficiente risolvere un sistema di due equazioni nelle

due incognite a1y e a2y.

Quando non si conoscono le superfici di rottura per un nuovo

calcestruzzo da modellare, per il qual è nota la sola resistenza a

compressione del cls non confinato f’c,new, i parametri a0n, a1n ed a2n, che

∞→→−Δ pquandoa

a2

01σ

mmy ppe σσσ Δ−=Δ=Δ355,045,0 '

⎟⎟⎠

⎞⎜⎜⎝

⎛+

+=Δpaa

pay21

045,0σ


130

definiscono la nuova superficie di rottura Δσn, possono essere calcolati a

partire dai parametri noti di un altro calcestruzzo f’c,old. In questo caso

viene suggerito di definire un fattore di scala r come il rapporto tra le

resistenze a compressione del calcestruzzo non confinato, del nuovo e

del vecchio calcestruzzo [13]:

I nuovi parametri sono quindi calcolati come:

Allo stesso modo si può pensare di calcolare i parametri che definiscono

l’andamento della Yeld Surface (a0y, a1y ed a2y), mentre i coefficienti che

determinano la Failure Surface (a1f ed a2f), si pongono uguali ad a1n ed

a2n [13].

A tal fine si fa riferimento alle caratteristiche meccaniche del

calcestruzzo impiegato per lo svolgimento di una prova sperimentale

condotta presso un centro di ricerca Svedese, per la quale si prevedeva

un blocco di calcestruzzo ad elevata resistenza, sottoposto ad un carico

da impatto [19]. Le prove condotte sono state in seguito modellate con il

codice di calcolo LS-Dyna Version 960. Note le caratteristiche

meccaniche del calcestruzzo modellato, e in particolar modo la sua

resistenza a compressione f’c,old, ed altresì note i parametri assegnati nella

Card del Material Model Concrete Damage, alla luce di quanto detto in

precedenza si sono calcolati i parametri relativi al calcestruzzo da

modellare f’c,new:

raaaa

raa

n

n

n

/22

11

00

===

',

',

oldc

newc

ff

r =


131

Caratteristiche dei Materiali Old Material New Material

Ec[Gpa] 58 Ec[Gpa] 17636643,45

Ec [Mpa] 58000 Ec [Mpa] 17636,64

f'c,old[Mpa] 153 f'

c,new [Mpa] 13,79

f'c,old[Pa] 153000000 f'

c,new [Pa] 13786202,43

Fattore di Scala

r 0,090105898

Parametri LS-Dyna Version 970

a0,old 5,06E+07 a0,new 4,56E+06

a1,old 0,465 a1,new = a1f,new 0,465

a2,old 6,57E-10 a2,new = a2f,new 7,29E-09

a0y,old 2,28E+07 a0y,new 2,05E+06

a1y,old 1,033 a1y,new 1,033

a2y,old 1,46E-09 a2y,new 1,62E-08

Tabella 6. 2 – Parametri numerici che definiscono le superfici di rottura d’input per LS-Dyna.

Durante la fase di carico, la superficie di rottura corrente Δσ, migra tra la

superficie di rottura massima e una tra quella di snervamento iniziale e di

rottura residua [13].

Quando λ ≤ λm, allora Δσmin coincide con Δσy, viceversa, quando λ > λm

in tal caso Δσmin coincide con Δσr. Allo stesso modo che per il Modello

Originale Pseudo Tensor, anche in questo caso η si esprime come una

funzione di λ, che indica la deformazione plastica effettiva. L’unica

differenza risiede nell’espressione di λ, che in questo caso varia a

seconda che la parte sferica del tensore delle tensioni p, sia maggiore o

minore di 0 [13].

( ) minminmax σσσησ Δ+Δ−Δ=Δ


132

Nell’originale Material Model (Mat_016), la misura della deformazione

plastica effettiva λ, è definita come segue [13]:

Dove l’incremento di deformazione plastica effettiva è dato da:

Nel Concrete Damage sono apportate due modifiche nel definire il

parametro λ. Innanzitutto, si introduce un fattore sperimentale rf, che

tiene conto del miglioramento di velocità (rate enhancement), inoltre il

parametro b1 è sostituito da b2 per valori di pressione minori di 0 [13]:

Quando la pressione varia tra 0 e –ft, dove con ft s’indica la resistenza a

trazione del calcestruzzo, il parametro λ rimane nullo, così come η.

L’equazione di stato riduce la pressione a –ft e la mantiene costante. Per

implementare la caduta di pressione, in seguito alla rottura per trazione,

può essere aggiunto al danno deviatorico preesistente un incremento di

danno volumetrico, ogni qual volta il percorso di tensione è chiuso al

percorso che si osserva in un test di trazione tri-assiale. Per limitare gli

∫⎟⎟⎠

⎞⎜⎜⎝

⎛+

=p

b

t

p

fp

dε ελ

01

1

pij

pijpd εεε )3/2(=

⎪⎪⎪⎪

⎩

⎪⎪⎪⎪

⎨

⎧

<

⎟⎟⎠

⎞⎜⎜⎝

⎛+

≥

⎟⎟⎠

⎞⎜⎜⎝

⎛+

=

∫

∫

p

p

p

frpr

d

p

frpr

d

b

tf

f

p

b

tf

f

p

ε

ε

ε

ε

λ

0

0

;0

1

;0

1

2

1


133

effetti di questo cambiamento, l’incremento del danno è moltiplicato da

un fattore fd, definito come segue [13]:

Dove è il rapporto che misura la ristrettezza di tale

percorso, in cui p è la tensione sferica, mentre , ed è

la superficie di rottura per le tensioni deviatoriche. Il termine J2, che ivi

compare nell’espressione di Δσ, è il secondo invariante del tensore delle

tensioni, così definito:

La misura della deformazione plastica effettiva modificata per il

Concrete Damage, è incrementata da un Δλ, definito come segue [13]:

Dove b3 è un parametro definito dall’utente, di cui si dirà in seguito, kd è

un fattore scalare interno, εv è la deformazione volumetrica, e infine

εv,yeld è la deformazione volumetrica quando si raggiunge lo snervamento

del materiale [13]:

⎪⎩

⎪⎨

⎧

≥

<≤−=

1,030

1,0/301,0

/31

2

22

Jse

pJsepJ

f d

( )yeldvvdd kfb ,3 εελ −=Δ

pJ /3 2

σΔ=23 J

2

23

22

2

21

σσσ ++=J

zvt

vd f

kk εεε31;

3−=−=


134

Per determinare l’andamento della funzione di danno η = η(λ), che

consente di definire la superficie di rottura corrente Δσ, è quindi

necessario calcolare i tre parametri b1, b2 e b3 [13].

I parametri b2 e b3 governano l’andamento del ramo di softening della

curva tensio-deformativa del materiale, per un calcestruzzo sottoposto ad

una prova di trazione rispettivamente monoassiale e tri-assiale, quando la

tensione si muove dalla superficie di rottura massima a quella residua.

L’apertura delle fessure sul modello, sono in ogni modo dipendenti dalle

dimensioni e dalle caratteristiche della Mesh. Un modo per eliminare

questa dipendenza è imporre che l’area sottesa alla curva σ−ε sia uguale

a Gf/h. Dove Gf è l’energia di frattura, ed h la dimensione del singolo

elemento che costituisce la Mesh. Quando la localizzazione della fessura

capita in un elemento, allora h coincide con wc, dove wc è l’ampiezza

della fessura, ed è tipicamente contenuta tra 1 e 6 volte la dimensione

massima degli inerti del calcestruzzo[13]. Invece, in accordo con i codici

Europei (CEB, FIP) l’energia di frattura varia tra 40 e 175 N/m [22].

Questo modello assicura un’indipendenza dalle dimensioni della Mesh,

che consente di apprezzare anche l’apertura delle fessure sul modello.

Ciononostante, la mappa topologica della Mesh potrebbe incidere

sull’orientamento delle fessure, sia in condizioni statiche sia dinamiche.

Il parametro b2 è definito facendo variare l’andamento della curva σ−ε,

definita per una prova di trazione monoassiale, fino a che l’area sottesa

dalla curva non coincide con Gf/h. Lo stesso procedimento iterativo è

eseguito anche nella determinazione di b3, con l’unica differenza che in

questo caso si fa affidamento ai dati sperimentali provenienti da una

prova di trazione tri-assiale idrostatica e non monoassiale. Questi due


135

parametri sono di fondamentale importanza, quando l’elemento

analizzato è a debole armatura.

Il parametro b1 è determinato allo stesso modo che per i precedenti, in

questo caso il coefficiente numerico è tarato su dati sperimentali

provenienti da test a compressione monoassiale al variare del livello di

confinamento. Alla luce di quanto detto, la superficie di rottura corrente

Δσ dipende indirettamente dai tre parametri b1, b2 e b3 [13].

È possibile, infatti, schematizzare le seguenti dipendenze:

Δσ = f(η) → η = η(λ) → λ = f(b1 o b2), mentre l’incremento del

parametro di danno Δλ = f(b3).

A loro volta b1, b2 e b3 si esprimono come una funzione di Gf e h.

Queste dipendenze sono implementate dallo stesso codice di calcolo LS-

Dyna Version 970, al quale è necessario assegnare tra gli input che

definiscono le caratteristiche del materiale, oltre ai parametri b1, b2 e b3,

anche in maniera tabulata la curva di danno η = η(λ).

Non avendo la possibilità di eseguire prove sperimentali finalizzate alla

caratterizzazione dei materiali, si è costretti ad assumere per questi tre

parametri b1, b2 e b3, tipici valori definiti per un calcestruzzo di scarsa

qualità, al quale compete un basso valore d’Energia di Frattura Gf, così

come per il calcestruzzo impiegato nello svolgimento delle prove

sperimentali modellate. A detti valori sarà associata la corrispondente

funzione di danno, definita così come descritto in precedenza:


136

Funzione di Danno η = η(λ)

0

0,2

0,4

0,6

0,8

1

1,2

0 0,01 0,02 0,03 0,04 0,05

Grafico 6. 1 – Funzione di Danno, calcolata per b1 = b2 = b3 = 0,1.

Tra i due Material Model in precedenza descritti, Pseudo Tensor e il

Concrete Damage, per la modellazione del calcestruzzo si è deciso di

impiegare questo ultimo. Essendo, infatti, il Concrete Damage un

evoluzione dello Pseudo Tensor.

Di seguito per semplicità d’esposizione si riporta la Card del Material

Model così come è stata definita per la modellazione della piastra:

*MAT_CONCRETE_DAMAGE_TITLE C.A. $# mid ro pr 1 2500.0000 0.200000 $# sigf a0 a1 a2 1.3790E+6 4.68E+06 0.465 7.11E-09 $# a0y a1y a2y a1f a2f b1 b2 b3 2.11E+06 1.033 1.58E-08 0.465 7.29E-09 0.100000 0.100000 0.100000 $# per er prr sigy etan lcp lcr 0.00528 1.999E+11 0.300000 3.8820E+8 38 39 $#lambda-1 lambda-2 lambda-3 lambda-4 lambda-5 lambda-6 lambda-7 lambda-8 0.000 0.02E-3 2.80E-03 41.0E-03 $#lambda-9 lambda-10 lambda-11 lambda-12 lambda-13 0.000 0.000 0.000 0.000 0.000 $# nu-1 nu-2 nu-3 nu-4 nu-5 nu-6 nu-7 nu-8 0.000 1.000000 0.150000 0.000

Tabella 6. 3 – Input per LS-Dyna del material Concrete Damage, definiti per la modellazione della piastra.


137

Dove, oltre ai parametri in precedenza descritti e calcolati, ao, a1, a2, aoy,

a1y, a2y, a1f, a2f, b1, b2 e b3, sono presenti altri parametri che definiscono le

caratteristiche del calcestruzzo:

• mid è il numero identificativo del Material Model;

• ro è il peso specifico del calcestruzzo;

• pr è il Rapporto di Poisson;

• sigf è la tensione di rottura a trazione del calcestruzzo, calcolata

come un decimo di quella a compressione f’c.

Le caratteristiche dell’armatura meccanica sono invece:

• per indica la percentuale di armatura nella sezione;

• er è il Modulo di Elasticità longitudinale di Young;

• prr è il Rapporto di Poisson;

• sigy è la tensione di snervamento dell’acciaio in condizione di

carico quasi statiche.

Con lcp (load curve principal) e lcr (load curve rinforcement), si

richiamano rispettivamente le curve di Strain-Rate relative al materiale

principale, il calcestruzzo, e al materiale di rinforzo, l’acciaio.

Queste curve sono definite così come di seguito descritto [12].

6.4 – Dynamic Increase Factor (DIF) Il Dynamic Increase Factor (DIF), è un coefficiente numerico

adimensionale, che moltiplicato per la resistenza meccanica del

calcestruzzo, calcolata in condizioni di carico quasi statico, fornisce una


138

misura delle caratteristiche dinamiche, al variare della velocità di

deformazione ε (Strain Rate). La velocità di deformazione non è altro

che la variazione della deformazione al variare del tempo [18].

Per la simulazione numerica, è di fondamentale importanza definire

l’andamento del DIF in funzione della velocità di deformazione.

Trattandosi di un elemento strutturale in cemento armato, sottoposto

all’azione dinamica prodotta da un’esplosione, si avrà una velocità di

deformazione molto elevate, comprese tra 10 sec-1 e 1000 sec-1 [17]. Per

valori così elevati di Strain-Rate, le evidenze sperimentali hanno

mostrato un incremento di resistenza del calcestruzzo sia a compressione

ma soprattutto a trazione. Il DIF a compressione in alcuni casi è

maggiore di 2, mentre per la resistenza a trazione può raggiungere valori

maggiori di 6 [17].

Per l’armatura metallica di rinforzo, si può arrivare ad un incremento di

resistenza maggiore del 50% [17].

Questo incremento di resistenza determina quindi un miglioramento del

comportamento meccanico d’entrambi i materiali. Inoltre, essendo la

risposta dell’elemento strutturale in c.a. molto dipendente dal

comportamento a trazione del calcestruzzo, l’aumento di resistenza a

trazione, enfatizza l’importanza che assume la descrizione del

comportamento dinamico del materiale.

Volendo quindi considerare il cosiddetto effetto di Strain-Rate, nel

comportamento a compressione e a trazione del calcestruzzo, sottoposto

ad un carico dinamico impulsivo, si è scelto di simulare numericamente

la prova sperimentale impiegando il software LS-Dyna [11]. Infatti, tra i

diversi legami costitutivi che si possono implementare con questo codice

di calcolo strutturale, alcuni di questi forniscono la possibilità di


139

descrivere il comportamento costitutivo di un materiale in funzione della

velocità di deformazione.

Al variare del Material Model, scelto per la modellazione del materiale,

il programma fornisce tre diverse opzioni nell’assegnare l’effetto di

Strain-Rate [12]:

• Cowper and Symonds.

Questa prima scelta si fonda sulla teoria di Cowper and Symonds. In

questo caso si assegnano nella Card che definisce le caratteristiche

del materiale, due costanti C e p. In funzione delle quali il programma

calcola un fattore amplificativo (DIF) del tipo:

p

C

1

1⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛+

•

ε

Al variare della velocità di deformazione ε , il fattore amplificativo

che ne consegue, è moltiplicato per la resistenza meccanica del

materiale, migliorandone il comportamento sotto l’azione di carichi

dinamici.

• Legame Tensione-Deformazione Parametrico.

In questo caso, durante la fase d’input, s’inseriscono i tratti plastici

del legame costitutivo del materiale, ricavati per diversi valori di ε ,

formando così una famiglia di legami costitutivi. In questo modo, al

variare della velocità di deformazione, il programma attinge da

legami costitutivi diversi, il valore della corrispondente tensione.


140

• Curva di Strain-Rate.

Il modo più semplice per assegnare l’effetto di Strain-Rate, consiste

nel definire una curva che lega alla velocità di deformazione il

corrispondente DIF. Definita e in seguito assegnata la Curva di

Strain-Rate alla Card del Material Model, al variare della velocità di

deformazione, il programma estrapola da questa curva il

corrispondente valore del DIF. Questo ultimo moltiplica la tensione

di rottura a trazione σf, calcolata in condizioni di carico Quasi-

Statico, ed assegnata nella Card che definisce le caratteristiche

meccaniche del materiale. Così facendo al variare della velocità di

deformazione, e quindi del corrispondente DIF, si avranno valori

variabili di resistenza a trazione.

Sia lo Pseudo-Tensor sia il Concrete-Damage, che si è scelto di

impiegare nella modellazione del calcestruzzo, utilizzano quest’ultima

possibilità nell’assegnare l’effetto di Strain-Rate. Inoltre, mediante

l’assegnazione di due curve ε –DIF, definite rispettivamente per il

calcestruzzo e per l’acciaio, si avrà la possibilità di considerare il

fenomeno in maniera disaccoppiata per i due materiali, tra i quali si

considera il calcestruzzo principale e l’armatura di rinforzo secondaria.

Ovviamente, ciò presuppone la conoscenza delle curve di Strain-Rate

per entrambi i materiali.

6.4.1 – Definizione della Curva di Strain-Rate Nel definire l’andamento del DIF in funzione della velocità di

deformazione, si è scelto di utilizzare il modello analitico proposto dal

CEB-FIP Model Code 1990 [20]. Questo modello consente di calcolare,


141

in funzione della velocità di deformazione, il valore numerico del DIF

per il calcestruzzo, sia a trazione sia a compressione.

6.4.2 – Origini della formulazione CEB-FIP La formulazione proposta per il calcolo del DIF in funzione della

velocità di deformazione, nasce ovviamente da osservazioni sperimentali

e modelli numerici. Queste prove, condotte per caratterizzare il

comportamento a compressione e a trazione del calcestruzzo, sotto

l’azione di carichi dinamici, hanno evidenziato come l’effetto di Strain-

Rate sul calcestruzzo, sia a trazione sia a compressione, si possa

esprimere attraverso un fattore adimensionale (DIF), inteso come il

rapporto tra lo sforzo dinamico e quello statico. Le prove sono state

condotte sottoponendo un provino di calcestruzzo cilindrico, di cui sono

note a priori le caratteristiche meccaniche, per una velocità di

deformazione Quasi-Statica, ad un carico dinamico-impulsivo. A diversi

intervalli di velocità di deformazione si misurano valori differenti di

sforzi di trazione o compressione dinamici. Il DIF è semplicemente

calcolato come rapporto tra le due caratteristiche meccaniche, quella

dinamica e quella statica. Operando un confronto tra i dati sperimentali

disponibili, e i risultati conseguiti nel definire la curva di Strain-Rate con

il modello teorico raccomandato dal CEB-FIP, si è giunti alla

calibrazione del Model Code 1990, secondo il quale il DIF per sforzi di

compressione per il calcestruzzo sarà dato da [20].

α

εε

026.1

⎟⎟⎠

⎞⎜⎜⎝

⎛=

scs

c

ff per ≤ε 30 s-1

31

⎟⎟⎠

⎞⎜⎜⎝

⎛=

ss

cs

c

ff

εεγ per ε > 30 s-1


142

Dove:

• fc è lo sforzo di compressione dinamico alla velocità di

deformazione ε;

• fcs è lo sforzo di compressione statico alla velocità di

deformazione sε ;

• cs

c

ff

è il Dynamic Increase Factor (DIF);

• ε è la velocità di deformazione compresa tra 30 x 10-6 sec-1 e

300 sec-1;

• sε = 30 x 10-6 sec-1, è la velocità di deformazione quasi statica a

compressione;

• 2156.6log −= αγ s ⇒ 492.0156.610 −= αγ s ;

•

⎟⎟⎠

⎞⎜⎜⎝

⎛+

=

co

cs

ff

95

1α ; fco = 10 Mpa

Questa formulazione analitica, valida per il solo calcestruzzo in

compressione, è stata ricavata su basi sperimentali.

Invece, il DIF per gli sforzi di trazione è dato da [20]:

δ

εε

016.1

⎟⎟⎠

⎞⎜⎜⎝

⎛=

sts

t

ff per ≤ε 30 s-1

31

⎟⎟⎠

⎞⎜⎜⎝

⎛=

sts

t

ff

εεβ per ε > 30 s-1


143

Dove:

• ft è lo sforzo di trazione dinamico alla velocità di

deformazione ε;

• fts è lo sforzo di trazione statico alla velocità di deformazione

sε ;

• ts

t

ff

è il Dynamic Increase Factor (DIF);

• ε è la velocità di deformazione compresa tra 3 x 10-6 sec-1 e

300 sec-1;

• sε = 3 x 10-6 sec-1, è la velocità di deformazione quasi statica a

trazione;

• 33.211.7log −= δβ ⇒ 33.211.710 −= δβ ;

•

⎟⎟⎠

⎞⎜⎜⎝

⎛+

=

co

ts

ff

610

1δ ; fco = 10 Mpa

Da queste espressioni si evince che, la relazione che sussiste tra il DIF e

la velocità di deformazioneε, riportate in un grafico e in scala

logaritmica, è di tipo bilineare, con un cambio d’inclinazione che si ha

per ε = 30 sec-1. Inoltre, più scadente è la qualità del calcestruzzo

adottato, e maggiore sarà il valore del DIF, essendo sempre più piccolo il

valore dello sforzo di trazione o di compressione, misurato per una

condizione di carico Quasi-Statica [20].


144

Questo modello è stato impiegato per la simulazione numerica, nelle

pagine a seguire si riportano in maniera schematica i risultati così

ottenuti:

Compres. εs' [sec-1] 30x10-6

Trazione εs' [sec-1] 3x10-6

Tabella 6. 4 – Velocità di deformazione quasi statica.

fcs [Mpa] 13,79

fcs [Mpa] 1,38

fco [Mpa] 10,00 α 0,057 γ 0,023 δ 0,092 β 0,021

Tabella 6. 5 – Caratteristiche statiche del calcestruzzo e parametri numerici utilizzati per il calcolo delle curve di Strain-Rate.

Si riportano di seguito le Curve di Strain-Rate impiegate nella

modellazione:

DIF Calcestruzzo

ε' [sec-1] DIF compressione DIF trazione

0,000003 ------ 1,000 3,00E-05 1,000 1,241 3,00E-04 1,145 1,541 3,00E-03 1,312 1,912 3,00E-02 1,503 2,373

0,3 1,721 2,946 3 1,971 3,656 6 2,053 3,902 9 2,103 4,053 12 2,139 4,164 15 2,167 4,252 18 2,191 4,326 21 2,211 4,389 24 2,228 4,444 27 2,244 4,493 30 2,258 4,538 32 2,307 4,670


145

34 2,354 4,766 36 2,399 4,857 38 2,443 4,946 40 2,485 5,031 42 2,526 5,114 44 2,565 5,193 46 2,603 5,271 48 2,640 5,346 50 2,677 5,420 52 2,712 5,491 54 2,746 5,560 56 2,780 5,628 58 2,812 5,694 60 2,844 5,759 62 2,876 5,822 64 2,906 5,884 66 2,936 5,945 68 2,966 6,005 70 2,994 6,063 72 3,023 6,120 74 3,050 6,176 76 3,078 6,231 78 3,104 6,285 80 3,131 6,339 82 3,156 6,391 84 3,182 6,443 86 3,207 6,493 88 3,232 6,543 90 3,256 6,593 92 3,280 6,641 94 3,304 6,689 96 3,327 6,736 98 3,350 6,782

100 3,372 6,828

Tabella 6. 6 – Curve di Strain-Rate per il calcestruzzo a Trazione e Compressione adottate nella modellazione.


146

DIF (Calcestruzzo)

1,0

10,0

1E-05 0,0001 0,001 0,01 0,1 1 10 100 1000

ε'

DIF Compressione

Trazione

Grafico 6. 2 – Funzione Bilineare in scala logaritmica del DIF per il calcestruzzo.

Dal grafico precedente si può facilmente osservare come l’aumento di

resistenza a trazione sia più del doppio di quello a compressione, ed

inoltre sì evidenzia l’andamento della curva bilineare, con il cambio

d’inclinazione intorno al valore di ε = 30 sec-1.

Per convenzione il codice di calcolo LS-DYNA considera positive le

deformazioni di trazione e negative quelle di compressione [11], e così

anche le rispettive velocità di deformazione. In conclusione, si è definita

un'unica curva di Strain-Rate, assegnando per valori negativi di velocità

di deformazione l’andamento del DIF a compressione, e per valori

positivi l’andamento del DIF a trazione:


147

Grafico 6. 3 – Input in LS-Dyna della Curva di Strain-Rate a Compressione e Trazione per il calcestruzzo.

Avendo già premesso che il Material Model scelto per la simulazione

numerica permette di considerare in maniera disaccoppiata gli effetti

benefici di Strain-Rate per l’acciaio e per il calcestruzzo, si dovrà

definire la curva ε -DIF anche per l’acciaio. Per gli acciai da armatura

lenta, il bollettino n° 187 del CEB-FIP suggerisce una formulazione

analitica per il calcolo del DIF, basato sul modello di Johnson and Cook.

Questo ultimo, infatti, è tra tutti i modelli teorici, quello che più si

avvicina ai risultati dei test sperimentali condotti su questo tipo di

materiale.

Si ottengono così, per un tipico acciaio da carpenteria laminato a caldo la

seguente formulazione, valida per la sola tensione di snervamento:

⎟⎟⎠

⎞⎜⎜⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛+=

0

ln0.61εε

ysys

yd

fff


148

Dove: fyd è la tensione di snervamento dinamica;

fys è la tensione di snervamento statica;

ε è la velocità di deformazione;

0ε è la velocità di deformazione Quasi-Statica, e si

assume pari a 5·10-10 s-1;

Il DIF adoperato per l’acciaio sarà quindi dato da:

DIF (Acciaio)

fys [Mpa] 388,17

ε0' [sec-1] 5E-10

ε' [sec-1] Ln(ε/ε0)

0,0000000005 0,000 1,000 3,00E-05 11,002 1,170 3,00E-04 13,305 1,206 3,00E-03 15,607 1,241 3,00E-02 17,910 1,277

0,3 20,212 1,312 2 22,110 1,342 6 23,208 1,359 10 23,719 1,367 14 24,055 1,372 18 24,307 1,376 22 24,507 1,379 26 24,675 1,381 30 24,818 1,384 35 24,972 1,386 40 25,105 1,388 45 25,223 1,390 50 25,328 1,392 55 25,424 1,393 60 25,511 1,394 65 25,591 1,396 70 25,665 1,397 75 25,734 1,398 80 25,798 1,399 85 25,859 1,400 90 25,916 1,401 95 25,970 1,401

100 26,022 1,402 104 26,061 1,403


149

108 26,099 1,403 112 26,135 1,404 116 26,170 1,405 120 26,204 1,405 124 26,237 1,406 128 26,268 1,406 132 26,299 1,407 136 26,329 1,407 140 26,358 1,407 144 26,386 1,408 148 26,414 1,408 152 26,440 1,409 156 26,466 1,409 160 26,492 1,409 164 26,516 1,410 168 26,540 1,410 172 26,564 1,411 174 26,575 1,411

Tabella 6. 7 – DIF a trazione adoperato per L’acciaio.

DIF (Acciaio)

1,0

1,1

1,2

1,3

1,4

1,5

1,6

0 20 40 60 80 100 120 140 160

ε'

DIF Trazione

Grafico 6. 4 – Curva di Strain-Rate a trazione adoperato per l’acciaio.

6.3 – Equation of State Form 8 (Tabulated Compaction) Nel momento in cui si sceglie di adottare i Material Model Pseudo

Tensor o Concrete Damage per la modellazione del calcestruzzo, in

entrambi i casi sarà necessario per la definizione del modello, associare a

questi materiali un equazione di stato [11].


150

Nella trattazione teorica dei due Material Model, si fa espressa richiesta

dell’impiego dell’Equazione di Stato Tabulated Compaction. Questa

equazione esprime la Pressione Media p in funzione della Deformazione

Volumetrica εv [11].

La pressione così calcolata sarà impiegata dal software per la definizione

delle varie superfici di rottura, in precedenza introdotte nella descrizione

dei Legami costitutivi del calcestruzzo.

Quella che indicata con p non è altro che la parte sferica del tensore delle

tensioni, o tensione sferica, definita come il valore medio delle tensioni

normali principali σ1, σ2 e σ3.

Compito dell’utente è quello di assegnare in maniera tabulata i parametri

C, γ, T ed E, da cui dipende la pressione.

Dove:

• C: Pressione espressa in Pa;

• γ: Coefficiente numerico adimensionale;

• T: temperatura espressa in gradi Kelvin;

• E: è l’energia interna.

Invece, la deformazione volumetrica εv è definita come il logaritmo

naturale del volume relativo V (LnV = εv).

La deformazione volumetrica εv di un elemento solido, quale ad esempio

un semplice prisma di dimensioni iniziali x0, y0 e z0, definite in un

sistema di riferimento (x, y, z), si può calcolare come:

( ) ( ) ( )1.6ETCp vv εγε +=

3210

0

0

εεεε ++=−

=Δ

=V

VVV

Vv


151

Dove:

• V0 è il volume iniziale, V0 = (x0 y0 z0);

• ΔV è la variazione di volume, ΔV = Δx Δy Δz

Δx = εx xo ; Δy = εy yo ; Δz = εz zo;

• V è il volume finale o relativo, V = V0 + ΔV;

Pertanto, il volume relativo V, da cui dipende la deformazione

volumetrica εv, si esprime come:

V = (x0 y0 z0) + [(εx xo)(εy yo)(εz zo)]

Attraverso semplici passaggi matematici, e semplificando i termini

d’ordine superiore, si giunge alla seguente conclusione:

V = V0 [1 + εx + εy + εz] ⇒ V = V0 [1 + εv]

Sapendo che la deformazione volumetrica si esprime come somma delle

tre deformazioni normali.

In condizioni di simmetria assiale le due deformazioni radiali εx e εy

sono tra loro uguali, quindi:

εv = εz + 2 εx

sapendo che il coefficiente di Poisson ν, esprime il rapporto che sussiste

tra deformazione assiale e quella radiale, quest’ultima si esprime in

funzione di quell’assiale come:

εx = −ν εz ⇒ εv = εz (1- 2ν)


152

avendo assunto per convenzioni, positive le deformazioni assiali se di

contrazione, e negative le rispettive deformazioni radiali di dilatazione.

In definitiva il volume relativo V si esprime come:

V = V0 [1 + εz (1- 2ν)]

Per semplicità si assumerà un volume iniziale unitario (V0 = 1).

Nella Card dell’EOS8 si possono tabulare fino ad un massimo di 10

valori di LnV. Ciononostante, anche soli due valori sono sufficienti per

ottenere un’estrapolazione della pressione [11].

L’andamento della pressione in funzione della deformazione

volumetrica, è illustrato nel successivo grafico [11]:

Grafico 6. 5 – Legame costitutivo Pressione-Deformazione Volumetrica. Durante la fase di carico, la pressione segue l’andamento dettato

dall’equazione 6.1, quando il materiale è invece scaricato, la pressione

segue un percorso deformativo, la cui inclinazione è data dal Bulk


153

Modulus di Unloading, che è a sua volta funzione della deformazione

volumetrica.

Il Bulk Modulus è una costante di elasticità del materiale, definita in

funzione del Modulo di Young e del Rapporto di Poisson:

Il Bulk Modulus viene calcolato in funzione del Modulo di Young

dinamico Ed, allo Strain-Rate di 400 sec-1, essendo questo il massimo

valore di velocità di deformazione che si consegue nella modellazione

numerica, secondo CEB-FIP Model Code 1990 [20].

ν 0,2

ε' [sec-1] 400,00

Es [Mpa] 17636,64

Ed [Mpa] 27018,69

Ku [Pa] 1,501E+10

Tabella 6. 8 – Bulk Modulus.

Ad ogni valore di εv in maniera tabulata è associato nella Card

dell’Equazione di Stato il valore costante di K, così come sopra definito.

Scaricando il materiale la pressione si riduce fino ad annullarsi, e si

osserva un valore residuo di deformazione volumetrica. Durante la fase

di ricarico, la pressione aumenta seguendo lo stesso percorso

deformativo, fino a ricongiungersi con la curva descritta dall’equazione

6.1, nello stesso punto in cui inizia la fase di scarico, per poi continuare

ad aumentare.

( )υ213 −=

EK


154

Assumendo γ = 0, l’Equazione di Stato si riduce in una forma più

semplice:

In questo modo per definire l’Equazione di Stato sarà sufficiente

assegnare, in maniera tabulata, ad ogni valore di εv il corrispondente

valore di pressione p.

Non avendo a disposizione dati provenienti da prove sperimentali,

finalizzate alla caratterizzazione dei materiali impiegati per lo

svolgimento della prova, si è costretti a ricavare numericamente questo

legame costitutivo del calcestruzzo (εv – p).

A tal fine si è utilizzato il modello numerico Microplane Model for

Brittle-Plastic Material, proposto da Bažant Z P, Professore presso il

Dipartimento d’Ingegneria Civile, della Northwestern University,

Evanston [21].

Il modello formulato da Bažant e Prat nell’1988, è stato calibrato e

verificato su base sperimentale, comparando i risultati numerici con i

risultati di diverse prove di compressione tri-assiali.

Le superfici di rottura sono definite dalle stesse funzioni utilizzate nel

Microplane Model di Bažant and Prat 1988, per la relazione tensio-

deformativa volumetrica e deviatorica. A noi interessa definire la

relazione di tipo volumetrica, che sussiste tra la tensione e la

deformazione (σv, εv).

Nella teoria classica della plasticità, una volta definita la superficie di

snervamento (Yeld Surface), la deformazione plastica inizia a crescere

nel momento in cui la tensione raggiunge detta superficie. Invece, nel

modello proposto le regioni elastiche e plastiche non sono separate, si

( )vCp ε=


155

assume quindi che la transizione continua della superficie di rottura

(Failure Surface) tra le due regioni, avvenga mediante una deformazione

plastica equivalente [21].

La relazione tra tensione e deformazione volumetrica si può definire

mediante la seguente relazione:

Dove:

• σv è la tensione volumetrica, che in condizione di carico

idrostatiche sarà equivalente alla sola σ3. Essendo σv = p = (σ1 +

σ2 + σ3)/3, in condizioni idrostatiche le tre tensioni principali sono

uguali tra loro, quindi p = σv = σ3.

• εpl è la deformazione plastica. Sapendo che la deformazione

volumetrica è somma delle due aliquote, plastica ed elastica, sarà

εpl = εv – εel. La deformazione elastica εel, viene calcolata come

rapporto tra la pressione p, incognita del problema, e il Bulk

Modulus, in precedenza calcolato.

• f’c è la resistenza a compressione del calcestruzzo non

confinato. Il vantaggio nell’utilizzare questo metodo risiede

proprio nella possibilità d’adimensionalizzare la tensione σv

rispetto alle caratteristiche meccaniche del calcestruzzo, espresse

mediante il parametro f’c.

• a, b, s e t sono coefficienti numerici adimensionali, definiti

mediante prove sperimentali. Nel momento in cui si definisce la

superficie di rottura relativa ad uno stato tensionale idrostatico, tali

coefficienti sono determinati mediante i risultati di prove

pl

tpl

splc

c

v

baE

fε

εευ

σ⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛+

−=

−

121'


156

idrostatiche tri-assiali. Si è scelto di utilizzare i coefficienti definiti

da Green and Swanson:

a 0,000025b 2,5 s 0,8 t 1,15

Tabella 6. 9 – Coefficienti numerici di Green and Swanson.

All’equazione così descritta compete la seguente curva tensio-

deformativa volumetrica [21]:

Figura 6. 5 – Legame costitutivo Tensione-Deformazione Volumetrica.

Fissato un valore di deformazione volumetrica εv, si calcola quella

plastica εpl come:

Dato che la pressione p, relativa alla deformazione volumetrica fissata, è

l’incognita, sarà necessario ipotizzare un valore iniziale di pressione

diverso da zero. Il valore di pressione p, per il quale si annulla la

seguente differenza, sarà la soluzione del problema:

Kpdove elelvpl =−= εεεε


157

Alla deformazione volumetrica εv corrisponde il valore di pressione p,

per il quale la precedente differenza si annulla.

Iterando questo procedimento per sette valori differenti di εv, si

ottengono i seguenti risultati numerici:

Legame p - εv

εv p [Mpa]

0,000000 0,000E+00

0,007009 4,430E+01

0,028545 1,131E+02

0,049701 2,241E+02

0,060903 2,959E+02

0,070039 3,595E+02

0,134452 896,27

Tabella 6. 10 – Risultati Numerici del legame costitutivo calcolato.

0121' =

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛+

−−

−

pl

tpl

splc

c baE

fp ε

εευ


158

Legame P - εv

0,0

100,0

200,0

300,0

400,0

500,0

600,0

700,0

800,0

900,0

1000,0

0,000 0,020 0,040 0,060 0,080 0,100 0,120 0,140 0,160

−εv

-P [M

pa

Grafico 6. 6 – Legame Tensio-Deformativo Volumetrico.

In funzione della deformazione volumetrica, si calcola il volume relativo

V, cosi come descritto in precedenza, se ne calcola il logaritmo naturale

LnV, e si associa a lui il corrispondente valore di pressione p, espressa in

Pa, così come richiesto nella Card dell’Equation of State [11].

V lnV P [Pa] 1,0000000 0,000000 0,000000000 0,9929907 -0,007034 44300397,98 0,9714553 -0,028960 113129376,66 0,9502986 -0,050979 224085069,19 0,9390965 -0,062837 295858999,46 0,9299607 -0,072613 359539685,41 0,8655484 -0,144392 896273020,00

Tabella 6. 11 – Valori Tabulati dell’equazione di stato.


159

Equazione di Stato 8

0,E+00

1,E+08

2,E+08

3,E+08

4,E+08

5,E+08

6,E+08

7,E+08

8,E+08

9,E+08

1,E+09

-0,200 -0,150 -0,100 -0,050 0,000

lnV

P[Pa

]

EOS8

Grafico 6. 7 – Equazione di Stato assegnata come input in LS-Dyna, per la modellazione della piastra in c.a.

Per una Card dell’Equation of State che si va così a comporre:

*EOS_TABULATED_COMPACTION $# eosid gama e0 vo 1 0.000 0.000 1.000000 $# ev1 ev2 ev3 ev4 ev5 0.000 -0.00703000 -0.02896000 -0.05098000 -0.06284000 $# ev6 ev7 ev8 ev9 ev10 -0.07261000 -0.14439000 $# c1 c2 c3 c4 c5 0.000 4.4300396e+007 1.1312938e+008 2.2408507e+008 2.9585901e+008 $# c6 c7 c8 c9 c10 3.5953968e+008 8.9627302e+008 $# t1 t2 t3 t4 t5 0.000 0.000 0.000 0.000 0.000 $# t6 t7 t8 t9 t10 0.000 0.000 0.000 0.000 0.000 $# k1 k2 k3 k4 k5 1.5010385e+010 1.5010385e+010 1.5010385e+010 1.5010385e+010 1.5010385e+010 $# k6 k7 k8 k9 k10 1.5010385e+010 1.5010385e+010

Tabella 6. 12 – Input per LS-Dyna dell’equazione di stato Tabulated Compaction.


160

Le Card così definite per il Concrete Damage, il DIF dei materiali

impiegati, e l’Equation of State, si ripeteranno identicamente nei diversi

modelli realizzati, sia in FEM sia in SPH.

CCAAPPIITTOOLLOO 77 –– RRIISSUULLTTAATTII DDEELLLLEE SSIIMMUULLAAZZIIOONNII NNUUMMEERRIICCHHEE

161

Capitolo 7 Risultati delle Simulazioni Numeriche

In conclusione si è svolto uno studio parametrico tra i modelli realizzati

con i due metodi descritti, FEM e SPH. Tale studio è reso possibile dalle

evidenze sperimentali a disposizione. Operando, infatti, un confronto tra

i risultati numerici forniti dalla modellazione, e quelli risultanti dalle

prove sperimentali realmente condotte, sarà possibile valutare quale dei

due metodi fornisce risultati numerici quanto più prossimi a quelli

sperimentali.

Per ogni prova sperimentale si avranno due modelli, uno realizzato in

FEM e l’altro in SPH, entrambi costruiti secondo le modalità descritte

nel corso della tesi.

Per ogni elemento in cui è stata suddivisa la piastra, che sia questo un

Element Solid o un Element SPH, si potranno leggere una serie

d’informazioni al variare del tempo.

Queste ad esempio possono essere le Tensioni, gli Spostamenti e le

Velocità del singolo elemento, misurate lungo gli assi del sistema di

riferimento cartesiano, o le rispettive risultanti di spostamento e velocità.

Il codice di calcolo consente anche di osservare dei risultati inerenti al

comportamento complessivo dell’intero elemento strutturale. A riguardo

si possono misurare l’Energia Cinetica, Interna e Totale, lo Spostamento,

la Velocità e l’Accelerazione del corpo rigido, nelle direzioni del sistema

di riferimento cartesiano, e le rispettive risultanti.

Per i singoli nodi si possono avere delle misure riguardanti le

Coordinate, gli Spostamenti, le Velocità, le Accelerazioni, e gli

Spostamenti Relativi.


162

7.1 – Prova Sperimentale I e II Per queste prime due condizioni di carico, le prove sperimentali non

hanno mostrato risultati apprezzabili o quantificabili numericamente.

In questi due casi sarà quindi sufficiente, per verificare la veridicità della

modellazione, osservare che non si verifichino spostamenti,

deformazioni o fessure di ampiezze apprezzabili sui modelli realizzati.

In entrambi i casi si è stabilito un tempo di analisi di 0,1 sec, suddiviso in

un numero complessivo di 52 Step. Ciò significa che il software

restituirà un output del calcolo effettuato ogni 0,002 sec.

Per i modelli realizzati in SPH, tenendo presente che questo metodo

impiega un tempo maggiore per convergere numericamente, ad

differenza del FEM, in questo caso si sceglierà un tempo di analisi di 0,2

sec, per un totale 102 Step, uno ogni 0,002 sec.

Ai fini della modellazione, e alla luce delle evidenze mostrate dalle

prove sperimentali, si è scelto di evidenziare la Deformazione Plastica

Effettiva (Effective Plastic Strain), come risultato conseguito dalla

modellazione delle prime due condizioni di carico. Non avendo la

possibilità di effettuare nessun confronto numerico con i risultati della

prova realizzata in laboratorio, in questo caso si preferisce mettere in

evidenza le differenze che sussistono tra gli stessi modelli realizzati in

FEM e in SPH. Abbiamo quindi confrontato allo stesso istante di tempo

t = 0,1 sec, l’andamento della deformazione plastica effettiva,

all’intradosso e all’estradosso della piastra.

In maniera sintetica si riportano di seguito i risultati conseguiti per i due

modelli:


163

• Modello FEM, suddivisione della piastra in 11520 Element_Solid,

Prova 1:

Figura 7. 1 – Effective Plastic Strain misurato all’estradosso della piastra, Prova 1.

Figura 7. 2 – Effective Plastic Strain misurato all’intradosso della piastra, Prova 1.


164

Grafico 7. 1 - Spostamento Residuo, misurato nel centro della piastra all’estradosso.


Prova 2:

Figura 7. 3 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 2.


165

Figura 7. 4 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 2.

Grafico 7. 2 - Spostamento Residuo, misurato nel centro della piastra all’estradosso.


166

• Modello SPH, suddivisione della piastra in 11520 Element_SPH,

Prova 1:



167


Grafico 7. 3 - Spostamento Residuo, misurato nel centro della piastra all’intradosso.


168


Prova 2:




169

Grafico 7. 4 - Spostamento Residuo, misurato nel centro della piastra all’intradosso.

7.2 – Prova Sperimentale III La terza condizione di carico, è quella che consente di osservare

sperimentalmente i risultati più interessanti. Infatti, in questo caso si

hanno misure effettuate al termine della prova, relative allo spostamento

residuo valutato nel centro della piastra, e sulla mezzeria del bordo

esterno non vincolato. Inoltre, si ha anche la misura della massima

ampiezza della fessura, che si manifesta all’intradosso della piastra, per

effetto dei carichi dinamici a cui è stata sottoposta.

Si ha quindi la possibilità di operare un confronto tra i risultati conseguiti

numericamente, mediante la modellazione, e quelli derivanti dalla prova

sperimentale, eseguita in laboratorio.


170

In questo modo si può verificare quale tra i due modelli, realizzati con i

metodi FEM e SPH, meglio approssima i dati sperimentali a

disposizione.

A tal fine, si misurerà lungo l’asse Z lo spostamento che la piastra

subisce nel centro e sul bordo, ed infine l’ampiezza della fessura che si

presenta centralmente all’intradosso dell’elemento strutturale.

La prova sperimentale ha evidenziato i seguenti risultati:

Figura 7. 9 – Immagine della piastra in seguito alla terza prova.

Wc [mm] Zr,c [cm] Zr,b [cm]

3 2,1 1,1

Tabella 7. 1 - Risultati numerici della terza Prova Sperimentale realizzata.

Dove con Wc s’indica l’ampiezza della fessura misurata, con Zr,c, Zr,b, si

indicano rispettivamente gli spostamenti misurati lungo l’asse Z in

prossimità del centro e della mezzeria del bordo esterno non vincolato

della piastra.


171


Prova 3:

Così come per le prime due prove, anche in questo caso si è scelto un

tempo d’analisi di appena 0,1 sec. Si può, infatti, osservare che il metodo

numerico utilizzato, converge asintoticamente ad un valore finito dello

spostamento in tempi brevissimi. Si riporta di seguito il risultato

conseguito mediante la modellazione all’istante di tempo t = 0,1 sec.

I nodi rispetto ai quali si sono misurati lo spostamento lungo l’asse Z,

sono rappresentati nelle seguenti immagini, nelle quali è possibile

osservare anche la deformazione che la piastra subisce, per effetto del

carico a cui è sottoposta all’istante di tempo t = 0,1 sec.

Per apprezzare anche la riduzione di spessore della piastra, si

misureranno gli spostamenti residui sia per elementi superiori sia

inferiori.

Figura 7. 10 - Nodi centrali della Mesh.


172

Figura 7. 11 - Nodi sul Bordo della Mesh.

Si mostrano in seguito le deformazioni plastiche effettive all’estradosso e

all’intradosso della piastra:


173




174

Graficamente si rappresenta di seguito l’andamento degli spostamenti

nodali al variare del tempo e in precedenza individuati sulla piastra.

Grafico 7. 5 – Spostamenti misurati nel centro della piastra all’intradosso e

all’estradosso.

Grafico 7. 6 – Spostamenti misurati sul bordo della piastra all’intradosso

all’estradosso.


175

Grafico 7. 7 – Ampiezza della fessura misurata all’intradosso della piastra in

posizione centrale.

Per misurare l’ampiezza della fessura, si calcola la distanza che sussiste

tra gli spostamenti di due elementi all’intradosso della piastra, adiacenti

tra loro. Tale differenza si misura lungo l’asse X, essendo questa la

direzione lungo la quale si ha l’apertura della lesione.

A tale scopo sono stati individuati gli elementi 1176 e 1177

all’intradosso della piastra e in posizione centrale.

In maniera sintetica si riportano di seguito i risultati numerici conseguiti:

Cond. Vincolo b1 b2 b3 Damping t [sec] Incastrata 0.1 0.1 0.1 5% 0,1

Tabella 7. 2 - Parametri che caratterizzano il modello in LS-Dyna Version 970.

Wc [mm] 0,82

Tabella 7. 3 – Ampiezza numerica della Fessura misurata nel centro della piastra all'intradosso.


176

Nodi centrali Zr,c [cm]

Estradosso 2,15 Intradosso 1,52

Nodi di bordo Zr,b [cm]


Tabella 7. 4 - Spostamenti residui, centrali e di bordo, misurati all’estradosso e all’intradosso della piastra.

Dove i nodi 13206 e 1201, sono quelli centrali, mentre i nodi 12030 e

2377, sono quelli di bordo, rispettivamente superiori ed inferiori.

Calcolando la differenza tra i due si può verificare di quanto la piastra si

riduce di spessore.


Prova 3:

Il metodo SPH, a differenza del FEM, presenta una maggiore instabilità.

D’altra parte questo suo limite è già stato messo in evidenza nel Capitolo

5. Ne consegue che il tempo necessario d’analisi affinché si possono

osservare dei risultati numerici asintotici ad un valore costante, deve

essere maggiore di quello imposto nel caso precedente. A tal fine si è

scelto di fare terminare l’analisi ad un istante di tempo t = 1,2 sec.

Si riportano di seguito i risultati conseguiti mediante la modellazione

all’istante di tempo t = 1,2 sec.

I nodi rispetto ai quali si è misurato lo spostamento lungo l’asse Z, sono

rappresentati nelle seguenti immagini, nelle quali è possibile osservare

anche la deformazione che la piastra subisce, per effetto del carico cui è

sottoposta all’istante di tempo t = 1,2 sec.


177

A differenza del modello realizzato in FEM, non si ha la possibilità di

individuare sia centralmente sia sul bordo un unico nodo, che sia questo

sull’estradosso o l’intradosso della piastra. In questo caso abbiamo

quindi scelto come elementi nodali di riferimento, quelli riportati nella

seguente tabella: Particelle

Centrali Bordo 10344 10345 9240 9241 5736 5737 4632 4633 1128 1129 24 25 10392 10393 5784 5785 1176 1177

Tabella 7. 5 - Elementi nodali centrali e di bordo sulla piastra discretizzata in SPH.

Di seguito si riportano le immagini della piastra all’istante di tempo

t = 1,2 sec, sulla quale sono evidenziate alcune delle particelle di centro

e di bordo, oltre alla deformazione della piastra stessa.

Figura 7. 14 – Elementi Nodali Centrali e di Bordo.


178

Figura 7. 15 – Elementi Nodali di Bordo.

Si mostrano in seguito le deformazioni plastiche effettive all’estradosso e

all’intradosso della piastra:



179


Per quanto concerne gli spostamenti misurati al centro e sul bordo della

piastra, e la conseguente misura dell’ampiezza massima della fessura, si

rimanda ai successivi grafici:

Grafico 7. 8 – Spostamenti misurati nel centro della piastra all’intradosso e

all’estradosso.


180

Grafico 7. 9 – Spostamenti misurati sul bordo non vincolato della piastra,

all’intradosso e all’estradosso.

Grafico 7. 10 – Ampiezza della fessura misurata nel centro della piastra

all’intradosso.

In maniera sintetica si riportano di seguito i risultati numerici conseguiti:

Cond. Vincolo b1 b2 b3 Damping t [sec] Incastrata 0,1 0,1 0,1 5% 1,2

Tabella 7. 6 - Parametri che caratterizzano il modello in LS-Dyna Version 970.


181

Wc [mm] 13

Tabella 7. 7 – Ampiezza numerica della Fessura misurata nel centro della piastra.

Particelle centrali Zr,c [cm]


Nodi di bordo Zr,b [cm]


Tabella 7. 8 - Spostamenti residui, centrali e di bordo, misurati all’estradosso e all’intradosso della piastra.

Dove le particelle 10344 e 1128, sono quelle centrali, mentre 9240 e 24,

sono quelli di bordo, rispettivamente superiori ed inferiori. Calcolando la

differenza tra i due si può verificare la riduzione di spessore che la

piastra subisce per effetto dei carichi cui è sottoposta.

7.3 – Prova Sperimentale IV Per la quarta prova sperimentale si sono realizzate delle condizioni di

carico tali da produrre la rottura dell’elemento strutturale. In questo caso

non si hanno quindi termini di confronto, tra l’evidenza sperimentale e la

modellazione numerica, se non la rottura effettiva dell’elemento.

Tuttavia, quest’ultima prova consente di apprezzare i vantaggi

nell’utilizzare un metodo di discretizzazione particellare, quale l’SPH

impiegato nella modellazione, a fronte degli insuccessi mostrati dalla

modellazione effettuata in FEM. In questo caso infatti, emerge uno dei

principali limiti del Finite Element Method, dato che le condizioni di

carico sono tali da indurre una deformazione eccessiva della Mesh,

realizzata per la discretizzazione dell’elemento continuo, che non


182

consente al metodo numerico di convergere. Per tale motivo non si ha la

possibilità di mostrare alcun risultato apprezzabile della modellazione

effettuata, in quanto il codice di calcolo impiegato, data l’incapacità di

gestire una distorsione eccessiva della maglia, non fornisce nessun

risultato d’output del calcolo realizzato, in quanto la computazione

numerica termina all’istante.

Lo Smoothed Particle Hydrodynamics, consente di superare questo

limite, come già annunciato nel capitolo dedicato alla descrizione dei

due metodi di discretizzazione. Uno dei principali vantaggi di tale

metodo risiede proprio nella capacità di trattare facilmente grandi

deformazioni, essendo le connessioni tra i nodi generate come una parte

della computazione stessa, che può subire variazioni nel tempo.

Di conseguenza si riporteranno in questo caso i soli risultati relativi alla

modellazione effettuata in SPH.

Nelle immagini successive si evidenzia l’andamento della deformazione

plastica effettiva, per due diversi istanti di tempo, per la piastra vista

dall’alto e lateralmente.

L’istante di tempo a cui si riferiscono le immagini successive è t = 0,02

sec. prima ancora che si consegua la condizione di rottura.


183

Figura 7. 18 – Effective Plastic Strain all’estradosso, all’istante t=0,02 sec.

Figura 7. 19 – Effective Plastic Strain nel piano (y, z), all’istante t = 0,02 sec.

Le immagini successive invece si riferiscono all’istante successivo

t = 0,04 sec. quando la piastra inizia a presentare segni di rottura.


184

Figura 7. 20 – Effective Plastic Strain all’estradosso all’istante t = 0,04 sec.

Figura 7. 21 – Effective Plastic Strain nel piano (y, z), all’istante t = 0,04 sec.

Si vuole mostrare in seguito, il risultato della modellazione che si

consegue all’istante di tempo t = 0,6 sec. in quest’immagine si possono

osservare le particelle che sono rimaste vincolate, e che quindi non

hanno subito spostamenti, e delle frange di calcestruzzo che restano a

loro ancorate.


185

Figura 7. 22 – Rottura della piastra all’istante t = 0,6 sec.

CCOONNCCLLUUSSIIOONNII

186

CONCLUSIONI

Alla luce dei risultati conseguiti mediante la modellazione numerica

delle quattro prove sperimentali, si possono trarre le successive

conclusioni.

Per quanto concerne lo studio parametrico condotto tra i due metodi

d’analisi adottati, FEM e SPH, per la modellazione d’elementi strutturali

sottoposti a carichi da esplosione, si può affermare che fin quando le

condizioni di carico inducono velocità di deformazioni sufficientemente

basse, in particolar modo inferiori a 300 sec-1, a fronte di quelle che si

possono ottenere per carichi da esplosione, come si evince dai risultati

numerici in precedenza esposti, entrambi i metodi ben si prestano ai fini

della modellazione. Tale dato emerge in particolar modo dai risultati che

si ottengono per la modellazione delle prime due prove sperimentali,

nelle quali si possono osservare risultati simili tra i due metodi, non

essendoci, infatti, spostamenti apprezzabili della piastra in entrambi i

casi.

All’aumentare delle condizioni di carico, e facendo un esplicito

riferimento alla terza prova simulata numericamente, si può apprezzare il

vantaggio nell’utilizzare il Finite Element Method, piuttosto che il più

innovativo Smoothed Particle Hydrodynamics. In virtù dei confronti

operati tra i risultati numerici e quelli messi in evidenza dalla prova

sperimentale, si osserva, infatti, che per la stessa condizione di carico, di

vincolo, a parità di Strain-Rate, di legami costitutivi, e facendo variare il

solo metodo d’analisi, il FEM restituisce risultati numerici più prossimi a

quelli sperimentali.


187

In questo caso si osserva la variabilità di risultati che si possono ottenere

al variare dell’infittimento della Mesh realizzata, e dei valori numerici

associati ai tre parametri b1, b2 e b3, già descritti nel Capitolo 6.

Per ognuno dei due metodi, FEM e SPH, fissata la dimensione della

Mesh realizzata, e facendo variare i soli valori numerici associati ai tre

parametri sopra indicati, si ottengono risultati in termini di spostamenti

residui differenti. Lo stesso risultato si può osservare mantenendo

inalterati i valori numerici dei tre parametri sopra indicati, e facendo

variare l’infittimento della Mesh in cui è stata suddivisa la piastra.

Per i metodi FEM e SPH, dagli spostamenti residui misurati

all’intradosso della piastra, rispettivamente sui nodi centrali e di bordo

1201 e 2377, e sulle particelle centrali e di bordo 1128 e 24, si riportano i

seguenti risultati:

F E M

Mesh b1 b2 b3 t [sec] Zr,c [cm] Zr,b [cm] Wc [mm] 0,1 0,1 0,1 0,1 1,2 0,6 0,45 6912 Element Solid 1 1 1 0,1 1,06 0,47 0,004

0,1 0,1 0,1 0,1 1,47 0,74 0,82 11520 Element Solid 1 1 1 0,1 0,99 0,45 0,36

Tabella 8. 1 – Variabilità dei risultati per il metodo FEM, al variare delle dimensioni della Mesh e dei parametri b1, b2 e b3.

S P H

Mesh b1 b2 b3 t [sec] Zr,c [cm] Zr,b [cm] Wc [mm] 0,1 0,1 0,1 1,81 5,8 2,8 1,16 6912 Element SPH 1 1 1 2,7 10 4 4,2

0,1 0,1 0,1 1,2 11,8 10,42 13 11520 Element SPH 1 1 1 1,94 9 3,4 9,9

Tabella 8. 2 - Variabilità dei risultati per il metodo FEM, al variare delle dimensioni della Mesh e dei parametri b1, b2 e b3.


188

Per il metodo SPH, data la sua maggiore instabilità numerica e la

necessità quindi di tempi d’analisi maggiore affinché si possano

osservare dei risultati finiti, avremo dei tempi d’analisi diversi rispetto a

quelli adottati per il metodo FEM. Inoltre, nelle precedenti tabelle si

evidenziano in grassetto i risultati ottenuti dalla modellazione.

Da queste tabelle si evince che il metodo FEM rispetto all’SPH, è meno

sensibile alla variazione sia dell’infittimento della Mesh, sia al variare

dei parametri numerici b1, b2 e b3.

Da quanto emerso si può quindi affermare che il metodo SPH, per questo

caso specifico, necessità di una calibrazione su base sperimentale di tali

parametri, i quali a loro volta come già evidenziato nel Capitolo 6,

assumono valori diversi al variare della dimensione della stessa Mesh.

Ai fini della modellazione, e per questo caso particolare, è quindi

preferibile impiegare il metodo FEM, data la sua minor sensibilità nei

confronti di tali fenomeni.

Dai risultati di questa stessa prova, si può altresì mettere in luce

l’attinenza del legame costitutivo Concrete Damage, proposto da

Malvar, ai fini della modellazione del calcestruzzo. Tenendo presente

che migliori risultati numerici si potrebbero conseguire, definendo i

parametri richiesti dall’autore, attraverso lo svolgimento di prove

sperimentali preposte alla caratterizzazione del calcestruzzo, così

com’esposto nel Capitolo 6.

Dai risultati della quarta prova modellata, emerge invece l’impossibilità

di impiegare il metodo d’analisi FEM, data l’eccessiva deformazione

indotta nella Mesh realizzata, e le conseguenti incapacità numeriche del

metodo di gestire tale distorsione.

In questo caso è quindi preferibile impiegare il metodo SPH, in quanto le

singole particelle non essendo legate tra loro da costrizioni di maglia,


189

non avranno problemi ad allontanarsi le une dalle altre simulando la

rottura dell’elemento, sotto l’azione dei carichi a cui sono sottoposte.

Un successivo risultato che questa prova mostra, è l’incremento di

resistenza del calcestruzzo per gli effetti di Strain-Rate. Si può, infatti,

osservare come il valore di tensione di trazione alla quale si ha la rottura

dell’elemento, sia maggiore di quell’assegnata in condizione di carico

quasi statica.

Sulla base dei risultati conseguiti, e delle conclusioni che da loro si

traggono, si potrà in seguito indirizzare la ricerca sui seguenti aspetti:

• Le stesse prove sperimentali andrebbero modellate utilizzando

altri metodi d’analisi, quali ad esempio quelli Euleriani,

l’Arbitrarian Lagrangian Eulerian (ALE), o altri metodi di tipo

Meshfree.

• Si potrebbero realizzare simulazioni numeriche di prove

sperimentali condotte su piastre interamente d’acciaio, o su

elementi strutturali iperstatici, con il fine di valutare gli effetti di

Strain-Rate, che potrebbero in questi casi essere più evidenti.

• Altri studi potranno essere condotti sulla base dei risultati

numerici conseguiti dalla simulazione della quarta prova,

provando a modellare in questo caso interventi di rinforzo con

materiali compostiti, messi tra loro a confronto, in modo tale da

evitare che si possa conseguire la condizione di rottura

dell’elemento, come già mostrato dalla modellazione e

dall’evidenza della prova sperimentale stessa.


190

• Si potranno calibrare mediante lo svolgimento di prove

sperimentali specifiche, i parameri numerici richiesti dal modello

Concrete Damage, per la caratterizzazione del calcestruzzo, a

differenza di quanto fatto nell’elaborato di tesi, in cui si è

pervenuti numericamente alla definizione di tali parametri.

IINNDDIICCII DDEELLLLEE FFIIGGUURREE,, DDEEII GGRRAAFFIICCII EE DDEELLLLEE TTAABBEELLLLEE

191

Indice delle Figure Figura 1. 1 – Andamento qualitativo di una Curva di Pressione prodotta da un’esplosione. ............................................................................................................ 21 Figura 2. 1 – Caratteristiche Geometriche della piastra, e dell’armatura metallica. 38 Figura 2. 2 – Setup di Prova. .................................................................................... 38 Figura 2. 3 – Immagine della piastra in seguito alla terza Prova.............................. 41 Figura 2. 4 – Immagini della piastra in seguito alla quarta Prova. ........................... 42 Figura 3. 1 – Curva di Pressione semplificata, usata dagli autori per la modellazione della prova. 44 Figura 3. 2 – Discretizzazione della superficie di carico.......................................... 46 Figura 3. 3 – Distanza effettiva dell’i-esimo elemento dal centro d’esplosione. ..... 49 Figura 4. 1 – Instabilità d trazione delle particelle. 92 Figura 4. 2 - Stress Point Method in uno spazio bi-dimensionale. ........................... 94

Figura 5. 1 - Modellazione di una Turbina con LS-Dyna. 98

Figura 5. 2 - Sezione della Piastra sottoposta a Carico distribuito. ........................ 109

Figura 6. 1 – Superficie di rottura massima e residua. 118

Figura 6. 2 – Card Material Model_016. ................................................................ 120 Figura 6. 3 – Andamento della superficie di rottura corrente Δσ, in funzione della tensione sferica p. .................................................................................................... 121 Figura 6. 4 – Superfici di Rottura. .......................................................................... 123 Figura 6. 5 – Legame costitutivo Tensione-Deformazione Volumetrica. .............. 156

Figura 7. 1 – Effective Plastic Strain misurato all’estradosso della piastra, Prova 1. 163

Figura 7. 2 – Effective Plastic Strain misurato all’intradosso della piastra, Prova 1.................................................................................................................................. 163 Figura 7. 3 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 2.................................................................................................................................. 164 Figura 7. 4 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 2.................................................................................................................................. 165 Figura 7. 5 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 1.................................................................................................................................. 166 Figura 7. 6 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 1.................................................................................................................................. 167 Figura 7. 7 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 2.................................................................................................................................. 168 Figura 7. 8 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 2.................................................................................................................................. 168 Figura 7. 9 – Immagine della piastra in seguito alla terza prova. ........................... 170 Figura 7. 10 - Nodi centrali della Mesh.................................................................. 171 Figura 7. 11 - Nodi sul Bordo della Mesh. ............................................................. 172 Figura 7. 12 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 3.................................................................................................................................. 173


192

Figura 7. 13 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 3.................................................................................................................................. 173 Figura 7. 14 – Elementi Nodali Centrali e di Bordo............................................... 177 Figura 7. 15 – Elementi Nodali di Bordo. .............................................................. 178 Figura 7. 16 - Effective Plastic Strain misurato all’estradosso della piastra, Prova 3.................................................................................................................................. 178 Figura 7. 17 - Effective Plastic Strain misurato all’intradosso della piastra, Prova 3.................................................................................................................................. 179 Figura 7. 18 – Effective Plastic Strain all’estradosso, all’istante t=0,02 sec.......... 183 Figura 7. 19 – Effective Plastic Strain nel piano (y, z), all’istante t = 0,02 sec...... 183 Figura 7. 20 – Effective Plastic Strain all’estradosso all’istante t = 0,04 sec......... 184 Figura 7. 21 – Effective Plastic Strain nel piano (y, z), all’istante t = 0,04 sec...... 184 Figura 7. 22 – Rottura della piastra all’istante t = 0,6 sec. ..................................... 185

Indice dei Grafici Grafico 3. 1 - Estrapolazione numerica del coefficiente sperimentale α.................. 51 Grafico 5. 1 – Andamento della Curva di Pressione, calcolata per l’area di carico A-1 e relativa alla prima prova modellata. 107 Grafico 6. 1 – Funzione di Danno, calcolata per b1 = b2 = b3 = 0,1. 136 Grafico 6. 2 – Funzione Bilineare in scala logaritmica del DIF per il calcestruzzo.................................................................................................................................. 146 Grafico 6. 3 – Input in LS-Dyna della Curva di Strain-Rate a Compressione e Trazione per il calcestruzzo. .................................................................................... 147 Grafico 6. 4 – Curva di Strain-Rate a trazione adoperato per l’acciaio. ................. 149 Grafico 6. 5 – Legame costitutivo Pressione-Deformazione Volumetrica. ............ 152 Grafico 6. 6 – Legame Tensio-Deformativo Volumetrico...................................... 158 Grafico 6. 7 – Equazione di Stato assegnata come input in LS-Dyna, per la modellazione della piastra in c.a.............................................................................. 159 Grafico 7. 1 - Spostamento Residuo, misurato nel centro della piastra all’estradosso. 164 Grafico 7. 2 - Spostamento Residuo, misurato nel centro della piastra all’estradosso.................................................................................................................................. 165 Grafico 7. 3 - Spostamento Residuo, misurato nel centro della piastra all’intradosso.................................................................................................................................. 167 Grafico 7. 4 - Spostamento Residuo, misurato nel centro della piastra all’intradosso.................................................................................................................................. 169 Grafico 7. 5 – Spostamenti misurati nel centro della piastra all’intradosso e all’estradosso. .......................................................................................................... 174 Grafico 7. 6 – Spostamenti misurati sul bordo della piastra all’intradosso all’estradosso. .......................................................................................................... 174 Grafico 7. 7 – Ampiezza della fessura misurata all’intradosso della piastra in posizione centrale. ................................................................................................... 175 Grafico 7. 8 – Spostamenti misurati nel centro della piastra all’intradosso e all’estradosso. .......................................................................................................... 179 Grafico 7. 9 – Spostamenti misurati sul bordo non vincolato della piastra, all’intradosso e all’estradosso.................................................................................. 180


193

Grafico 7. 10 – Ampiezza della fessura misurata nel centro della piastra all’intradosso............................................................................................................ 180

Indice delle Tabelle Tabella 1. 1 – Parametri che caratterizzano fisicamente alcuni materiali esplosivi . 35 Tabella 2. 1 – Proprietà elastiche dei materiali impiegati per la costruzione della piastra, allo Strain-Rate di 100 sec-1. 39 Tabella 2. 2 – Proprietà dinamiche del calcestruzzo allo Strain-Rate di 100 sec-1. .. 39 Tabella 2. 3 – Tensione di Snervamento dell’acciaio, per due diversi valori di Strain-Rate. ........................................................................................................................... 39 Tabella 2. 4 – Proprietà del calcestruzzo, per velocità di deformazione quasi statica.................................................................................................................................... 40 Tabella 2. 5 – Condizioni di carico ed Evidenze Sperimentali ad esse associate. .... 41 Tabella 3. 1 – Calore specifico degli esplosivi e Coefficiente d’Equivalenza. 48 Tabella 3. 2 - Distanza Effettiva Ri, Carico Equivalente WTNT, Distanza Ridotta Zi.................................................................................................................................... 52 Tabella 3. 3 - Picco di Pressione Δp(Z), calcolato secondo le varie formulazioni, e per le singole prove sperimentali. .............................................................................. 53 Tabella 3. 4 - Picco di Pressione medio, Coefficiente α........................................... 53 Tabella 3. 5 - Durata della fase positiva t+, e costante numerica B. ......................... 53 Tabella 3. 6 - Curve di Pressione P(Z, t), calcolate per le quattro prove sperimentali.................................................................................................................................... 54

Tabella 5. 1 – Unità di Misure adottate dal LS-Dyna Version 970. 101

Tabella 5. 2 – Element Solid................................................................................... 102 Tabella 5. 3 – Part. .................................................................................................. 103 Tabella 5. 4 – Section Solid. ................................................................................... 103 Tabella 5. 5 – Set_Node_List Vincoli..................................................................... 105 Tabella 5. 6 – Boundary Condiction. ...................................................................... 105 Tabella 5. 7 – Curva di Pressione Tabulata. ........................................................... 107 Tabella 5. 8 - Set_Segment. .................................................................................... 108 Tabella 5. 9 – Load Set_Segment. .......................................................................... 109 Tabella 5. 10 – Set_Part_List.................................................................................. 110 Tabella 5. 11 – Carico da Peso Proprio................................................................... 110 Tabella 5. 12 – Element SPH.................................................................................. 112 Tabella 5. 13 – Section SPH. .................................................................................. 112 Tabella 5. 14 – Set_Node_List Carico. ................................................................... 114 Tabella 5. 15 – Load_Node_Set.............................................................................. 115 Tabella 5. 16 – Damping......................................................................................... 116

Tabella 6. 1 – parametri che definiscono la Card del Material Model. 120

Tabella 6. 2 – Parametri numerici che definiscono le superfici di rottura d’input per LS-Dyna................................................................................................................... 131 Tabella 6. 3 – Input per LS-Dyna del material Concrete Damage, definiti per la modellazione della piastra. ...................................................................................... 136 Tabella 6. 4 – Velocità di deformazione quasi statica. ........................................... 144


194

Tabella 6. 5 – Caratteristiche statiche del calcestruzzo e parametri numerici utilizzati per il calcolo delle curve di Strain-Rate................................................................... 144 Tabella 6. 6 – Curve di Strain-Rate per il calcestruzzo a Trazione e Compressione adottate nella modellazione. .................................................................................... 145 Tabella 6. 7 – DIF a trazione adoperato per L’acciaio. .......................................... 149 Tabella 6. 8 – Bulk Modulus................................................................................... 153 Tabella 6. 9 – Coefficienti numerici di Green and Swanson. ................................. 156 Tabella 6. 10 – Risultati Numerici del legame costitutivo calcolato. ..................... 157 Tabella 6. 11 – Valori Tabulati dell’equazione di stato.......................................... 158 Tabella 6. 12 – Input per LS-Dyna dell’equazione di stato Tabulated Compaction.................................................................................................................................. 159

Tabella 7. 1 - Risultati numerici della terza Prova Sperimentale realizzata. 170

Tabella 7. 2 - Parametri che caratterizzano il modello in LS-Dyna Version 970. .. 175 Tabella 7. 3 – Ampiezza numerica della Fessura misurata nel centro della piastra all'intradosso. ........................................................................................................... 175 Tabella 7. 4 - Spostamenti residui, centrali e di bordo, misurati all’estradosso e all’intradosso della piastra. ...................................................................................... 176 Tabella 7. 5 - Elementi nodali centrali e di bordo sulla piastra discretizzata in SPH.................................................................................................................................. 177 Tabella 7. 6 - Parametri che caratterizzano il modello in LS-Dyna Version 970. .. 180 Tabella 7. 7 – Ampiezza numerica della Fessura misurata nel centro della piastra.................................................................................................................................. 181 Tabella 7. 8 - Spostamenti residui, centrali e di bordo, misurati all’estradosso e all’intradosso della piastra. ...................................................................................... 181

Tabella 8. 1 – Variabilità dei risultati per il metodo FEM, al variare delle dimensioni della Mesh e dei parametri b1, b2 e b3. 187

Tabella 8. 2 - Variabilità dei risultati per il metodo FEM, al variare delle dimensioni della Mesh e dei parametri b1, b2 e b3. ..................................................................... 187


195

Bibliografia

1. Camillo Belgrano. Gli esplosivi. Ed. Hoepli, Milano, 1973.

2. C. Giorgio. Tecnica degli esplosivi-Impiego degli esplosivi (2

vol.). Ed. Bianco, Udine, 1964.

3. A.Bacci, L’Esplosivistica in Campo Civile, pubblicato sugli Atti

del “I° seminario Nazionale sugli Ordigni Esplosivi”, svolto

presso il Comando Generale dell’Arma dei Carabinieri - III°

Reparto Armamenti ed Equipaggiamenti Speciali.

4. H. von Dach, Der Totale Widerstand, Ed.SUOV, Biel, 1972

5. Pedro F.Silva, Binggeg. Lu & Antonio Nanni. Prediction of

blast loads on the expected damage level by using displacement

based method. Department of Civil, Architectural and

Environmental Engineering, University of Missouri-Rolla, USA.

6. Francesco Cesari. Introduzione al Metodo degli Elementi Finiti.

Ed. Pitagora, Bologna, ISBN 88-371-0013-2, 1989.

7. Shaofan Li, Wing Kam Liu, Meshfree Particle Methods.

Springer, 2005.

8. Shaofan Li, Wing Kam Liu, Meshfree and Particle Methods and

their applications. American Society of Mechanical Engineers,

Appl Mech Rev vol 55, no 1, January, 2002.

9. M.B. Liu, G.R. Liu, K.T. Lam, A one-dimensional Meshfree

particle formulation for simulating shock waves. Shock Waves,

vol 13, pag 200-211, 2003.

10. Lu Young, Xu Kai, Numerical Simulation of Concrete Break-up

under Blast Loading. Civil Engineering Research, January, 2004.

11. LS-DYNA, Theory Manual 2006 Version 970, Livemore Software

Technology Corporation, LSTC Report.


196

12. LSTC (1999). LS-DYNA keyword user’s manual, non-linear

dynamic analysis of structures, version 970. Livermore Software

Technology Corporation, LSTC report.

13. Malvar L J, Crawford J E, Wesevich J W, Simons D: A

Plasticity Concrete Material Model for Hydrocodes. International

Journal of Impact Engineering, Vol. 19, No. 9-10, pp 847-873.

Elsevier Science Ltd. 1997.

14. Leonard E Schwer, L. Javier Malvar, Simplified Concrete

Modeling with Mat_Concrete_Damage_Rel3. JRI LS-Dyna User

Week, August, 2005.

15. Bažant Z P, Fracture Mechanics of Concrete Structures, State-of-

the-Art and Proc. 1st Int. Conf. Breckenridge, CO, 1992, Elsevier

Applied Science, Amsterdam.

16. Honggun Park, Jae-Yo Kim. Plasticity model using multiple

failure criteria for concrete in compression. International Journal

of Solids and Structures , Volume 42, Issue 8 , April 2005, Pages

2303-2322.

17. Ezio Cadoni, Kamil Labibes, Mario Berra, Marco Ginagrasso,

Carlo Albertini. Influence of Aggregate size on Strain-Rate

Tensile, Behavior of concrete. ACI Material Journal, no 98-M24,

May-June, 2001.

18. Malvar L J, Crawford E J, Dinamic Increase Factor for

Concrete. Twenty-Eighth DDESB Seminar. Orlando, FL, August

1998.

19. Mattias Unosson, Numerical simulations of penetration and

perforation of high performance concrete with 75mm steel

projectile. Defence Research Establishment Weapons and


197

Protection Division SE-147 25 TUMBA, SWEDEN. November

2000. ISSN 1104-9154.

20. Comité Euro-International du Béton, CEB-FIP Model Code

1990, Redwood Books, Trowbridge, Wiltshire, UK, 1993.

21. Bažant Z P, Prat C, Microplane Model for Brittle-Plastic

Material: Theory and Verification. Journal of Engineering

Mechanics, Vol. 114, No. 10, October, 1988. ASCE, ISSN 0733-

9399/88/0010. Paper No. 22823.

22. Comité Européen du Béton, CEB-FIP Model Code 90, Thomas

Telford Service, London, 1993.

RRIINNGGRRAAZZIIAAMMEENNTTII

198

Ringraziamenti

Il conseguimento di un obiettivo importante è difficile da raggiungere,

quale una Laurea in Ingegneria, è una gioia che va condivisa con le

persone più care. Appaiono scontati i ringraziamenti ai propri genitori

per gli sforzi profusi nel corso di questi anni, per consentirmi di tagliare

quest’importante traguardo, più appropriati sono invece i

ringraziamenti nei loro confronti per la pazienza e la capacità di

sopportare i tanti momenti di nevrosi, di stress ed esaurimento, che

hanno pervaso il mio corso di studi. A tal proposito mi sembra più che

opportuno ringraziare le mie care ed adorate sorelline, Rossana e

Ylenia, costrette a subire la medesima sorte. Come non poter ringraziare

la mia fidanzata Fabiana, per tutte le volte che mi è stata vicina

ascoltando le mie telefonate e leggendo i miei messaggi pregni di ansie e

paure prima di ogni esame. Un amico sempre fedele, capace di farmi

ritrovare il sorriso nei momenti difficili, e che per tanto mi sento di

ringraziare è il mio Beagle, nonostante le tavole dei progetti, gli

appunti, le matite e le penne da lui distrutte.

Lunga sarebbe la lista di amici e parenti da ringraziare, compagni di

studi con i quali in questi anni ho condiviso le gioie per gli esami

superati,e qualche immancabile amarezza.

Ma in particolar modo ringrazio gli amici di sempre, i cari vecchi

compagni del liceo, Mauro, Ferruccio, Gennaro, Luca, Marcello, Rocco,

Mino, Fonzie e tutti gli altri che non me ne vogliano se non nominati,

per tutte le serate passate insieme, le partite di pallone, le feste, le notti

brave, le zingarate, e quanto altro fatto insieme che è servito a

distogliermi dalle fatiche universitarie. Amici con i quali spero di

RRIINNGGRRAAZZIIAAMMEENNTTII

199

tornare a vivere momenti di spensieratezza, che tanto mi sono mancati

in questi ultimi mesi di sacrifici.

Un ringraziamento di vero cuore lo devo alle mie care nonne, a loro che

hanno sempre pregato nella speranza di potermi presto vedere laureato,

sono consapevole di dare loro e alla mia famiglia un enorme gioia. Il

mio più grande rammarico è che i miei due meravigliosi nonni, entrambi

di nome Giuseppe come me, non siano partecipi di questa stessa gioia. A

loro va la mia gratitudine, per tutte le volte che mi hanno protetto ed

aiutato, per tutte le preghiere che hanno sempre ascoltato.

Infine, ma non per questo meno importante, non posso esimermi dal

ringraziare il mio amico non che correlatore Ing. Domenico Asprone,

per la pazienza, il tempo e la disponibilità che mi ha concesso. È stata

per me una vera fortuna incontrarlo e conoscerlo nel momento più

importante del mio percorso universitario. Insieme a lui ringrazio il

Prof. Ing. Gaetano Manfredi e il Prof. Ing. Andrea Prota, per

l’opportunità concessami di prendere parte a questo interessantissimo

studio di ricerca. Ringrazio l’intero Dipartimento di Ingegneria

Strutturale, dottorandi, ricercatori e tecnici di laboratorio, per la

simpatia dimostratami.

2007.DellaPorta

Documents

Transcript of 2007.DellaPorta