Analisi di Buckling

Come detto in precedenza esistono due procedure attivabili per il calcolo dei carichi di collasso

La prima si basa sulla ricerca del punto di biforcazione mediante l’analisi modale

La seconda utilizza il calcolo non lineare e determina il carico di collasso direttamente

valutando il valore massimo raggiungibile fino ad un cambio drastico di configurazione

deformata

Metodo non lineare: si opera in controllo di

carico (carico crescente) e il calcolo non

converge più a Fc

Oppure in controllo di deformazione e si

esamina l’intera curva di risposta

Fc

Calcolo autovalori determinando il punto di

biforcazione (come si vede il calcolo non è in

sicurezza in quanto si linearizza il

comportamento fino al punto di collasso)

Impostazione agli autovalori

0 *

geomelast KK

Il problema si risolve mediante un’analisi agli autovalori delle matrici

È quindi necessaria un’analisi preventiva statica non lineare che consente di calcolare la

K*geom (opzione di stress stiffening)

Viene di seguito presentato come si forma tale matrice in una trave in flessione

Per piccoli spostamenti, si può scrivere l’equazione di equilibrio dei momenti attorno al punto di

sinistra

P P

u

v

2

dxdM dT dvM M dx q T dx dx P dx 0

dx 2 dx dx

Eliminando i termini in dx2 , e derivando rispetto ad x: 0dx

vd P

dx

dT

dx

Md2

2

2

2

Non considerando le deformazioni a

taglio si ha dalla teoria della trave 2

2

dx

vd I EM q

dx

dT

qdx

vd P

dx

vd I E

2

2

4

4

Si può ora applicare la formulazione debole per la determinazione delle matrici di

rigidezza a partire dalla formulazione della deformazione troncata al II ordine

2 2 21

2xx

u u v w

x x x x

2

2

2

xxdx

dv

2

1

dx

vd y

dx

du

(termine flessionale) (grandi deformaz.)

eV

2

xx dV E2

1UEnergia di deformazione: Sostituendo e svolgendo prodotti

22 4 2 22 2 22

2 2 2

1 1 2

2 4L A

du d v dv du d v d v dv du dvU A y y y E dA dx

dx dx dx dx dx dxdx dx dx

Che è la classica equazione differenziale dalla quale si può

dedurre l’instabilità delle colonne, valida per ogni C.C.

1 2 3 4sin cosv x C kx C kx C x C + Condizioni al Contorno

Il precedente integrale si può ridurre ad un integrale

di linea ricordando le proprietà delle sezioni:

22 4 22

2

1

2 4L

du d v A dv du dvU E A I A dx

dx dx dx dx dx

E ricordando chedx

duAEP

20 ; 0 ;A A A

dA y dA y dA I

Trascurando il termine di potenza più alta

41

02 4

L

A dvE dx

dx

22 22

2

1

2x

L

du d v dvU EA EI P dx

dx dx dx

Si perviene alla seguente

espressione della energia elastica

L’energia elastica è suddivisibile in due componenti, il primo associato a deformazioni

assiali e il secondo a deformazioni flessionali

ax flU U U

21

2ax

L

duU EA dx

dx

2 22

2

1

2fl x

L

d v dvU EI P dx

dx dx

L

0

22

2

22

dx dx

dv P

dx

vd I E

dx

duA E

2

1U

WU 1 1

2 2

all all fl fl all

TT T

gf K f f K f f P

Nel problema discretizzato si può esprimere lo spostamento mediante le funzioni di forma, per

cui l’energia totale diviene:

1

2

ax ax ax

Tf K f

1

2fl fl fl

T

f K f 1

2

fl fl

T

gf K f

1 , flx

dv

dx

N f 1 1, ,0

A L T

x xdx

*

gK N N P *

g gK K

2 3 2 3 2 3 2 3

1 11 12 13 14 2 3 2 2 3 2, ,

,

N N N N 1 3 2 , 2 , 3 2 ,x x

x

x x x x x x x xx

L LL L L L L L

,xN

1 2 3 2 2 3 2,

6 4 6 212 , 6 , 12 , 6

x

x x x x

L LL L L L L L

N

Derivando le precedenti si ha:

2 3

2

2 3 2 2 3 20

2 3

2

612

46

6 4 6 2 12 6 12 6 A

612

26

L

x

L L

x

x x x xL Ldx

x L LL L L L L L

L L

x

L L

*

gK

Per ricavare non resta che risolvere l’integrale

6 1 6 1

5 10 5 10

1 2 1

10 15 10 30

6 1 6 1

5 10 5 10

1 1 2

10 30 10 15

L L

LL

L L

LL

*

gK

1 1 2 2 v v

1

1

2

2

v

v

*

gK

Si noti che sono presenti solo gdl flessionali

Esistono due modalità di evoluzione non lineare, la prima facendo crescere il carico e la

seconda lo spostamento

Soluzione non lineare

In realtà ne esiste anche una terza, combinazione delle due precedenti, che prende il

nome di Arch-length od analisi di Riks

Il calcolo non lineare si

muove su questo arco

Il calcolo a collasso non lineare è a tutti gli effetti un’analisi non lineare tout court e

quindi si rimanda l’approfondimento a quello del calcolo non lineare in genere

Risoluzione di problemi strutturali non lineari

Dal punto di vista strutturale, i problemi sono non lineari se, nota la risposta del sistema ad un

determinato livello di carico, incrementando tale livello i risultati non sono proporzionali

all’incremento stesso

Se chiamiamo step l’intera soluzione, dal momento in cui non è applicato il carico al suo valore

finale, è quasi sempre necessario suddividere lo step in un certo numero di substep, ove si

imponga la soddisfazione dell’equilibrio (forze residue piccole)

Anche all’interno dei substep il calcolo è non lineare, quindi sarà necessario operare un certo

numero di iterazioni all’interno (ciascuna costituita da un calcolo lineare) finché non si avrà al

corrispondente substep una soluzione con forze residue minori del voluto

tempo

risposta

Substeps

STEP

iterazioni

In ogni caso, il calcolo non lineare si sviluppa come un’analisi incrementale, del tipo:

1) In controllo di carico – i carichi sono applicati con una legge preimpostata

2) In controllo di spostamento – spostamenti imposti indipendentemente dalla risposta sistema

3) Controllo indiretto sulla risposta – controllo basato su combinazioni della risposta, come ad

esempio spostamenti relativi tra gdl

3) Controllo di tipo Arc length – si realizza mediante una combinazione di spostamento e carico

applicato la cui combinazione deve muoversi entro un raggio funzionale delle due grandezze

Il modello costitutivo del materiale fa ricorso alla matrice di rigidezza tangente del materiale

* ε=ε

σD ε

ε

ε=ε

σσ ε ε σ ε ε

ε

Attraverso la quale si calcola la matrice di rigidezza tangente di ogni elemento della struttura

*int T T T

T d d d

f σ ε

K B σ B B D Bu u ε u

8) Fine del substep con soluzione

convergente e pronti ad un nuovo

incremento del carico

5) Si trova 2 1 1 n n nu u u

6) Si determinano le forze interne, si

calcola il nuovo residuo associato allo

spostamento estrapolato

int

2

2

3 2 2

n

T

R

n 1 n n

f

K

u u u

7) . . . Si prosegue finché il

residuo risulta essere più grande

di un valore considerato

accettabile

Il metodo di Newton-Raphson o

della rigidezza tangenziale ricerca la

soluzione dall’inizio di un substep

alla sua fine calcolando ad ogni

iterazione la matrice tangente ed

estrapolando su di essa

Per semplificare la trattazione, si espone la procedura per un sistema a risposta non lineare

che però possiede un solo grado di libertà

3) Si determina la nuova matrice

di rigidezza

1

TK

2) dal substep precedente è nota 1

last

1

n nu u

1

1 1 1

nR

n TAN

u K4) Si calcola per estrapolazione

lineare l’incremento spostamento

1

last

1

n nu u

1) l’incremento di carico imposto

definisce il valore del residuo iniziale1

n nR P

Il metodo di Newton-Raphson modificato differisce dal precedente in quanto non si ricalcola la

matrice tangente ad ogni iterazione ma si utilizza quella relativa alla fine del substep precedente

Il metodo di Newton Raphson può

entrare in crisi se la struttura è hardening

- la convergenza diviene lentissima

P

x0 x1x2x3x4

Se il sistema ha comportamento non monotono si

può mancare totalmente la convergenza anche col

metodo non modificato

Miglior tendenza con N.R.

modificato

P

Caso di sistema hardening

Non convergenza con Newton-

Raphson modificato

Convergenza lenta con Newton-

Raphson

Nel metodo della secante si effettua una prima iterazione secondo il metodo classico di Newton-

Raphson e poi si determina la matrice secante mediante

La determinazione della matrice secante è immediata in un problema scalare, come

quello rappresentato in figura, ma molto meno banale per rigidezze matriciali

1

1 1 2 1 2

SEC n nd R R

nu K

Esiste anche un metodo della secante modificato (Iterazione diretta di Picard) secondo il quale la

secante è computata sempre considerando la soluzione iniziale

Il miglior funzionamento delle varianti presentate dipende … dal problema stesso, e quindi è

di difficile determinazione aprioristica

L’adaptive descent può fare uso di due matrici di rigidezza incrementali, una è la matrice

tangente Kt e l’altra la matrice secante Ks

La matrice in uso è definita dalla espressione ts KKK 1

L’iterazione ha inizio con = 0, se la convergenza è difficile da raggiungere si fa crescere

in modo da avvicinarsi alla matrice secante, più rigida se il sistema è hardening

Se nell’iterazione il residuo cresce (tendenziale divergenza):

Si porta ξ = 1 e si continua ad iterare così

Se nell’iterazione il residuo decresce (tendenziale convergenza):

Si riduce progressivamente ξ da 1 fino a riportarlo a 0 (matrice tangente)

Se durante l’iterazione si verifica un pivot negativo (matrice mal condizionata):

Si riduce riutilizza la matrice secante (ξ = 1) ed eventualmente si riduce la

dimensione del substep di avanzamento

Il line search è un metodo che prova a determinare la soluzione della iterazione modificando

i salti calcolati da Newton Raphson con l’adozione di un coefficiente sk appropriato

1 i i i

n n k ns u u u con0 1.0s

Noto il valore di sk si associa ad esso uno scalare che indica l’energia coinvolta nell’iterazione

1,T

i i k

k n ng u R

Ove naturalmente il residuo è calcolato dall’equilibrio delle forze esterne ed interne

1, int 1,i k ext i k

n n

R f f u

A questo punto si determina la nuova approssimazione di

sk+1 estrapolando verso g=0 la retta indicata nella figura 0g

ks

kg

1ks

s

g

01

0

k k

k

gs s

g g

Il processo termina dopo un numero prefissato di

iterazioni (e.g. 5) oppure quando si verifica una delle due

0 0.5k

g

g 1k k

k

g gerr

g

L’arc length è un metodo molto adatto quando i sistemi sono softening-hardening, per esempio

nelle condizioni di post-buckling

A seconda del metodo adottato, si realizzano perdite di convergenza diverse

La curva degli equilibri successivi in pratica si determina da una combinazione sia degli

spostamenti delle variabili incognite, sia del fattore scalare di carico che agisce su tutti i carichi

applicati

Il metodo non è invece adatto nei sistemi ove discontinuità di carico si realizzano

frequentemente durante gli step di carico (analisi con contatti tra superfici, …)

Dato però che spostamenti e forze applicate sono dimensionalmente differenti, occorre definire

una costante di comparazione c:

2Tl c u u

Per gli elementi in cui sono presenti sia gdl

traslazionali che rotazionali si introducono

anche altri fattori di comparazione

In un substep, si parte da 1 1 , n n u

0u

e si vuole determinare i nuovi , n nu

In pratica si ha una nuovo parametro con la condizione vincolare sulla lunghezza din l

Riscriviamo la linearizzazione dell’equilibrio della generica iterazione i :

1 1 1

0 int, 0i i i i i

n n n n nf f f f K u

1 1 1

0 int,

i i i i i

n n n n nK u f f f f

INCOGNITE

Ora, pensando che una iterazione porta da una soluzione ad una possiamo

linearizzando imporre l’equilibrio nella condizione iniziale (0) e finale (f)

1

0=i

n u u fu

1 1 1

0 0 int,

1

i i i

n n n

i

n f

K u f f f

K u f

Esprimendo la

variazione di

spostamento come

Sia che sono note e quindi, sostituendo nella l’unica incognita presente èfu

0 i i

n n f u u u

i

n

2

1 1 2 1 2T

i i i i i i

n n n n n nc l u u u uIl sistema, di II grado, fornisce due

soluzioni: una avanza nel percorso,

l’altra retrocede

Si separa

l’equazione

precedente in due

Per quanto riguarda i criteri di convergenza essi possono essere basati sia sulle forze

(e momenti) residui sia sullo spostamento della soluzione dall’iterazione precedente

I primi sono senz’altro da preferire in quanto hanno un preciso significato fisico, i secondi

invece possono mascherare enormi errori nell’equilibrio per strutture particolarmente rigide

Infine, in genere i codici controllano se da una iterazione alla successiva si ha una nuova

condizione di contatto ed in tale caso effettuano sempre una iterazione aggiuntiva,

fermandosi quando non si altera più il quadro complessivo dei contatti

La ricerca della soluzione può anche diventare estremamente difficile quando sono presenti

condizioni non continue, come contatto o bruschi cambi di pendenza curva plastica del

materiale, passaggi di stato, ...

In particolare i problemi di contatto, peggio se associati a plasticità che invece produce

softening sono spesso difficoltosi perché instaurano delle condizioni di hardening

all’insorgere della penetrazione sul target

Comportamento plastico dei materiali (metallici)

La plasticità è caratterizzata da una correlazione non

biunivoca tra il valore della deformazione accumulata e

quello della tensione, in pratica si può vedere se in un ciclo

di carico e scarico si ripercorre o no il medesimo percorso

Nei materiali elasto-plastici perfetti esiste un valore di

tensione di yield oltre il quale le deformazioni sono

indeterminate

carico

scarico

y

I materiali che incrudiscono sono invece caratterizzati

da un valore di yielding che dipende da qualche

parametro di controllo

Il più accreditato è la

deformazione plastica

L’equazione della superficie di Yield è 0,F σ

Nei materiali che presentano incrudimento la superficie

tende ad espandersi mano a mano che si accumula

plasticità

1

2

Analisi del flusso plastico

Nel caso di presenza di flusso plastico il materiale è caratterizzato da un legame tensione –

deformazione che si modifica durante l’evento plastico.

In una semplificazione comunemente adottata, il materiale all’interno di una superficie di

snervamento si comporta elasticamente, al di fuori plasticizza

Se incrudente, si ha una estensione della superficie di snervamento stessa determinata dall’entità

della deformazione plastica stessa.

Chiamando con F la superficie di snervamento, essa sarà funzione dello stato di tensione e di un

parametro che tiene in conto dell’incrudimento subito

, 0F σ

Il tensore delle tensioni è naturalmente il seguente:

xx xy xz

xy yy yz

xz yz xx

σ

Considerando la tensione media, si può scrivere un tensore

deviatorico che si costruisce eliminando la componente idrostatica: 1 3m x y z

xx m xy xz

xy yy m yz

xz yz xx m

s

(1)

Si rammentano anche i tre invarianti delle tensioni che sono così definiti, una volta che è stato

individuato il sistema di riferimento principale

1 1 2 3I σ 2 1 2 2 3 3 1I σ 3 1 2 3I σ

Nel caso si faccia uso del tensore deviatorico, i tre invarianti assumono i valori

1 0J σ 2 2 2

2 1 2 2 3 3 1

1

6J

σ 3 1 2 3J s s sσ

Il criterio di snervamento di Von Mises è imediatamente legato a J2

2 2 2

1 2 2 3 3 1

1

2eq 23 eq J σ

Qualora ci si trovasse in un sistema di riferimento non principale, il J2 diverrebbe

2 2 2 2 2 2

2

16

6x y y z z x xy yz zxJ

σ

La deformazione plastica equivalente viene anche essa in

genere conteggiata mediante un criterio alla Von Mises 2

33

2

pl pl pl pl

eq ij ijJ ε

Un metodo per determinare le direzioni principali deviatoriche, alternativo alla soluzione del

problema agli autovalori, è dato dalle seguenti espressioni (Kachanov):

2 2

2 1cos

33s J

3 2

2cos

3s J

3

32

2

3 3 cos3

2

J

J

;

;

Avendo posto

A questo punto le tensioni principali altro non sono che

1 1 ms

2 2 ms

3 3 ms

;

;

Tornando alla superficie di yielding F, se si adotta la tensione equivalente di Von Mises, l’eq (1) è

, 0eq yF σ 2 2 2

1 2 2 3 3 1

10

2y

In questa ultima equazione l’effetto dell’incrudimento viene conteggiato mediante la variazione

della tensione di snervamento, che sarà funzione di un parametro da definire ma che tiene in

considerazione il livello di plasticità raggiunto

Incrudimento isotropico

1 2

2 1cos

33s J

Superficie di snervamento

Materiale Elasto-Plastico perfetto

La superficie di snervamento definisce l’innesco della condizione di plasticità

L’estensione di tale superficie si attualizza con la deformazione plastica accumulata

, 0plF σ ε

Incrudimento isotropico

Incrudimento cinematico

In questo caso la condizione di plasticizzazione è totalmente indipendente

dalla deformazione plastica

0F σ

La condizione di plasticizzazione rimane la stessa nel corso della deformazione plastica, sia nella

estensione che nella posizione originaria

La superficie di snervamento può crescere ma non cambia posizione

La superficie di snervamento può solo traslare ma non

cambiare di forma

La superficie di snervamento segue l’evoluzione dello stato tensionale in modo che la disequazione

sia sempre rispettata

Condizione di consistenza

, 0plF σ ε

Un cambiamento dello stato tensionale si sviluppa elasticamente (scarico elastico) se , 0pldF σ ε

Si ha invece ulteriore plasticizzazione se si realizza una fase di carico

ulteriore con incremento della deformazione plastica , 0pldF σ ε

Materiali incrudenti

La superficie cresce in estensione dal

primo snervamento fino alla superficie

di rottura (il cammino è irreversibile)

Principio di normalità

Esso asserisce che in condizioni di incipiente uscita dalla curva limite di snervamento, la

deformazione plastica fluisce secondo una direzione che risulta normale alla curva limite stessa

In termini di equazioni, si può scrivere che, per ciascuna delle 6 componenti indipendenti

dell’incremento di deformazione plastica, si ha la condizione di parallelismo tra l’incremento plastico

di deformazione e la normale alla superficie di snervamento

pl

ij pl

ij

Fd

(2)

La costante è per ora niente di più che un fattore di proporzionalità che andrà determinato.

Nella plasticità associativa, di cui qui si discute, si ha che F risulta essere proprio la superficie di

snervamento prima definita

Se si considera plasticità non associativa (valida per esempio per materiali porosi che risentono della

tensione media) , allora occorre sostituire a F una funzione che definisca più appropriatamente il

potenziale plastico

tot el pl

ij ij ijd d d Per ottenere un legame costitutivo coerente, nel caso elasto-plastico, si

suddivide la deformazione in una somma di contributo elastico e plastico

Ricordando il principio di normalità ed il legame tensioni-

deformazioni in campo elastico e plastico1 F

d d

ε D σ

σ(3)

Dove D - legame elastico - vale, nel caso tridimensionale

1 - - 0 0 0

- 1 - 0 0 0

- - 1 0 0 01

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 0 0 1

xx xx

yy yy

zz zz

xy xy

yz yz

xz xz

d d

d d

d dd d

d dE

d d

d d

ε D σ

Come detto, quando si realizzano le condizioni di flusso plastico, la superficie limite si incrementa

(incrudimento) e fa si che essa continui ad inglobare lo stato tensionale raggiunto.

Questo vuol dire che il suo differenziale totale (incremento assoluto), composto da 7 termini, rimane

sempre nullo

0xx yy z xy yz zx

xx yy zz xy yz zx

F F F F F F FdF d d d d d d d

(4)

Definendo ora con 1 F

A d

(5)

Utilizzando il vettore a 6 componenti di , l’equazione scalare (4) assume una forma semplificata

0

TF

d A

σσ

(6)

Se ora si raggruppano le eqq. (3) e (6) si perviene ad un sistema composto di 7 eq. in 7 incognite:

1

0T

Fd

d

Fd A

D σσε

σσ

(7)

Vediamo ora di discutere il significato del termine A. Chiaramente, se non si ha incrudimento,

esso si annulla in quanto F rimane lo stesso al crescere di

Il punto di partenza è la determinazione del parametro che misura il lavoro

incrementale di hardening, per cui la sua crescita vale

T pld d σ ε

... pl pl pl

xx xx yy yy zx zxd d d d

Ricordando la legge di flusso pl

ij

ij

Fd

(2) e sostituendola in quest’ultima T F

d

σσ

Che ci consente di eliminare il fattore d presente nella equazione (5)

TF FA

σ

σ(8)

Von Mises presenta il notevole vantaggio di poter esprimere il gradiente di F in modo agevole

2 2 2 2 2 21

3 3 3 02

xx yy yy zz zz xx xy yz zx yF σ

33 3; ; ;

2 2 2

yx z

xx eq yy eq zz eq

ss sF F F

3 3 3; ; ;

xy yz zx

xy eq yz eq zx eq

F F F

Con semplici passaggi si scrivono le 6 derivate

Se si dispone di una semplice prova di trazione (monodimensionale), si ha a

disposizione una curva del tipo m - m , spostandosi dalla condizione di

incipiente snervamento, il lavoro incrementale di deformazione plastica è: y md d

Si esprime - in questo caso monodimensionale - il valore della prima derivata scalare presente

nella eq. (8) avendo bene in mente che il differenziale totale (4) risulta nullo

0s

s

F Fd d

1

y y

y m y

d dF H

d d

=1

y

m

dH

d

Avendo definito H non è altro che la pendenza istantanea

della curva di incrudimento monotona.

Se ora si sostituisce nella (8) ricordando anche che nel caso monodimensionale la

seconda derivata ivi presente vale semplicemente 1, si perviene ad una definizione di

A H

TF FA

σ

σ

Si può ora tornare al sistema (7)

cercando di darne una soluzione il più

possibile generale. In particolare

occorre approntare una opportuna strategia risolutiva che

rimanga valida anche quando A = 0, ossia il materiale

presenti un comportamento elasto-plastico perfetto

Lo scopo sarà quello di fare scomparire nel sistema (7) il termine λ.

Si comincia col premoltiplicare il primo set di equazioni (7) per il termine T

F

Dσ

1

T T TF F F F

d d

D ε D D σ D

σ σ σ σ

Semplificando, e mettendo in evidenza il primo termine a destra dell'uguale, si ottiene il seguente

T T TF F F F

d d

σ D ε Dσ σ σ σ

termine che può essere poi sostituito nel II set di equazioni (7) - che è poi

una semplice eq. Scalare - si elimina così il dσ0

TF

d A

σσ

0

T TF F F

d A

D ε Dσ σ σ

1 Fd d

ε D σσ

Quest’ultima permette di eliminare il fattore moltiplicativo

T

T

Fd

F FA

D εσ

Dσ σ

(9)

(6)

Che viene infine eliminato nella prima delle (7) 1 Fd d

ε D σσ

1

T

T

Fd

Fd d

F FA

D εσ

ε D σσ

Dσ σ

(10)

Ora si scrive la precedente equazione in termini espliciti, indicando cioè la variazione dello

stato tensionale per effetto di un incremento di deformazione totale (si premoltiplica per D)

T

T

F F

d d dF F

A

D Dσ σ

σ D ε ε

Dσ σ

(11)

In pratica si è trovato un legame incrementale tra il tensore di deformazione e quello della tensione

che si scrive noto che siano la matrice D del materiale elastico, l’ipotesi di rottura attraverso la

definizione di F e quindi delle sue derivate, la pendenza incrementale di incrudimento plastico (A).

*d dσ D ε *

T

T

F F

F FA

D Dσ σ

D D

Dσ σ

Questo modo di procedere fin qui descritto ha lo svantaggio di considerare una D* indipendente

dallo stato di tensione che di fatto fuoriesce istantaneamente dalla curva di snervamento

Se gli incrementi di deformazione sono sufficientemente piccoli l’errore che si commette non

è grande, tuttavia se così non è l’errore accumulato può diventare intollerabile.

Per il momento il sistema non viene descritto, si ritiene utile solo accennare al fatto che la

risoluzione comporta la ricerca di zeri della funzione attraverso l’algoritmo di Newton-Raphson.

In questo caso si può far ricorso all’algoritmo del return mapping, proposto nel 1964 da

Maenchen e Sacks. La tecnica prevede uno scaling della tensione in modo da permanere

sempre sulla superficie di incrudimento

Qualora la plasticità non fosse asssociativa, ugualmente la ricerca del punto finale sulla superficie

di snervamento va risolta in forma numerica

Esempio su deformazione 3D

Si ipotizza la presenza di un determinato stato di tensione e di

deformazione, giacenti sulla superficie limite di snervamento:

275.04 59.01 -202.86

59.01 260.36 -27.52

-202.86 -27.52 -65.40

σ

La quale presenta, come si può risolvere, tre tensioni principali pari a:1 400 MPa

2 230 MPa

3 160 MPa

Secondo i criterio di Von Mises, essendo la tensione equivalente sulla curva di snervamento, si ha

una tensione di snervamento 2 2 2400 230 160 400 230 230 160 160 400 497.29y

Si calcola il tensore di

deformazione elastico, ipotizzando

di essere in regime elastico e

quindi di non aver ancora compiuto

lavoro plastico (=0):

1 - - 0 0 0

- 1 - 0 0 0

- - 1 0 0 01

0 0 0 1 0 0

0 0 0 0 1 0

0 0 0 0 0 1

xx xx

yy yy

zz zz

xy xy

yz yz

xz xz

E

Ponendo per il materiale le seguenti caratteristiche: E= 2.06e11 Pa, =0.3, sn=500 Mpa

H= 1/1000*E

1 - - 0 0 0 275.04 1.0512

- 1 - 0 0 0 260.36 0.9586

- - 1 0 0 0 -65.40 -1.09721

0 0 0 1 0 0 59.01 0.3724

0 0 0 0 1 0 -27.52 -0.1737

0 0 0 0 0 1 -202.86 -1

xx

yy

zz

xy

yz

xz

E

310

.2802

Da un calcolo complessivo, a partire da questo stato, si ha un incremento di deformazione totale

(di primo tentativo) ad esempio di semplice dilatazione unidirezionale in x , pari a 0.001:

3

1.0512 1.0000 2.0512

0.9586 0.0000 0.9586

-1.0972 0.0000 -1.097210

0.3724 0.0000 0.3724

-0.1737 0.0000 -0.1737

-1.2802 0.0000 -1.2802

new old d

ε ε ε310

La prima verifica da fare è che l’incremento di deformazione totale faccia sforare la sollecitazione

equivalente oltre il limite di snervamento attuale:

3

1- 0 0 0 2.0512 552.39

1- 0 0 0 0.9586

1- 0 0 0 -1.0972 10

0 0 0 1 2 4 0 0 0.37241 1 2

0 0 0 0 1 2 4 0 -0.1737

0 0 0 0 0 1 2 4 -1.2802

xx

yy

zz

xy

yz

xz

E

379.21

53,44

59.01

27.52

202.86

573.31 MPa 497.29 MPaeq La sollecitazione esce dalla zona elastica

Con il che si evince che si è entrati in campo plastico e bisogna utilizzare la equazione (11) per il

calcolo del tensore delle tensioni.

Per poter applicare la (11) occorre calcolare le derivate del gradiente di F rispetto alle tensioni –

nell’ultimo punto noto in elasticità:

+ 3=156.67 m x y z MPa

=118.37 x x ms MPa

=103.70 y y ms MPa

= -222.07 z z ms MPa

= 59.01 xy xys MPa

= 27.52 yz yzs MPa

= -202.86 xz xzs MPa

3 222.073 118.37 3 103.700.3571; 0.3128; 0.6698;

2 497.29 2 497.29 2 497.29xx yy zz

F F F

3 27.52 3 202.863 59.010.3560; 0.1660; 1.2238;

497.29 497.29 497.29xy yz zx

F F F

0.3571 0.3128 -0.6698 0.3560 0.1660 1.2238

TF

σ

Calcolando ora la (12) si ottiene

* 5

2.6863 1.1124 1.3513 -0.0865 0.0404 0.2975

1.1358 2.7065 1.3311 -0.0758 0.0354 0.2606

1.4816 1.3063 2.4676 0.1624 -0.0757 -0.55821 10

-0.1559 -0.0627 0.3485 1.4983 0.0402 0

T

T

F F

F FA

D Dσ σ

D D

Dσ σ

.2966

-0.2739 -0.1101 0.6124 -0.3256 1.5659 -0.1383

0.0727 0.0292 -0.1625 0.0864 0.1518 0.5649

Per cui, osservando che *d dσ D ε

* 268.63 111.24 135.13 8.66 4.04 29.75 T

Td d MPa σ D ε

Si determina anche il nuovo stato di tensione

543.67 371.61 69.73 50.36 23.11 173.11 T

d MPa σ σ σ

Con una nuova tensione di snervamento 521.38 y MPa

Evidentemente più bassa di quella estrapolata al primo tentativo elastico (573.31 MPa)

Il parametro λ può essere ora calcolato dalla (9): 4

51.53 10

T

T

Fd

F FA

D εσ

Dσ σ

E il lavoro incrementale di hardening (per unità di volume) assume il valore

'4 = 51.53 10 543.67 371.61 69.73 50.36 23.11 173.11

0.3571 0.3128 -0.6698 0.3560 0.1660 1.2238 0.0763

T Fd

Nmm

σσ

Il punto debole del calcolo precedente è la determinazione del gradiente di tensione, preso nello stato

precedente alla deformazione aggiuntiva

Qualora tale gradiente dovesse modificarsi in modo significativo per effetto della deformazione

imposta dallo step, occorrerà ridurre gli intervalli di crescita di deformazione, in modo da inseguire

l’evoluzione dell’incrudimento restando sempre sulla superficie (aggiornata) di snervamento stessa

Questo modo di procedere risponde al nome di Return Mapping ed esistono svariati modi di

procedere in funzione anche del metodo utilizzato per calcolare la tensione equivalente

Nel calcolo agli elementi finiti, il modo di precedere si inserisce negli step di calcolo non lineare,

come ad esempio nello svolgimento del metodo di Newton Raphson.

Tale metodo prevede lo svolgimento di iterazioni, e quindi l’assunzione di un campo di spostamenti

di tentativo, che viene via via affinato finché non si raggiunge un equilibrio soddisfacente tra le

forze esterne applicate e le forze interne che si sviluppano

Il problema del calcolo elasto-plastico del materiale nasce proprio per la valutazione delle forze

interne, la cui definizione richiede la soluzione del legame non lineare tra deformazioni e tensioni

Pertanto, ad ogni iterazione si considera un sottoincremento dello spostamento

In modo da ottenere

i i

u u

ii ε ε e quindi ricercare le conseguenti variazione di tensione

Questa operazione è ripetuta per ogni elemento, in corrispondenza ad ognuno dei suoi punti

di Gauss, utilizzati per l’integrazione numerica

In pratica, la tensione, dipendente dal valore assunto al passo precedente, si risolverà da una

sommatoria che contiene integrali del tipo

*

1

0

i

n n n n iid

ε

σ σ σ σ D ε

Il modo esplicito più semplice di affrontare questo problema è di approssimare l’integrale

considerando un valore costante nell’arco di integrazione. Questo metodo è molto rapido, ma

decisamente poco preciso a meno che non si suddivida in intervalli molto piccoli e di dimensione

difficilmente prevedibile

In alternativa, si può determinare ciascun integrale con un approccio alla Runge-Kutta del II

ordine, che conduce ad una doppia stima della matrice tangente *

iD

Si impone una deformazione pari alla metà di quella effettiva e si

calcola la tensione relativa*

1 2 0

1

2 σ D ε

Con questo incremento di tensione si ridetermina la matrice tangente

che viene utilizzata per il calcolo dell’intero stato tensionale*

1 0 1 2 σ D ε

Questo metodo risolutivo consente anche di avere una stima

dell’errore mediante l’espressione1 1 22 err σ σ

Il metodo può comportare il continuo scostamento del tensore della tensione dalla curva di Yielding

In alternativa è stato sviluppato il metodo del Return-Mapping (1964), del quale si può dare una

semplice interpretazione geometrica dalla figura sottostante

In pratica si impone di riportare lo stato di tensione sulla curva di snervamento, mediante tratti di

ritorno tutti ortogonali e di livelli di incrudimento decrescente. L’algoritmo impone che lo stato di

tensionale finale corrisponda ad una tensione equivalente pari a quella di snervamento attualizzata

I dettagli operativi dipendono dal modello costitutivo utilizzato per la valutazione della plasticità

Dal punto di vista del codice di calcolo agli elementi finiti, si possono comunque prevedere i

seguenti passaggi

Assunzione della deformazione totale: allo step n secondo un’iterazione del

procedimento non lineare di calcolonε1)

Calcolo della deformazione di trial eliminando la parte già

plastica al passo precedente

2)1

tr plast

n n n ε ε ε

3) Determinazione della tensione di trial che si avrebbe se il

materiale reagisse linearmente 1

tr plast

n n n σ D ε ε

4) Calcolo della tensione equivalente per valutare se si rimane in regime elastico o se

si rientra in un tratto di plasticità

Risoluzione, mediante la formulazione esplicita (oppure il return mapping) del

moltiplicatore plastico λn e quindi della parte plastica dell’incremento di deformazione

5)

plast

n n

F

ε

σ

Determinazione della componente elastica e plastica della deformazione totale6)

elast tr plast

n n n ε ε ε 1

plast plast plast

n n n ε ε ε

Calcolo dello stato tensionale presente (dovuto solo alle deformazioni elastiche)7)

elast

n nσ Dε

Aggiornamento dell’incremento del lavoro plastico di incrudimento e dell’incremento plastico

equivalente della deformazione (per utilizzare la curva monotona)

8)

T

n n n

Fd

σ

σ

1n n nd 2

3

pl plT pl

eq ε ε