Una formulazione FEM per la verifica a shakedown

Una formulazione FEM per la verifica a shakedown
Antonio, Salvatore Petrolo
Dipartimento di Strutture, Universit`
a della Calabria, 87030 Arcavacata di Rende (Cosenza), Italy
tel. +39 0984 494031 fax +39 0984 494045 e.mail petrolo @ unical.it
Indice
Introduzione
1 La
1.1
1.2
1.3
1.4
1.5
1.6
1.7
2
teoria classica dello shakedown
Il legame costitutivo elastoplastico perfetto . . .
Assunzioni sull’andamento dei carichi, l’inviluppo
Deﬁnizione di shakedown elastico . . . . . . . . .
Il teorema statico dell’adattamento . . . . . . . .
Il teorema cinematico dell’adattamento . . . . . .
Il fattore di sicurezza a shakedown . . . . . . . .
Alcuni commenti sul moltiplicatore di sicurezza .
. . . . . . . . . . . . . .
di sollecitazione elastica
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
. . . . . . . . . . . . . .
2 Riformulazione del problema in termini della funzione di shakedown
2.1 La funzione di snervamento a shakedown e il dominio elastico ridotto . . .
2.2 Il fattore di sicurezza a shakedown . . . . . . . . . . . . . . . . . . . . . .
2.3 Il problema di shakedown . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Lo schema return mapping per la funzione di shakedown . . . . . . . . . .
2.5 Propiet`
a notevoli dello schema . . . . . . . . . . . . . . . . . . . . . . . .
3 Un
3.1
3.2
3.3
3.4
3.5
metodo iterativo di soluzione
Generica formulazione in elementi ﬁniti . . . . .
Soluzione iterativa del problema dello shakedown
Convergenza dello schema iterativo . . . . . . . .
Convergenza alla soluzione di adattamento . . . .
Alcuni commenti . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4 L’algoritmo di rientro sul dominio ammissibile
4.1 Il dominio ridotto per una generica funzione di snervamento . . . . . . .
4.2 L’algoritmo di return–mapping in plasticit`
a multifalda . . . . . . . . . .
4.3 Il problema di rientro come problema di programmazione quadratica QP
4.4 Cenni ai metodi di soluzione per problemi QP e superiorit`
a dei metodi
proiezione nel caso in esame . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Il metodo di Golfbard e Idnani . . . . . . . . . . . . . . . . . . . . . . .
4.6 Cenno alla soluzione del problema senza linearizzare i vincoli . . . . . .
1
.
.
.
.
.
.
.
.
.
.
.
.
5
5
6
8
9
11
12
14
.
.
.
.
.
16
16
17
18
18
19
.
.
.
.
.
22
22
24
26
29
30
. .
. .
. .
di
. .
. .
. .
32
32
33
34
35
35
36
2
5 Esempio di implementazione: i telai in c.a.
5.1 L’elemento ﬁnito asta 3d . . . . . . . . . . . . . . . . . . . . .
5.1.1 La teoria di De Saint Ven`
ant . . . . . . . . . . . . . .
5.1.2 Generazione dell’elemento asta 3d . . . . . . . . . . .
5.1.3 Soluzione numerica delle equazione di De Saint Ven`ant
5.1.4 Analisi di sezioni in campo elastico . . . . . . . . . . .
5.2 Deﬁnizione del dominio elastico delle sezioni . . . . . . . . . .
5.3 Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 La trave appoggiata . . . . . . . . . . . . . . . . . . .
5.3.2 Telaio piano . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
38
38
39
41
46
49
55
58
58
59
6 Esempio di implementazione: stati piani di tensione e deformazione
6.1 L’elemento ﬁnito Simplex . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 La discretizzazione del dominio nel caso di stati piani di tensione e deformazione
6.3 Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3.1 Lastra quadrata con foro centrale circolare . . . . . . . . . . . . . .
6.3.2 Ulteriori test con concentrazione di tensioni . . . . . . . . . . . . . .
6.3.3 Alcuni test con collasso incrementale . . . . . . . . . . . . . . . . . .
62
62
65
66
67
70
71
Considerazioni conlcusive
74
Bibliograﬁa
75
Introduzione
La possibilit`
a, data dalla normativa, di eﬀettuare il dimensionamento delle strutture sulla base di veriﬁche di sicurezza condotte in campo nonlineare ha dato rilevanza pratica
all’analisi elastoplastica, in passato spesso conﬁnata in ambito unicamente accademico.
Nel caso di carichi crescenti monotonamente, il moltiplicatore di collasso, deﬁnito dai
teoremi dell’analisi limite [1]–[5], fornisce una valutazione eﬃcace della sicurezza della struttura. Procedure per l’analisi limite, generalmente di tipo path–following [6]–[7], sono ampiamente presenti nei codici FEM attualmente in commercio, e pertanto la veriﬁca risulta
relativamente agevole.
Tuttavia, le strutture sono soggette in generale ad azioni cicliche variabili deﬁnite solo
in termini di inviluppo complessivo, come nel caso delle azioni sismiche. In questi casi
l’analisi limite, anche considerando separatamente tutte le possibili combinazioni all’interno
dell’inviluppo, non fornisce da sola una valutazione aﬃdabile di sicurezza. Infatti il continuo
insorgere di deformazioni plastiche, conseguente alle variazioni cicliche dei carichi, potrebbe
portare alla perdita di funzionalit`
a o al crollo della struttura per degrado di resistenza,
per valori del moltiplicatore inferiori a quello teorico di collasso. E’ necessaria quindi
l’ulteriore richiesta progettuale che la formazione di deformazioni plastiche resti conﬁnata
ad un transitorio iniziale. Tale problematica costituisce la classica teoria dell’adattamento
o shakedown ed `e governata dai teoremi di Melan [8] e Koiter [9].
L’uso dei risultati teorici soﬀre, ancora oggi, della mancanza di algoritmi computazionali
eﬃcienti capaci di calcolare il moltiplicatore che assicura l’adattamento elastico di strutture reali. I metodi di soluzione proposti sembrano essere limitati al solo uso accademico
poich`e adattabili ad un problema speciﬁco. La maggior parte dei metodi numerici proposti valutano il moltiplicatore di adattamento direttamente dall’applicazione dei teoremi dell’adattamento (metodi diretti). In questi casi, l’analisi di adattamento `e formulata
come un problema di ottimizzazione vincolata che si presta ad essere risolto con i metodi
di programmazione matematica senza sfruttare le particolari caratteristiche del problema
strutturale. Ci`
o implica una bassa eﬃcienza computazionale del processo di soluzione e
impedisce che l’analisi di adattemento sia facilmente perfezionata per i codici di calcolo
strutturale commerciali.
Recentemente un nuovo metodo per la valutazione del fattore di sicurezza all’adattamento di strutture a comportamento elasto–plastico `e stato proposto da Casciaro–Garcea [10].
Tale metodo `e basato su una tecnica iterativa analoga all’algoritmo di tipo path–following
alla Riks, usato attualmente in analisi elasto–plastiche per valutare il percorso di equilibrio
di una struttura, ed oﬀre le stesse caratteristiche di robustezza, eﬃcienza e richiede un
analogo impegno computazionale. Nel presente lavoro di tesi saranno discussi gli aspetti
teorici e computazionali del metodo di analisi proposto in [10]. L’approccio `e esteso a due
3
Introduzione
4
nuove categorie di strutture: telai spaziali in cemento armato e strutture bidimensionali.
L’analisi a shakedown dei telai in cemento armato `e eﬀettuata mediante metodo FEM
il cui modello discreto di telaio tridimensionale `e ottenuto utilizzando l’elemento ﬁnito di
asta proposto da Petrolo–Casciaro [11]. Il comportamento elastoplastico delle sezioni `e
modellato mediante una rappresentazione lineare a falde del dominio limite per pressoﬂessione deviata, ottenuta in modo automatico a partire dalla geometria della sezione e dalle
armature presenti (Malena et al [12]). Per ogni sezione si assumono 26 falde, ciascuna delle
quali `e deﬁnita in funzione di un possibile cinematismo plastico della sezione. L’algoritmo `e
descritto in dettaglio e sono riportati e discussi una serie di risultati numerici che mostrano
l’eﬃcacia del metodo proposto come concreto strumento di veriﬁca strutturale.
L’analisi a shakedown di strutture bidimensionali piane caricate nel proprio piano `e realizzata mediante una formulazione ad elementi ﬁniti ottenuta con un elemento triangolare
misto che prevede tensioni costanti nell’area di inﬂuenza di ciascun nodo e spostamenti lineari all’interno dell’elemento. Tale elemento `e confrontato con elementi triangolari
compatibili a deformazione costante e deformazione lineare. Il comportamento elastoplastico del materiale `e descritto mediante una appropriata rappresentazione lineare del dominio
elastico, resa necessaria per ottenere un algoritmo di rientro applicabile in contesti del tutto
generici. L’algoritmo `e descritto in modo dettagliato `e sono riportati una serie di test che
ne dimostrano la robustezza e l’eﬃcienza.
La tesi `e organizzata come segue: nel capitolo 1 si riportano gli elementi della teoria
classica dell’adattamento, mentre nel capitolo 2 `e descritta una riformulazione della teoria
dell’adattamento appropriata per l’analisi numerica. Nel capitolo 3 `e decritto il metodo
iterativo di soluzione proposto in [10]; nel capitolo 4 si descrive il problema del rientro nel
dominio ammissibile ridotto che pu`
o essere realizzato in contesti del tutto generici. Inﬁne,
i capitoli 5 e 6 riportano le applicazioni fatte per i telai spaziali in c.a. e per le strutture
bidimensionali descrivendone il modello discreto utilizzato ed i risultati ottenuti.
Capitolo 1
La teoria classica dello shakedown
Per rendere meglio comprensibile l’argomento della tesi, in questo capitolo sar`
a trattata
brevemente la teoria elasto–plastica con particolare attenzione al problema di adattamento
plastico. Per trattazioni approfondite dell’argomento ci si pu`
o riferire alla presentazione
generale del problema fornita da Koiter [9] oppure al libro di K¨
onig [13].
1.1
Il legame costitutivo elastoplastico perfetto
Le ipotesi su cu si basa il legame elastico perfettamente plastico stabiliscono che il materiale
o subire incrementi all’aumentare
`e capace di raggiungere una tensione limite σ y che non pu`
dell’entit`
a delle azioni sollecitanti. La tensione σ[x] `e conﬁnata ad appartenere al dominio
ammissibile
E := {σ : f [σ] ≤ 0}
(1.1.1)
dove f [σ] (yield function) `e convessa in d , e tale che f [0] < 0. L’insieme E `e chiuso e
convesso in d e la sua frontiera ∂E `e caratterizata dalla condizione f [σ] = 0.
Per i corpi composti da questo tipo di materiale, ad una prima fase elastica della resistenza segue, al crescere delle azioni esterne sollecitanti, una seconda fase detta elastoplastica,
caratterizzata dalla presenza di zone dove il materiale, pervenuto alle condizioni di tensione
limite, diviene sede di deformazioni plastiche.
La natura irreversibile della deformazione plastica comporta che il legame possa essere
scritto solo in termini incrementali. Nell’ipotesi di piccoli spostamenti le deformazioni
elastiche e plastiche contribuiscono additivamente alla deformazione totale. Supponendo
che le propriet`
a elastiche non siano inﬂuenzate dalla storia di carico si pu`
o scrivere:
ε˙ = ε˙e + ε˙p
(1.1.2)
di cui, solo la componente elastica ε˙e `e linearmente dipendente dalla tensione attraverso il
legame elastico:
(1.1.3)
σ˙ = E ε˙e = E(ε˙ − ε˙p)
o
dove E `e il tensore elastico, simmetrico e deﬁnito positivo. La componente plastica ε˙p pu`
essere diversa da zero solo se σ appartiene a ∂E ed `e deﬁnita dalla legge di ﬂusso plastico:
0
if f < 0
p
(1.1.4)
ε˙ =
˙
µ˙ g[σ] , µ˙ ≥ 0 , µ˙ f = 0 if f = 0
5
La teoria classica dello shakedown
6
essendo il vettore g contenuto nel subdiﬀerenziale1 ∂f [σ] della funzione f in σ.
Le ipotesi di convessit`
a del dominio elastico e legge di ﬂusso plastico hanno anche una
formulazione di natura meccanica sul comportamento del materiale che si basa sul postulato
di Drucker, il quale si esprime con la diseguaglianza:
(σ y − σ)T ε˙p ≥ 0
∀σ ∈ E
(1.1.5)
in cui σ y e ε˙p sono legate dalla legge di ﬂusso (1.1.4). Per ogni valore di σ interna al
dominio elastico E la diseguaglianza (1.1.5) `e veriﬁcata in senso stretto. Inoltre, essendo
σ = 0 contenuto in E per deﬁnizione, si ha:
σTy ε˙p > 0
∀ε˙p = 0
(1.1.6)
` facile rendersi conto (v. ﬁg. 1.1) che il postulato di Drucker risulterebbe violato se la
E
deformazione plastica non fosse normale alla superﬁcie del dominio elastico o se il dominio
E non fosse convesso.
Figura 1.1: Condizioni che non veriﬁcano il postulato di Druker
1.2
Assunzioni sull’andamento dei carichi, l’inviluppo di sollecitazione elastica
Si assume che il carico esterno p[t] sia variabile nel tempo e sia espresso come una combinazione di azioni base pi , che possono variare in maniera arbitraria in funzione di un
parametro evolutivo che per comodit`a indicheremo con il tempo t, restando per`
o conﬁnate
nel dominio ammissibile che nel seguito indicheremo come dominio dei carichi
p
αi [t]pi : αmin
≤ αi [t] ≤ αmax
, ∀t
(1.2.1)
P := p[t] :=
i
i
i=1
1
Il subdiﬀerenziale della funzione f nel punto σ `e il cono deﬁnito dai vettori gradiente della funzione
calcolati in tutti i punti adiacenti a σ posti a distanza inﬁnitesima in tutte le direzioni. Nei punti regolari
il subgradiente coincide con il classico gradiente calcolato nel punto. Ulteriori dettagli di analisi convessa
possono essere tratti dal libro di Rockafellar [14]
La teoria classica dello shakedown
7
dove i fattori αi [t] sono contenuti in un poliedro in p. Tale modo di modellare i carichi `e
in eﬀetti conforme al modo di descrivere le azioni esterne in ambito ingegneristico per come
stabilito dai vari codici normativi europei. In eﬀetti tale descrizione riﬂette il fatto che
mentre `e incognita l’eﬀettiva evoluzione nel tempo dei carichi, in maniera statistica `e invece
possibile conoscere le escursioni minime e massime delle varie azioni base pi , o almeno si
puo pensare di conoscere il valore dell’azione per un preﬁssato periodo di ritorno in qualche
modo legato alla vita di esercizio della struttura che si sta progettando. Il dominio dei
carichi P deﬁnito nella forma (1.2.1) risulta essere un politopo convesso (v. ﬁg. 1.2).
pmin
= αmin
1
1 p1
min
p2 = αmin
2 p2
min p
pmin
=
α
3
3
3
pmax
= αmax
p1
1
1
max
max
p2 = α2 p2
pmax
= αmax
p3
3
3
Figura 1.2: Dominio dei carichi
Se si indica con σei la soluzione elastica in tensione equilibrata con le singole azioni base
pi `e possibile deﬁnire il dominio Se delle tensioni elastiche σ e [t] associate a p[t] come segue
p
αi [t]σei : αmin
≤ αi [t] ≤ αmax
, ∀t
(1.2.2)
Se := σ e [t] :=
i
i
i=1
L’insieme Se rappresenta l’inviluppo delle tensioni elastiche. Esso rappresenta i valori che
punto per punto del corpo vengono prodotti ai vari istanti di tempo e pertanto durante
la vita della struttura per ogni carico p[t] contenuto in P. Ovviamente Se `e anche chiuso
e convesso. In particolare nello spazio di deﬁnizione delle tensioni Rd esso `e un politopo
convesso, cio`e un poliedro chiuso e limitato.
Nella trattazione presente la variabile temporale t `e assunta come una variabile evolutiva
mentre nel corso della trattazione considereremo sempre trascurabili gli efetti dinamici
dovuti ai carici bench`e la trattazione possa essere facilmete generalizzata alla presenza di
eﬀetti dinamici (vedi Ceradini [15, 16] e Polizzotto [17, 18]).
Il problema dello shakedown pu`
o, in maniera molto sempliﬁce, essere aﬀermato come
la condizione per cui una struttura, soggetta a carichi varibili nel tempo e limitati solo da
un dominio di ammissibilit`
a possa, a parte un transitorio iniziale in cui il comportamento
pu`
o essere di tipo elastoplastico e quindi comportare anche l’insorgere di incrementi di
deformazione plastica, ritornare ad essere puramente elastico.
Sia λc il moltiplicatore di collasso della struttura, la condizione λ ≤ λc non implica
l’assenza di deformazioni plastiche, ﬁnch`e le stesse sono ancora contenute all’interno delle
La teoria classica dello shakedown
8
restanti zone elastiche della struttura e siano pertanto tali da non poter realizzare un
meccanismo di collasso. Le deformazioni plastiche, sotto l’azione di carichi ripetuti possono
per`
o sommarsi e quindi portare al collasso la struttura per perdita di funzionalit`
a. Non
meno grave e problematico `e il caso in cui le deformazioni, invece di sommarsi all’interno di
cicli di carico e scarico si elidono a vicenda sommandosi con segni opposti. L’evidenza ﬁsica
mostra infatti come in questi casi il modello matematico utilizzato di materiale elasto–
plastico perfetto sia assolutamente inadeguato a descrivere la realt`a del comportamento
strutturale, che presenta invece, all’aumentare dei cicli di carico e scarico cui e sottoposto,
un decadimento delle caratteristiche di resistenza e pertanto nella nostra schematizzazione
matematica una contrazione del dominio elastico del materiale stesso. Il fenomeno di cui
abbiamo appena parlato prende comunemente il nome di “fatica” ed `e molto pericoloso
perch`e comporta un abbattimento del valore del moltiplicatore di collasso λc e quindi la
possibilit`
a di collasso per valori del carico inferiori a quelli teorici.
La problematica precedentemente evidenziata mette in luce come una ulteriore richiesta
necessaria per poter operare con un certo margine di sicurezza `e che la struttura in qualche
modo si adatti elasticamente. In altri termini si vuole progettare la struttura in modo che,
dato il carico esterno p[t] esista un istante di tempo t¯ tale che per ogni t ≥ ¯t la struttura
si comporti in maniera elastica o, per usare la terminologia tecnica, si adatti elasticamente.
Ci`
o ovviamente non esclude che ci siano deformazioni plastiche accumulate ma esclude che
ci siano incrementi di deformazioni plastiche per istanti di tempo t ≥ t¯ e pertanto che in
qualche modo la deformazione plastica risulti limitata. In termini pi`
u formali ci`
o signiﬁca
che, per una struttura soggetta ad un sistema di azioni variabili nel tempo, superato un
periodo transitorio iniziale in cui `e possibile l’insorgere di deformazioni plastiche irreversibili,
la risposta della struttura in tensioni debba essere del tipo
¯
σ[t] = σe [t] + σ
f [σ[t]] < 0
∀t > t¯
(1.2.3)
¯ lo
con σ e [t] termine che si avrebbe se la riposta della struttura fosse puramente elastica e σ
stato di coazione indotto dalle deformazioni plastiche che non varia con il tempo per istanti
successivi a t¯. Ovviamente il termine elastico `e, per deﬁnizione, equilibrato con i carichi
esterni, ne discende che il contributo in σ
¯ debba essere autoequilibrato e quindi rappresenta
uno stato di autotensione. L’equazione precedente ci dice che, per una struttura soggetta
ad una storia di carico variabile e dipendente dal tempo, una volta superato un transitorio
iniziale identiﬁcato dalla condizione t < t¯, il termine in autotensioni non dipende piu dal
tempo e pertanto le variazioni in tensione sono legate al solo contributo elastico.
1.3
Definizione di shakedown elastico
Si dice che, sotto azioni esterne variabili e ripetute nel tempo, una struttura si adatta in
campo elastico se ci`o avviene per ogni possibile carico contenuto nel dominio dei carichi P.
Ne scaturisce la seguente deﬁnizione:
Deﬁnizione 1.3.1 (Shakedown). Si dir`
a che una struttura si adatta ad uno stato elastico, o semplicemente che avviene l’adattamento se, dopo una fase iniziale durante la quale `e
possibile il formarsi di deformazioni plastiche, la risposta strutturale, per ogni carico contenuto nel dominio di ammissibilit`
a p[t] ∈ P, tende ad essere puramente elastica, cio`e la
La teoria classica dello shakedown
9
dissipazione di lavoro plastico si esaurisce. Si avr`
a:
∞ T p
σ[t] ε˙ [t] dv dt < ∞
t=0
(1.3.1)
B
essendo σ[t] e εp [t] la tensione e la deformazione plastica prodotte durante il processo di
carico. Il punto indica la derivazione rispetto al tempo, mentre t = 0 indica lo stato iniziale
vergine caratterizzato da εp[0] = 0.
Una ovvia conseguenza di questa deﬁnizione `e che la precedente condizione di shakedown
¯ indipendente dal
o di adattamento implica l’esistenza di un campo di autotensioni σ
¯ ∈S
tempo tale che
¯] < 0
∀σ e [t] ∈ Se
(1.3.2)
f [σ e [t] + σ
Ci`o `e facilmente dimostrabile. Infatti, per la positivit`
a della quantit`
a sotto segno di integrale, conseguente alla condizione di stabilit`
a di Drucker eq.(1.1.6), dovr`a necessariamente
essere
˙ − ε˙e [t]) = 0 , ∀p[t] ∈ P
lim ε˙p [t] ≡ lim (ε[t]
t→∞
t→∞
Per l’unicit`
a della soluzione incrementale elastica si ha:
σ[t]
˙ = E(ε[t]
˙ − ε˙p [t])
cio`e
ε˙p [t] = ε[t]
˙ − E −1 σ[t]
˙ = E −1 (E ε[t]
˙ − σ[t])
˙
⇒
˙
ε˙p[t] = E −1 (σ˙ e [t] − σ[t])
Si ottiene pertanto
˙ − σ˙ e [t]) = 0
lim (σ[t]
t→∞
⇒
lim σ[t] = σ e [t] + σ
¯
t→∞
,
∀σ e [t] ∈ Se
Inoltre, poich`e σ[t] e σe [t] sono equilibrati dallo stesso carico p[t], per deﬁnizione la loro
diﬀerenza σ
¯ `e un campo di autotensione. Si osservi che la forza del risultato sta proprio nel
fatto che la condizione deve valere per ogni tensione elastica contenuta nel dominio elastico
e quindi, nella nostra formulazione, per ogni storia di carico possibile per la struttura e per
ogni t.
1.4
Il teorema statico dell’adattamento
Il teorema statico fu dimostrato da Bleich–Melan, [8] [19], nel 1936 `e si enuncia nel modo
seguente:
Teorema 1.4.1 (Teorema statico di adattamento). Condizione necessaria e suﬃciente
aﬃnch`e una struttura soggetta a carichi esterni variabili nel tempo o ripetuti si adatti in
campo elastico `e che esista uno stato di autotensioni equilibrato σ
¯ , non dipendente dal
tempo, tale che:
¯
σ ∗ [t] = σ e [t] + σ
,
f [σ ∗ [t]] < 0
,
∀σ e [t] ∈ Se , t > ¯t
(1.4.1)
avendo indicato con σ e [t] il campo tensionale corrispondente alle azioni esterne nel caso
di comportamento puramente elastico e con σ ∗ [t] una soluzione ﬁttizia che indicheremo
come nominale, corrispondente alla soluzione elastica σ e [t] ad eccezione di uno stato di
autotensioni σ
¯ indipendente dal tempo e tale che σ∗ [t] sia all’interno del dominio elastico
del corpo.
La teoria classica dello shakedown
10
Dimostrazione. Che la condizione sia necessaria `e evidente. Infatti se si veriﬁca l’adattamento deve esistere lo stato di coazione deﬁnito nell’enunciato.
Che la condizione sia suﬃciente si dimostra nel seguente modo. Siano σ[t] e ε[t] le
tensioni e deformazioni nella struttura per eﬀetto dei carichi p[t] variabili nel tempo:
εe [t] = F σ e [t]
¯
σ[t] = σ e [t] + σ[t]
,
¯[t] + εp [t]
ε[t] = εe [t] + ε
¯[t] = F σ[t]
¯
ε
¯ un campo di tensioni auto–equilibrate e con ε
¯[t] la deformazione
in cui si `e indicato con σ[t]
associata. In particolare la diﬀerenza tra la tensione reale e quella ﬁttizia sar`
a il campo di
autotensioni
¯ −σ
σ[t] − σ ∗ [t] = σ[t]
¯
avendo sfruttato il fatto che la soluzione elastica `e ovviamente la medesima per i due
processi. Si introduca la quantit`
a, deﬁnita positiva:
1
(σ − σ ∗ )TF (σ − σ ∗ ) dv ≥ 0
Ψ[t] :=
2
B
per derivazione rispetto al tempo si ottiene:
d
˙ =
¯ − σ)
¯ dv = (σ − σ ∗ )TF σ
(σ − σ ∗ )TF (σ[t]
¯˙ [t] dv = (σ − σ ∗ )T¯ε˙ [t] dv
Ψ[t]
dt
B
B
B
avendo fatto uso dell’equazione di legame elastico che pone ¯ε˙ [t] = F σ
¯˙ [t]. Ricordando che
˙ − ε˙ e [t] − ε˙ p [t]
¯ε˙ [t] = ε[t]
si pu`
o ancora porre:
∗
˙ =
Ψ[t]
(σ − σ ∗ )Tε˙ p [t] dv
e
˙ − ε˙ [t]) dv −
(σ − σ ) (ε[t]
T
B
B
e
˙ e ε˙ [t] sono deformazioni cinematicamente compatibili. Il primo integrale rapprein cui ε[t]
senta il lavoro di un campo di autotensioni per un campo di deformazioni cinematicamente
compatibili per cui `e nullo2 . Il termine restante
˙ < 0 se ε˙p = 0)
˙ = − (σ − σ ∗ )T ε˙ p dv ≤ 0
(Ψ
Ψ[t]
B
`e sempre minore di zero se la deformazione plastica incrementale ε˙p = 0. In tal caso, infatti,
l’instaurarsi di deformazioni plastiche incrementali comporta che σ[t] sia sulla frontiera del
dominio elastico del materiale, essendo la tensione associata dalla legge di normalit`
a alla
deformazione incrementale. Ne segue che, per il postulato di Drucker, il termine sotto
2
¯ un campo di autotensioni, cio`e equilibrato a carichi nulli (f = 0, b = 0) e sia ε˙ un campo di
Sia σ
deformazioni cinematicamente compatibile tale da rispettare le condizioni al contorno rese omogene (u˙ = 0).
Il princio dei lavori virtuali fornisce:
¯ T εdV
˙
˙
˙
¯ T εdV
˙
bT udV
f T udA
σ
−
−
=0
⇒
σ
=0
B
B
∂Bf
B
La teoria classica dello shakedown
11
l’integrale `e strettamente maggiore di zero ogni qualvolta l’incremento di deformazione
plastica `e diverso da zero. Ne deriva che Ψ[t] `e una funzione non–negativa (quindi limitata
inferiormente) e sempre descrescente durante il processo plastico poich`e ha una derivata
prima negativa. Ci`
o implica necessariamente:
˙ =0
lim Ψ[t]
t→∞
⇒
lim ε˙p = 0
t→∞
(1.4.2)
Il processo al limite utilizzato, comunque, non da informazioni sul tempo t¯ di adattamento
della struttura che potrebbe essere imprecisabilmente elevato.
In eﬀetti la dimostrazione rigorosa di shakedown, che richiede appunto che la condizione di ammissibilit`a sia soddisfatta in senso stretto, comporta la possibilit`
a di deﬁnire
un numero suﬃcientemente piccolo ν > 0 tale che
f [(1 + ν)σ ∗ [t]] ≤ 0 ,
ν>0
Il postulato di Drucker si scrive come:
(σ[t] − (1 + ν)σ ∗ [t])Tε˙ p [t] ≥ 0
σ[t]Tε˙ p [t] ≥ (1 + ν)σ ∗ [t]Tε˙ p
⇒
sfruttando la relazione
σ[t]Tε˙ p [t] = (1 + ν)σ[t]T ε˙ p [t] − νσ[t]T ε˙ p[t]
ed integrando sull’intero volume, si ottiene:
T p
˙
σ[t] ε˙ [t] dv ≤ (1 + ν) (σ[t] − σ ∗ [t])Tε˙p [t] dv ≡ −(1 + ν)Ψ[t]
ν
B
B
Ricordando la condizione di adattamento (1.3.1) si ottiene la richiesta condizione di energia
plastica dissipata ﬁnita:
∞ 1+ν
1+ν
1+ν
T p
(Ψ[0]−Ψ[∞]) ≤
Ψ[0] =
σ[t] ε˙ [t] dv dt ≤
σ
¯ T E −1 σ
¯ dv < ∞
ν
ν
2ν B
t=0
B
¯[0] = 0 si ha σ[0] = σe [0], mentre Ψ[∞] ≥ 0 per deﬁnizione. La condizione di
poich`e ε
adattamento risulta soddisfatta.
1.5
Il teorema cinematico dell’adattamento
Il teorema cinematico dell’adattamento, dimostrato per la prima volta da Koiter [20],
deﬁnisce una condizione suﬃciente di mancato adattamento. Anche in questo caso sar`a
conveniente in qualche modo svincolarsi dal legame con la storia di carico.
Sia ε˙ p un campo di incrementi di deformazioni plastiche non necessariamente congruenti nei singoli istanti in cui si veriﬁcano e σ y [t] le tensioni associate alle deformazioni ε˙ p [t]
attraverso la legge di ﬂusso (1.1.4). Diremo che tale distribuzione costituisce un ciclo ammissibile se il suo integrale in un intervallo di tempo [t1 , t2 ] `e cinematicamente compatibile.
Sia εp [t] tale integrale si avr`a:
t2
p
p
(1.5.1)
ε˙ p [t]dt
ε = ε [t] :=
t1
La teoria classica dello shakedown
12
⎧
⎨ εp = 1 ∇u + ∇uT 2
⎩u = 0
tale che
in
B
in ∂Bu
avendo indicato con B il volume del solido in esame e con ∂Bu la porzione di frontiera su
cui sono assegnati gli spostamenti. Si dimostrer`
a il seguente teorema:
Teorema 1.5.1 (Condizione di non adattamento). Condizione suﬃciente per il mancato adattamento `e che esista un ciclo ammissibile ε˙ p per cui sia valida la seguente condizione
t2 (σe [t] − σ y [t])T ε˙ p [t]dv dt > 0
(1.5.2)
∃σ e ∈ Se :
t1
B
Dimostrazione. Il teorema si dimostra per assurdo. Assumendo che la struttura presenti
¯ tale che:
shakedown, per la (1.3.2), esister`a un campo di autotensioni σ
¯ ≤0
f [σ e [t] + σ]
∀σ e ∈ Se
¯ `e un campo di
Poich`e la tensione σy [t] `e associata alla εp [t] dalla legge di ﬂusso e σ e [t] + σ
tensioni staticamente ammissibile, `e possibile scrivere, utilizzando il postulato di Drucker:
¯ Tε˙ p [t] ≥ 0
(σy [t] − σ e [t] − σ)
∀σ e ∈ Se
Integrando l’intero ciclo sul volume del corpo B, si ottiene:
t2
t2 T p
T
(σ y [t] − σe [t]) ε˙ [t] dv dt −
σ
¯
ε˙p [t] dt dv ≥ 0
t1
B
B
t1
Il secondo integrale `e nullo poich`e σ
¯ `e una tensione autoequilibrata e εp `e una deformazione
cinematicamente compatibile, per cui diventa:
t2 (σ e [t] − σy [t])T ε˙p [t] dv dt ≤ 0
∀σ e ∈ Se
t1
B
che pertanto contraddice l’assunto del teorema.
1.6
Il fattore di sicurezza a shakedown
La sicurezza rispetto allo shakedown pu`o essere valutata come massimo moltiplicatore del
dominio dei carichi P oppure, in forma equivalente del dominio delle tensioni elastiche Se ,
compatibile con le condizioni di shakedown.
Possiamo deﬁnire questo valore, che chiameremo fattore di sicurezza a shakedown λa,
riferendoci al dominio ampliﬁcato λσe , come l’estremo superiore dei moltiplicatori “strettamente sicuri” λs (safe) che soddisfano i requisiti del teorema statico, cio`e
¯ : f [λsσ e + σ
¯] < 0
∃¯
σ∈S
∀σ e ∈ Se
(1.6.1)
La teoria classica dello shakedown
13
Analogamente, possiamo deﬁnire λa come l’estremo inferiore dei moltiplicatori “strettamente insicuri” λu (unsafe) che soddisfano i requisiti del teorema cinematico, cio`e
t2
t2 p
p
ε˙ [t] ∈ K ∃σ e [t] ∈ Se :
(λu σ e [t] − σ y [t])Tε˙p [t] dv dt > 0
(1.6.2)
∃ε :=
t1
t1
B
essendo σ y [t] e ε˙p [t] associati con la legge di ﬂusso (1.1.4). Si dimostra che
λa = sup λs = inf λu
(1.6.3)
Dimostrazione. Assumiamo λa deﬁnito dalla condizione
¯] = 0
σ e ∈ Se , x ∈ B
λa : min max f [λaσe + σ
σ
¯ ∈¯S σ e , x
che signiﬁca riferirsi al punto x ∈ B per cui si veriﬁca che la funzione f [λσe + σ
¯ ] assume il
¯
valore massimo e scegliere tra i possibili campi di autotensioni σ
¯ ∈ S quello minimo capace
¯.
di ripristinare l’ammissibilit`
a della tensione λaσ e + σ
Per deﬁnizione di funzione convessa si ha:
f [cσ] < cf [σ] + (1 − c)f [0] ≤ 0
e l’assunzione f [0] < 0 implica che per 0 ≤ c < 1 si veriﬁca semplicemente:
f [cσ] < cf [σ] ≤ f [σ]
se
f [σ] ≤ 0
Sfruttando la propriet`
a appena introdotta per la funzione f [σ] si ha:
f[
λs
λs
(λaσ e + σ
¯ )] <
f [λaσ e + σ
¯ ] < f [λaσ e + σ
¯]
λa
λa
se
λs < λa
poich`e si pu`
o scrivere:
¯ ) = (λs σe +
c(λaσ e + σ
λs
σ
¯) =
λa
dove
c=
λs
λa
con c < 1 per λs < λa. Si ha, quindi:
λs
¯ ] < min f [λaσ e + σ
¯] = 0
min f [λs σ e + σ
λa
σ
¯ ∈¯S
σ
¯ ∈¯S
∀σ e ∈ Se
se
¯ s < λa
λ
¯ s < λa `e strettamente sicuro secondo la deﬁnizione (1.6.1). Inoltre, per la (1.1.5),
Pertanto, λ
si ha
¯ )Tε˙p [t]
min max (σy [t] − λaσ e [t] − σ
σ
¯ ∈ S¯ σ e ∈ Se
≥0
∀ε˙p
dove σ y [t] e ε˙p [t] sono associati dalla legge di ﬂusso (1.1.4) e il segno di eguale `e raggiunto
¯ ∗ e ε˙p∗ [t], per deﬁnizione. Tale soluzione `e caratterizzata dalla
per σ ∗y [t] := λaσ ∗e [t] + σ
condizione di minimo:
t2 (σ y [t] − λaσ ∗e [t] − σ
¯ )Tε˙p [t] dv dt = 0
min
¯ ∈¯S t1 B
ε˙p , σ
La teoria classica dello shakedown
14
¯ implica:
ottenuta integrando sull’intero ciclo [t1 , t2 ]. Il minimo rispetto a σ
t2
T
¯
δσ
¯
ε˙p∗ [t] dt dv = 0 ∀δ σ
¯∈S
B
t1
ci`o implica che εp∗ debba essere una deformazione cinematicamente compatibile. Inoltre
per il postulato di Drucker (1.1.6), si ha
t2 t2 ∗ T p∗
σ e [t] ε˙ [t] dv dt =
σ ∗y [t]Tε˙p∗ [t] dv dt > 0
λa
t1
e di conseguenza
t2
t1
B
B
t1
B
(λu σ ∗e [t] − σ ∗y [t])Tε˙p∗ [t] dv dt > 0
se
λu > λa
¯ u > λa `e strettamente insicuro secondo la deﬁnizione (1.6.2).
Pertanto λ
Riferendoci, quindi, ai moltiplicatori strettamente sicuri λs e strettamente insicuri λu ,
il moltiplicatore di adattamento `e deﬁnito semplicemente dalla relazione:
λa := max λs = min λu
1.7
(1.6.4)
Alcuni commenti sul moltiplicatore di sicurezza
Per quanto visto nei paragraﬁ precedenti, si noti come il teorema statico dello shakedown `e
una generalizzazione del teorema statico dell’analisi limite per una combinazione dei carichi
esterni che variano in un dominio ammissibile P. Infatti, il teorema statico dell’analisi
limite `e la forma particolare del teorema statico dello shakedown ottenuta per una singola
condizione di carico, di intensit`
a unitaria e moltiplicatore proporzionale al fattore di carico λ
(p = 1 , 0 ≤ α1 ≤ λ). Analogamente, il teorema cinematico dello shakedown corrisponde ad
una generalizzazione del teorema cinematico dell’analisi limite e quest’ultimo corrisponde
alla sua specializzazione per singola condizione di carico. Come nell’analisi limite, uno stato
autoequilibrato di pretensioni non ha nessuna inﬂuenza sul moltiplicatore di adattamento.
Inoltre, il teorema statico dello shakedown assicura che il moltiplicatore al limite elastico
λe = max λ : f [λσe ] ≤ 0
,
∀σ e ∈ Se
(1.7.1)
`e minore (o quanto meno non maggiore) del moltiplicatore di shakedown λa . Infatti, ad¯ = 0, si soddisfano i
dizionando alla tensione elastica λσe [t] la tensione autoequilibrata σ
requisiti del teorema statico. D’altra parte, il teorema cinematico dello shakedown assicura
che, per ciasuna condizione di carico p ∈ P, il moltiplicatore di collasso
(σy − λσe )T ε˙p dv = 0
∀σ e ∈ Se , ∀ε˙p ∈ K
(1.7.2)
λc = min λ :
B
`e maggiore (o quanto meno non minore) di λa. Il moltiplicatore di shakedown `e quindi
limitato dalle condizioni:
(1.7.3)
λe ≤ λa ≤ λc
La teoria classica dello shakedown
15
L’osservazione precedente fornisce chiare motivazioni per l’uso del moltiplicatore al limite
elastico nella pratica della progettazione strutturale. Infatti, cos`ı facendo `e possibile evitare
l’analisi a shakedown pur soddisfacendone implicitamente i requisiti. Tuttavia, si deve tener
presente che una esplicita analisi a shakedown diventa necessaria se la progettazione `e basata
su metodi nonlineari, essendo la semplice determinazione del moltiplicatore di collasso non
sicura.
Capitolo 2
Riformulazione del problema in
termini della funzione di
shakedown
In questo capitolo `e presentata una riformulazione del problema dello shakedown adatta
per una implementazione di tipo FEM, ottenendo anche dei risultati preliminari utili per
la discussione delle propriet`
a di convergenza del metodo proposto.
2.1
La funzione di snervamento a shakedown e il dominio
elastico ridotto
` conveniente introdurre la funzione di snervamento a shakedown
E
fs [σ, λ] := max {f [λσe + σ]}
(2.1.1)
σe ∈ Se
ed il relativo dominio ammissibile a shakedown
Es [λ] := {σ : fs [σ, λ] ≤ 0}
(2.1.2)
che rappresenta l’insieme di tutte le possibili traslazioni σ del dominio λSe all’interno di E
a vuoto per valori positivi
(v. ﬁg. 2.1). Poich`e si `e assunto f [0] < 0, il dominio Es [λ] non sar`
suﬃcientemente piccoli di λ. Inoltre, essendo S e E entrambi domini chiusi e convessi,
¯ il valore:
risulter`a anche Es [λ] un dominio chiuso e convesso. Indichiamo con λ
¯ := sup{λ : Es [λ] = ∅} .
λ
(2.1.3)
risulta:
Es [λ1 ] = ∅ ,
¯
λ2 < λ1 < λ
αEs [λ1 ] ⊂ Es [λ2 ] se
,
α :=
λ2
<1
λ1
(2.1.4)
La deﬁnizione (2.1.2) implica l’equivalenza
(λσe + σ) ∈ E
⇐⇒
16
σ ∈ Es [λ] .
(2.1.5)
Riformulazione in termini del dominio elastico e della funzione di shakedown
17
Figura 2.1: Domini ammissibili E, Es [λ]
¯ e σ y ∈ ∂Es [λ], esister`a una o pi`
Pertanto, per λ ≤ λ
u tensioni σ yk = λσ ek + σ y ∈ ∂E che
possono essere associate σ y , come mostrato in ﬁgura 2.1. Per la condizione (1.1.5) si ha:
(σyk − σ k )Tε˙pk ≥ 0 ∀σ ∈ Es [λ]
(2.1.6)
(σy − σ)T εp =
k
in cui σ k := λσek + σ ∈ E,
deﬁnito dalla combinazione
ε˙pk
sono le deformazioni plastiche associate alle σ yk , e εp ,
εp :=
ε˙ pk
(2.1.7)
k
`e l’incremento di deformazione plastica associato a σ y . La condizione (2.1.6) implica la
convessit`a della funzione fs [σ, λ] e la legge di normalit`a tra σ y e εp :
εp = µ g
2.2
,
g ∈ ∂fs [σ y ; λ] ,
µ≥0
,
µfs = µf˙s = 0
(2.1.8)
Il fattore di sicurezza a shakedown
Utilizzando le deﬁnizioni introdotte nel paragrafo precedente si pu`
o derivare una caratterizzazione semplice del fattore di sicurezza a shakedown. Dal teorema statico, ricordando
la (2.1.1), si ottiene:
¯ : fs [¯
¯∈S
σ , λs ] ≤ 0
(2.2.1)
λs ≤ λa se ∃ σ
Riformulazione in termini del dominio elastico e della funzione di shakedown
18
in cui λs `e un moltiplicatore “strettamente
p sicuro”. Dal teorema cinematico, ricordando le
p
assunzioni σ y = σ yk − λσek e ε = k ε˙k , si ottiene:
p
σ Ty εp dv ≤ 0 , εp = 0
(2.2.2)
λu ≥ λa se ∃ ε ∈ K :
B
con λu moltiplicatore “strettamente insicuro”.
Il fattore di sicurezza a shakedown pu`
o essere deﬁnito come l’estremo superiore dei moltiplicatori “strettamente sicuri” secondo la (2.2.1) e come estremo inferiore dei moltiplicatori
“strettamente insicuri” secondo (2.2.2).
2.3
Il problema di shakedown
Tenendo conto delle deﬁnizioni introdotte nel paragrafo precedente, il problema dello shakedown pu`
o essere deﬁnito come:
Deﬁnizione 2.3.1 (Il problema dello shakedown). Determinare il fattore di sicurezza
allo shakedown
(2.3.1)
λa := max λs = min λu
in cui λs e λu sono moltiplicatori “strettamente sicuri” e “strettamente insicuri” secondo
la (2.2.1) e (2.2.2).
2.4
Lo schema return mapping per la funzione di shakedown
Si assume che la struttura sia stata discretizzata ad elementi ﬁniti, e sia u ∈ N il vettore
che raccoglie i parametri liberi di spostamento, u0 il vettore noto degli spostamenti iniziali,
a D[x]
ε[u] := D(u − u0 ) le deformazioni associate attraverso la matrice di compatibilit`
all’incremento di spostamento u − u0 , σ 0 la tensione nello stato iniziale. Partendo da un
predittore elastico σ∗ = σ0 + Eε[u] non necessariamente contenuto in Es [λ], la tensione
¯ si ottiene dal seguente algoritmo di rientro
ammissibile σa [ε, λ] corrispondente a u e λ ≤ λ
nel dominio ammissibile:
σ a [ε, λ] ≡ σ ∗ − Eεp
,
σa [ε, λ] ∈ Es [λ]
(2.4.1)
che ci permette di costruire una tensione associata σ a , contenuta in Es [λ], corrispondente
alla deformazione εp attraverso la legge di ﬂusso (2.1.8):
=0
if fs [σ ∗ , λ] < 0
(2.4.2)
εp = µg , g ∈ ∂fs [σ ; λ] , µ :
≥ 0 , fs [σ, λ] = 0 if fs [σ ∗ , λ] ≥ 0
Lo schema (2.4.1-2.4.2) corrisponde al noto “return mapping by closest–point projection”
che pu`
o essere ottenuto anche minimizzando la funzione di Haar–Karman sotto particolari
vincoli. Tale descrizione, che si presta ad essere applicata in ambito FEM, dove le leggi
costitutive sono deﬁnite localmente, `e ottenuta per ciascun punto di Gauss attraverso la
condizione:
1
∗
(σ − σ ∗ )TE −1 (σ − σ ∗ ) dv = min , σ = σ[ε, λ] ∈ Es [λ] (2.4.3)
φ[σ − σ ] :=
2 Bg
(σ )
Riformulazione in termini del dominio elastico e della funzione di shakedown
19
in cui Bg `e il volume del punto di Gauss (v. ﬁg. 2.2). Ovviamente, deﬁnita come minimo
di una funzione strettamente convessa su un dominio convesso, σ a `e funzione univoca di ε
e λ.
¯
Figura 2.2: Dominio elastico E[λ]
per λ2 > λ1
2.5
Propiet`
a notevoli dello schema
In questa sezione alcuni lemmi illustreranno le propriet`
a caratteristiche dello schema di
rientro (2.4.1-2.4.2). Ci`
o render`
a meglio comprensibile l’algoritmo di soluzione che sar`a
proposto nel prossimo capitolo e di cui saranno dimostrate le propriet`
a di convergenza.
Lemma 2.5.1. La funzione σ a [ε, λ] ha derivate direzionali rispetto a ε ed a λ. Esiste
inoltre una funzione potenziale convessa Ψ[ε, λ] tale che
σ a [ε, λ] =
per cui l’operatore tangente
E t[ε, λ] :=
∂Ψ[ε, λ]
∂ε
∂2Ψ
∂σ a
=
∂ε
∂ε2
(2.5.1)
`e autoaggiunto: E t = E Tt .
Dimostrazione. Questo `e un classico risultato della teoria della plasticit`
a incrementale ([21],
[22]) e deriva direttamente dalla deﬁnizione di σ a[ε, λ] attraverso la condizione di Haar–
K´
arm´an. Infatti, per le assunzioni fatte su f [σ] e la deﬁnizione (2.1.1), fs [σ, λ] `e una
trasformazione limitata, univoca e continua d+1 → , pertanto `e direzionalmente difˆ con t ≥ 0, risultano deﬁnite le
ferenziabile. Quindi, per ogni percorso {ε + tˆ
ε , λ + tλ},
ˆ e E t [ε, λ; ε
ˆ] tale che
derivate direzionali σ t [ε, λ; λ]
dσ[ε, λ] = E t dε + σ t dλ
,
σ t :=
∂σ a [ε, λ]
∂λ
Riformulazione in termini del dominio elastico e della funzione di shakedown
20
L’esistenza del potenziale Ψ `e assicurata dalla condizione
σ a [ε, λ]Tdε = 0
C
su ogni curva chiusa C nello spazio delle ε. Ricordando la (2.4.1), si pu`
o porre dε =
E −1 dσ a + dεp che consente di scrivere
1 T −1
T
T
−1
T
p
T p
σ E σ a + σ a ε − dσ Ta εp
σ a dε = σ a E dσ a + σ a dε = d
2 a
che implica
C
σ a dε +
T
T
C
p
dσ a ε =
d
C
C
=0
da cui
1 T −1
σ E σ a + σTa εp
2 a
σTa dε = −
C
dσTa εp
Per la legge di ﬂusso (2.4.2), si ha εp = 0 se σa `e interna a Es [λ]; si ha dσ a = 0 se σ a
corrisponde ad uno spigolo di ∂Es [λ]e εp `e interna al cono delle normali; si ha dσ Ta εp = 0
negli altri casi. Pertanto si ottiene c dσ Ta εp ≡ 0, con ci`o provando la seconda parte del
lemma.
Lemma 2.5.2. Sia {σ 0 , ε0 } uno stato iniziale e {σ1 , ε1 } e {σ 2 , ε2 } due stati diﬀerenti
ottenuti dall’algoritmo di rientro (2.4.3):
σ 1 = σ a [ε1 , λ] = σ 0 + E(ε1 − εp1 )
(2.5.2)
σ 2 = σ a [ε2 , λ] = σ 0 + E(ε2 − εp2 )
in cui εp1 e εp2 sono associati a σ 1 e σ 2 dalla legge di ﬂusso (2.1.8). Valgono le seguenti
aﬀermazioni:
0 ≤ (σ 2 − σ 1 )T(ε2 − ε1 ) ≤ (ε2 − ε1 )T E(ε2 − ε1 )
e il veriﬁcarsi di (σ 2 − σ 1 )T(ε2 − ε1 ) = 0 implica
σ1 = σ2
Dimostrazione. La condizione di Druker (1.1.5) fornisce
p
(σ2 − σ1 )T ε2 ≥ 0
quindi, addizionando:
,
p
(σ 1 − σ 2 )Tε1 ≥ 0
(σ 2 − σ 1 )T(εp2 − εp1 ) ≥ 0
Dalla (2.5.2) si ottiene inoltre
(σ 2 − σ 1 ) = E(ε2 − ε1 ) − E(εp2 − εp1 )
Otteniamo quindi:
(σ 2 − σ 1 )T (ε2 − ε1 ) = (ε2 − ε1 )TE(ε2 − ε1 ) − (εp2 − εp1 )T E(εp2 − εp1 ) − (σ 2 − σ 1 )T(εp2 − εp1 )
≤ (ε2 − ε1 )TE(ε2 − ε1 ) − (εp2 − εp1 )T E(εp2 − εp1 ) ≤ (ε2 − ε1 )TE(ε2 − ε1 )
(σ 2 − σ 1 )T(ε2 − ε1 ) = (σ 2 − σ 1 )TE −1 (σ 2 − σ 1 ) + (σ 2 − σ 1 )T (εp2 − εp1 )
≥ (σ 2 − σ 1 )TE −1 (σ 2 − σ 1 )
Che prova il lemma.
Riformulazione in termini del dominio elastico e della funzione di shakedown
21
Lemma 2.5.3. La matrice tangente E t [ε, λ] soddisfa le condizioni
0 ≤ εT E t ε ≤ εT Eε
∀ε
,
Dimostrazione. La dimostrazione deriva direttamente dal lemma 2.5.2 se si pone ε := ε2 −ε1
per ε2 → ε1 .
Lemma 2.5.4. Sia σ1 := σ a [ε1 , λ1 ] e σ 2 := σ a [ε2 , λ2 ], possiamo scrivere la loro diﬀerenza
nella forma
σ 2 − σ 1 = σ s [ε1 , λ1, λ2 ](λ2 − λ1 ) + E s [λ2 , ε1 , ε2 ](ε2 − ε1 )
Inoltre, risulta
E s = E Ts
,
0 ≤ εT E s ε ≤ εT E ε
e
σ s [ε1 , λ1 , λ2]T εp < 0
p
,
se εp :=
,
∀ε
εp11 if λ1 ≤ λ2
p
ε12 if λ2 < λ1
= 0
p
in cui ε11 e ε12 sono rispettivamente gli incrementi di deformazione plastica associati a
σ a [ε1 , λ1] e σ a [ε1 , λ2 ].
Dimostrazione. La prima parte del lemma si ottiene direttamente deﬁnendo E s and σs
attraverso i rapporti secanti
1
1
∂σ[ε1 , λ[t]]
dt =
σ t [ε1 , λ[t]]dt , λ[t] := t λ2 + (1 − t)λ1
σ s [ε1 , λ1, λ2] :=
∂λ
0
0
1
1
∂σ[ε[t], λ2]
dt =
E s [λ2 , ε1 , ε2 ] :=
E t [ε[t], λ2]dt , ε[t] := t ε2 + (1 − t)ε1
∂ε
0
0
La seconda parte si ottiene osservando che E s pu`
o essere ricavata come media, sul segmento
t ∈ [0 · · ·1], della matrice E t ed utilizzando i lemmi 2.5.1 e 2.5.3. La terza parte, inﬁne, si
ottiene osservando che, per le (1.1.5) e (2.1.4), risulta
(σ a [ε1 , λ1] − ασ a[ε1 , λ2])Tεp11 > 0
(σ a [ε1 , λ2] − ασ a[ε1 , λ1])Tεp12 > 0
se λ1 − λ2 < 0 , εp11 = 0
se λ2 − λ1 < 0 , εp12 = 0
a,
e che, per λ2 → λ1 si ottiene, per continuit`
σ s [ε1 , λ1, λ1 ]εp11 < 0
se εp11 = 0
Capitolo 3
Un metodo iterativo di soluzione
Nel presente capitolo `e descritto un metodo di soluzione del problema dello shakedown
illustrato nel capitolo 2, appropriato per una analisi agli elementi ﬁniti. Il metodo `e basato
su un processo iterativo incrementale che produce una sequenza λ(k) , k = 1, 2, · · · di
moltiplicatori che, a meno dell’errore tipico legato alla discretizzazione, sono strettamente
sicuri secondo la (2.2.1) e converge monotonicamente al fattore di sicurezza a shakedown λa.
3.1
Generica formulazione in elementi finiti
Assumendo una discretizzazione in elementi ﬁniti della struttura, sia u ∈ N il vettore
degli spostamenti nodali liberi, u0 il suo valore iniziale e ε[u] := D(u − u0 ) l’incremento di
deformazione, cinematicamente compatibile, associato all’incremento di spostamento u−u0
¯ possiamo deﬁnire il vettore della risposta
dalla matrice di compatibilit`
a D[x]. Per λ ≤ λ
nodale s[u, λ] ∈ N associato ad u:
DT σ[u, λ]dv , σ[u, λ] := σa [σ 0 + Eε, λ]
(3.1.1)
s[u, λ] :=
B
e la matrice elastica K e ∈ N × N simmetrica e deﬁnita positiva
DT EDdv , K e = K Te > 0
K e :=
(3.1.2)
B
mediante le equivalenze energetiche:
δεT σ[u, λ]dv , δuTK e δu ≡
δεTEδεdv
δuTs[u, λ] ≡
B
,
∀δu , δε := Dδu
B
(3.1.3)
in cui E `e la matrice elastica che collega il vettore della deformazione ε a quello della
tensione σ e σ 0 `e lo stato di tensione iniziale. La tensione σ[u, λ] ∈ Es [λ] per deﬁnizione.
Inoltre, essendo σ[u, λ] uno stato di tensione auto–equilibrato si ha:
s[u, λ] = 0
(3.1.4)
Con queste notazioni, il problema dello shakedown pu`
o essere riformulato in forma discreta.
La determinazione del fattore di sicurezza a shakedown pu`o essere ricondotta alla soluzione
del seguente problema:
(3.1.5)
λa := max λ : ∃u : s[u, λ] = 0
22
Un metodo iterativo di soluzione
23
Si osservi che, cos`ı formulato, il problema dello shakedown si presenta simile alla formulazione discreta del teorema statico dell’analisi limite. L’unica diﬀerenza `e il ruolo svolto
dal fattore di sicurezza λ che agisce come parametro interno della funzione fs [σ, λ] invece
di essere un moltiplicatore esterno dei carichi. Tuttavia, questa diﬀerenza non pu`
o essere
considerata signiﬁcativa: infatti, in caso di carico monotono, possiamo formulare l’analisi
limite nella stessa forma (3.1.5) se si usa il fattore di sicurezza come ampliﬁcatore della
tensione elastica invece che dei carichi esterni.
La stretta analogia tra i due problemi suggerisce che i metodi disponibili per risolvere
i problemi di analisi limite possano essere direttamente estesi alla soluzione dei problemi
di shakedown. Il metodo proposto in [10] e di seguito riportato, in eﬀetti, segue questa
linea. Esso pu`
o essere considerato come una specializzazione allo shakedown dell’algoritmo
incrementale–iterativo strain driven utilizzato in analisi elastoplastica [22] e corrisponde ad
una estensione diretta del metodo arc–length path–following descritto in [6].
Prima di entrare nei dettagli del metodo di soluzione proposto in [10], `e utile fare alcuni
commenti. Per prima cosa, essendo le leggi costitutive deﬁnite localmente, l’algoritmo di
rientro
(3.1.6)
σ := σ a [σ0 + Eε, λ] , ε := D(u − u0 )
pu`
o essere applicato separatamente a ciascun elemento oppure, se l’elemento `e deﬁnito
mediante integrazione numerica, a ogni punto di Gauss. Pertanto, la valutazione del vettore
s[u, λ] per assegnati u e λ attraverso la (3.1.1) risulta in eﬀetti molto semplice e rapida.
Alcune cautele tuttavia devono essere usate per ridurre l’errore di discretizzazione ed evitare
incoerenze. E’ ben noto che discretizzazioni strettamente compatibili possono produrre
locking (ved. [23], [24]) e che un uso poco accorto di procedure di collocazione possono far
perdere la connotazione energetica nel modello discreto.
Rimandando alla letteratura per un approfondimento dell’argomento, in generale si
suggerisce che il campo di tensione σ sia deﬁnito estendendo la condizione di Haar–Karman
(2.4.3) all’intero elemento, attraverso la condizione
1
(3.1.7)
(σ − σ ∗ )TE −1 (σ − σ ∗ ) dv = min , σ[x] ∈ Es [λ] ∀x ∈ X
φe :=
2
Be
dove Be `e il volume dell’elemento, σ ∗ := σ 0 + Eε. L’insieme X := {x1 , · · · xn } contiene i
punti di controllo che consentono di discretizzare, in termini approssimati, la condizione di
ammissibilit`
a plastica σ ∈ Es [λ] , ∀x ∈ Be , e σ `e auto equilibrato sull’elemento o, quanto
meno, il suo lavoro compiuto in Be su ogni δε := Dδu sia nullo. L’implementazione della
condizione (3.1.7) risulta in generale relativamente semplice.
Si noti, inoltre, che per il lemma 2.5.1 e 2.5.4, se si considerano due stati diﬀerenti
o scrivere
{u1 , λ1 } e {u2 , λ2 }, si pu`
σ[u2 , λ2 ] − σ[u1 , λ1] = E s [λ2 , ε1 , ε2 ](ε2 − ε1 ) + σ s [ε1 , λ1 , λ2](λ2 − λ1 )
dove ε1 = Du1 , ε2 = Du2 . Integrando sul volume B ed introducendo gli operatori secanti
⎧
⎪
⎪
DT E s [λ2 , ε1 , ε2 ]D dv
⎨ K s [λ2 , u1 , u2 ] :=
B
(3.1.8)
⎪
⎪
⎩ y s [u1 , λ1, λ2] :=
DT σs [ε1 , λ1 , λ2] dv
B
Un metodo iterativo di soluzione
24
si ottiene
s[u2 , λ2] − s[u1 , λ1] = K s [λ2 , u1 , u2 ](u2 − u1 ) + (λ2 − λ1 )ys [u1 , λ1, λ2 ]
(3.1.9)
Per il Lemma 2.5.4, la matrice K s , simmetrica per costruzione, soddisfa le condizioni
0 ≤ δuT K s δu ≤ δuT K e δu ,
∀δu
(3.1.10)
¯ , ∀j. Il moltiplicatore λ
¯ pu`
o essere
Inﬁne, l’algoritmo di rientro (3.1.6) richiede λj ≤ λ
ottenuto semplicemente come il pi`
u piccolo valore di λ che rende vuoto l’insieme Es su
almeno uno dei punti di controllo xi ∈ B. Ci`
o si ottiene facilmente, una volta per tutte,
all’inizio dell’analisi.
3.2
Soluzione iterativa del problema dello shakedown
Determineremo lo stato limite di shakedown, cio`e il moltiplicatore λa , le autotensioni associate σ a ed i corrispondenti spostamenti ua , attraverso una sequenza di stati sicuri x(k) :=
{λ(k) , σ (k) , u(k) } , k = 1, 2, · · ·, a partire dallo stato al limite elastico x(0) := {λe , 0 , 0}.
La sequenza comporta s(k) := s[u(k) , λ(k)] = 0 , ∀k e risulta monotona non decrescente in
λ(k) e viene arrestata non appena `e raggiunto lo stato limite di shakedown che si ottiene
quando λ(k+1) = λ(k).
In ciascun passo del processo, il nuovo stato x(k) si ricava dal precedente x(k−1) mediante
un algoritmo iterativo che, a partire da x0 := x(k−1) , produce una sequenza convergente
xj := {λj , σ j , uj } , j = 1, 2, · · · aggiornando ricorsivamente il vettore spostamento ed il
moltiplicatore dei carichi:
uj+1 := uj + u˙ j
(3.2.1)
λj+1 := λj + λ˙ j
dove le correzioni u˙ j e λ˙ j sono deﬁnite in modo da soddisfare, almeno approssimativamente, la condizione di equilibrio s[uj+1 , λj+1 ] = 0 richiesta dalla (3.1.5) e σ j+1 `e deﬁnito
dall’algoritmo di rientro (3.1.6).
Ottenere σ j+1 da λj+1 e uj+1 `e, per quanto detto, relativamente semplice, ci concentreremo quindi sull’aggiornamento di uj e λj , ottenuto utilizzando le correzioni u˙ j e λ˙ j
deﬁnite come soluzione del sistema lineare
u˙ j
K e yj
−s[uj , λj ]
(3.2.2)
=
y Tj 0
0
λ˙ j
in cui K e `e la matrice di rigidezza elastica deﬁnita dalla (3.1.3), il vettore sj `e il residuo
all’equilibrio
(3.2.3)
sj := s[uj+1 , λj+1 ]
ed il vettore yj `e deﬁnito, dalla (3.1.8), come
1
DT σs [εj , λj , λj+1 ] dv =
(s[uj , λj+1 ] − s[uj , λj ])
yj :=
λj+1 − λj
B
(3.2.4)
˙ j } sono ottenuti facilmente dal sistema (3.2.2):
Gli incrementi {λ˙ j , u
λ˙ j = −
y Tj K −1
e sj
y Tj K −1
e yj
,
−1
˙
u˙ j = −K −1
e sj − λj K e y j
(3.2.5)
Un metodo iterativo di soluzione
25
Tabella 3.2. Schema dell’algoritmo proposto
1. Assemblaggio e decomposizione della matrice K E
2. Calcolo delle soluzione elastiche {uei , σei } per ciascu carico pi
ui := K −1
E pi
,
σ i := EDui
3. Calcolo del moltiplicatore al limite elastico λe ed inizializzazione del processo incrementale
u(0) := 0 , σ (0) = 0 , λ(0) := 0 , u(1) := 0 , σ (1) = 0 , λ(1) := λe
4. Ripetere, per k = 2, 3, ...
(a) Inizializzazione del passo mediante la (3.2.8):
λ1 := λ(k−1) + β(λ(k−1) − λ(k−2) )
u1 := u(k−1) + β(u(k−1) − u(k−2))
(b) Ripetere per j=1, 2, ...
- Calcolo σ j mediante la (3.1.6)
- Calcolo dei vettori sj and y j mediante la (3.2.2) e (3.2.7)
- Correzione iterativa del vettore xj+1 :
λj+1 = λj + λ˙ j
y T K −1 s
˙j = − j E j
,
λ
−1
yT
uj+1 = uj − K −1 (sj + λ˙ j y )
j K E yj
j
E
Fino a che ||sj || ≤ tol1
(c) Aggiornamento della soluzione corrente
λ(k) = λj
,
u(k) := uj
,
σ (k) := σ j
Fino a che λ(k) − λ(k−1) ≤ tol2
5. Al termine si ottiene la soluzione del problema dello shakedown
λa := λ(k)
ua := u(k)
,
,
σ a := σ (k)
Il valore λj+1 = λj + λ˙ j `e richiesto nella (3.2.4) per ottenere y j , pertanto la (3.2.4) `e
accoppiata alla prima delle (3.2.5). Tuttavia λj+1 potrebbe essere ottenuto per iterazione
come limite della sequenza
˜ Ti K −1
e sj
˜ i := λj − y
λ
T
˜ i K −1
˜i
y
e y
,
˜ i :=
y
1
˜ i] − sj )
(s[uj , λ
˜ i − λj
λ
(3.2.6)
che pu`
o essere inizializzata assumendo come primo valore di y j la tangente iniziale
˜ 1 := y[uj , λj , λj ]
yj ≈ y
(3.2.7)
La sequenza implementa il metodo classico di iterazione secante che converge sempre e anche
in modo rapido. In eﬀetti, considerando anche che la sua soluzione `e utilizzata all’interno
di uno schema iterativo esterno, `e suﬃciente fermarsi alla prima iterazione.
Un metodo iterativo di soluzione
26
Lo schema iterativo dell’algoritmo `e descritto nella tabella 3.2. Il processo iterativo
(3.2.1) e (3.2.2), al passo k, `e inizializzato assumendo:
λ1 := λ(k−1) + β (k) (λ(k−1) − λ(k−2))
u1 := u(k−1) + β (k)(u(k−1) − u(k−2) )
(3.2.8)
o essere collegato al numero delle iterazioni
dove β (k) `e un fattore di ampliﬁcazione che pu`
(k−1)
eﬀettuate nel passo precedente
n
β (k) =
n(k−1) − n
¯
n(k−1) + n
¯
en
¯ il suo valore medio desiderato (tipicamente n
¯ ≈ 4 · · · 8). La stessa formula `e usata nel
primo passo del processo incrementale a partire dai vettori:
{u(0) = 0 , σ (0) = 0 , λ(0) = 0} ,
{u(1) = 0 , σ (1) = 0 , λ(0) = λe }
e assumendo β1 abbastanza piccolo (β1 ≈ 0.01).
L’iterazione `e arrestata e si considera concluso il passo quando
sj ≤ tol1
,
sj 2 := sTj K −1
e sj
(3.2.9)
Inﬁne, il processo incrementale `e arrestato quando λ(k)) = λ(k−1), all’interno di una tollera che ci`o corrisponde al raggiungimento del limite di shakedown.
anza assegnata tol2 . Si vedr`
Dimostreremo che:
1. La sequenza xj prodotta dallo schema iterativo (3.2.1), (3.2.2) converge ad una nuova
soluzione x(k) := {λ(k) , σ (k) , u(k) } che soddisfa le σ (k) ∈ Es [λ(k)] e s[u(k) , λ(k)] = 0
e quindi fornisce λ(k) ≤ λa.
2. La successione x(k) soddisfa la condizione λ(k) > λ(k−1) se λ(k−1) < λa . Inoltre, il
caso λ(k) = λ(k−1) implica λ(k) = λa .
3.3
Convergenza dello schema iterativo
Per sempliﬁcare la discussione sulla convergenza dello schema iterativo proposto, `e conveniente riferirsi alla soluzione {µi , z i } del problema ad autovalori
[µK e − K j ]z = 0
Essendo entrambe le matrici K j e K e simmetriche e la prima deﬁnita positiva, la soluzione
a caratterizzata dalle ben note condizioni:
µi and z i sar`
1
if
i
=
j
µi if i = j
, z Ti K j z j =
(3.3.1)
z Ti K e z j =
0 if i = j
0 if i = j
Inoltre, per la (3.1.10), si ha
0 ≤ µi ≤ 1 ∀i
(3.3.2)
Un metodo iterativo di soluzione
27
Teorema 3.3.1 (convergenza dello schema iterativo). La sequenza xj prodotta dallo schema (3.2.1) e (3.2.2) `e convergente ed il suo limite x(k) := {λ(k) , σ (k) , u(k) } `e
caratterizzato dalle condizioni σ (k) ∈ Es [λ(k)] e s[uk , λk ] = 0.
Dimostrazione. Per deﬁnizione (3.2.5) risulta
sj+1 − sj = K j u˙ j + λ˙ j y j
Si ha quindi
λ˙ j := −
,
y Tj K −1
e sj
˙
u˙ j := −K −1
e (sj + λj y j )
,
y j K −1
e yj
T
˙
sj+1 = sj + λ˙ j yj − K j K −1
e (sj + λj y j )
Pertanto, espandendo yj , sj e sj+1 in termini di z i
yj = K e
βi z i , sj = K e
αi z i
i
,
sj+1 = K e
i
α
˜izi
i
si ottiene, per la (3.3.1),
α
˜i = (1 − µi )(αi + λ˙ j βi )
e quindi
sj+1 2 =
α
˜2i =
i
(1 − µi )2 α2i − λ˙ 2j βi2 + 2λ˙ j βi (αi + λ˙ j βi )
i
Per la (3.1.1), si ha
˙
yTj u˙ j = −y Tj K −1
e (sj + λy j ) = 0
=⇒
βi(αi + λ˙ j βi) = 0
i
Risulta quindi
sj+1 2 =
(1 − µi )2 (α2i − λ˙ 2j βi2 ) ,
sj 2 =
i
α2i
i
e per la (3.3.2) si ha sempre (1 − µi )2 ≤ 1, e quindi rsulta anche:
sj+1 <1
sj se
∃ i : µi αi = 0
oppure
λ˙ j = 0 .
Pertanto la sequenza {||sj ||} `e monotonamente non crescente. Essa `e anche limitata inferiormente (sj ≥ 0) per cui `e convergente ad un valore ﬁnito positivo. Tuttavia si dimostra
che non pu`
o essere limj sj > 0. Infatti, l’assunzione:
lim sj > 0
j→∞
implicherebbe
lim λ˙ j = 0 ,
j→∞
lim µi αi = 0 ∀i
j→∞
⇒
lim (sj+1 − sj ) = 0
j→∞
Ci`
o signiﬁca che la sequenza {sj } converge ad una valore ﬁnito positivo ¯s
¯ = lim sj = 0
s
j→∞
Un metodo iterativo di soluzione
28
Pertanto, per la (3.2.5) e (3.2.1) risulta:
¯ := K −1
¯ = 0
lim (uj+1 − uj ) = u
e s
j→∞
ed inoltre, si veriﬁca che:
(σj+1 − σ j )T (εj+1 − εj ) dv ≡ lim (sj+1 − sj )T (uj+1 − uj ) = 0
lim
j→∞ B
j→∞
Ci`
o, per il lemma 2.5.2, implica limj→∞ (σ j+1 − σ j ) = 0, ed inoltre:
p
p
¯ = 0
lim (εj+1 − εj ) := lim (εj+1 − εj ) = ¯ε := Du
j→∞
j→∞
ci`o comporta
ε
lim εpj+k = εpj + k¯
j→∞
1 p
¯ = Du
¯ = 0
ε
=ε
j,k→∞ k j+k
=⇒
lim
Da questa si ottiene
¯ y j+k
lim (uj+k+1 − uj+k ) yj+k = u
T
T
j,k→∞
1
=
k
B
εpj+k σ s [εj+k , λj+k , λj+k ] dv < 0
che comporta un assurdo essendo in contrasto con la (3.2.1, 3.2.2) che fornisce (uj+k+1 −
uj+k )T yj+k = 0 , ∀j, k. Risulta pertanto
lim sj = 0
j→∞
Da questa si ottiene
lim (uj+1 − uj ) = lim K −1
e sj = 0
j→∞
j→∞
lim (λj+1 − λj ) = lim
,
j→∞
j→∞
yTj K −1
e sj
yTj K −1
e yj
=0
e quindi anche le sequenze {uj } e {λj } sono convergenti e possiamo deﬁnire
u(k) = lim uj
j→∞
,
λ(k) = lim λj
j→∞
Inoltre, essendo σ a [u, λ] continua in u e λ, σ j := σa [uj , λj ] `e anche convergente. Pertanto,
risultando σj ∈ Es [λj ] , ∀j, per costruzione, si ottiene
σ (k) = lim σ j ∈ Es [λ(k)]
j→∞
Come si voleva dimostrare.
Un metodo iterativo di soluzione
3.4
29
Convergenza alla soluzione di adattamento
Il teorema che segue dimostrer`a che la sequenza x(k) prodotta dal processo incrementale
converge alla soluzione del problema dello shakedown.
Teorema 3.4.1 (convergenza del processo incrementale al limite di shakedown).
La sequenza λ(k) `e caratterizzata da
λ(k−1) ≤ λ(k) ≤ λa
Inoltre, il caso λ(k) = λ(k−1) con u(k) = u(k−1) implica il raggiungimento dello stato limite
di shakedown:
λk = λa
¯ e σ (k) ∈ Es per ogni k, per costruzione. Pertanto la
Dimostrazione. Risulta σ (k) ∈ S
(k)
condizione λ ≤ λa deriva direttamente dalla (2.2.1). Poich`e le σ (k) sono delle autotensioni
e le ε(k) delle deformazioni cinematicamente compatibili, per ogni valore di k, anche una
qualsiasi loro combinazione `e un’autotensione o un campo di deformazione cinematicamente
compatibili. Pertanto, tenendo conto della (3.1.4), si ha:
(σ(k) − βσ (k−1) )T ε(k) dv = (s(k) − βs(k−1) )T(u(k) − u(k−1) ) = 0
B
Posto
ε(k) := ε(ke) + ε(kp)
e ricordando la (3.1.1), si ottiene:
T
(k)
(k−1) T (kp)
(σ − βσ
) ε dv = −
ε(ke) Eε(ke) dv ≤ 0
B
B
Questa implica λ(k) ≥ λ(k−1) . Assumiamo, infatti, λ(k) < λ(k−1). Poich`e il dominio elastico
E `e convesso e contiene la tensione tensione nulla, `e valida la relazione:
λ(k−1)σ e + σ (k−1) ∈ E ⇒ β λ(k−1)σ e + σ (k−1) ∈ E con 0 ≤ β < 1
Posto β tale che risulti λ(k) := βλ(k−1), si ottiene:
λ(k)σ e + βσ (k−1) ∈ E
Per la deﬁnizione di dominio ammissibile a shakedown (2.1.2) si ha:
σ (k−1) ∈ Es [λ(k−1)]
,
βσ (k−1) ∈ Es [λ(k)]
da cui deriva la condizione:
βEs [λ(k−1)] ⊆ Es [λ(k)]
se
λ(k−1) > λ(k)
Utilizzando la condizione di Druker (2.1.6) si ricava:
(σ(k) − βσ (k−1) )Tε(kp) dv ≥ 0
B
Un metodo iterativo di soluzione
30
in contrasto con l’assunto che tale quantit`
a debba essere minore di zero e quindi λ(k) ≥ λ(k−1).
Il caso λ(k) = λ(k−1) `e possibile solo se si veriﬁca la condizione ε(ke) = 0, per cui si ha:
ε(kp) = ε(k) = D(u(k) − u(k−1) )
¯ si ottiene
cio`e ε(kp) ∈ K e ε(kp) = 0. Ricordando che σ (k) ∈ S,
T
σ(k) ε(kp) dv = 0
B
dove σ(k) e ε(kp) sono associati dalla legge di ﬂusso (2.4.2), per costruzione. Tale condizione
implica λ(k) ≥ λa, per la (2.2.2). Inoltre, dovendo essere λ(k) ≤ λa, si ottiene λ(k) = λa,
come si voleva dimostrare.
3.5
Alcuni commenti
Il processo di soluzione descritto pu`
o essere visto come un processo di analisi incrementale
step–by–step che simula il caso di un’evoluzione di carico in aumento tale che, in ogni passo
k, il carico `e inteso come tutti i possibili valori di λ(k) P per cui sia possibile conseguire
adattamento elastico.
Cos`ı interpretato, il processo di soluzione proposto corrisponde al processo di soluzione
standard path–following descritto in [6]. Lo schema di iterazione (3.2.1)–(3.2.2) corrisponde
ad una implementazione del metodo arc–length proposto da Riks ([25], [26]). Tale schema,
pu`
o essere riscritto nella forma [27]:
˜
u˙ j
−sj
˜
K
y
(3.5.1)
=
0
λ˙ j
∆uTM γ
˜ `e una approssimazione dell’Hessiano K t := ∂s/∂u e il vettore y
˜ `e una
in cui, la matrice K
a ∆u := uj − u0 , M e γ esprimono
approssimazione del vettore y t := ∂s/∂λ. Le quantit`
la metrica utilizzata per misurare la distanza nello spazio {u, λ}, nella forma:
T ∆u
M ·
∆u
2
(3.5.2)
{∆u, ∆λ} :=
∆λ
· γ
∆λ
La (3.2.2) si ottiene dalla (3.5.2) assumendo:
yt := yj
,
˜ := K e
K
,
M ∆u := y j
,
γ=0
(3.5.3)
Ovviamente possono essere fatte anche scelte diverse. Le particolari assunzioni eﬀettuate
sono convenienti poich`e rendono agevole la dimostrazione della convergenza dello schema
di iterazione e del processo di soluzione. Una simile assunzione `e stata proposta in [6] nella
implementazione del processo di soluzione path–following adoperato in problemi di analisi
limite.
` utile notare che, oltre a λ(k) e σ (k) , il processo di soluzione produce gli spostamenti
E
(k)
u e la deformazione plastica totale (cumulata):
k
i=0
ε(ip) :=
k Du(i) − E −1 σ (i) .
i=0
(3.5.4)
Un metodo iterativo di soluzione
31
Queste quantit`
a sono strettamente correlati all’evoluzione del carico. Oltre ad oﬀrire informazioni utili sul processo di adattamento elastico, possono essere presi come valori
di riferimento per essere usati, all’interno di un criterio di veriﬁca di sicurezza convenzionale, per valutare fattori di sicurezza per stati limite di spostamento o stati limite di
deformazione plastica (uno stato di sicurezza imposto sulla limitazione della deformazione
plastica, [28]–[29], sembra essere appropriato per scopi tecnici).
Inﬁne, possiamo osservare che il processo di soluzione proposto in [10] pu`
o essere usato
direttamente in problemi di analisi di limite. Infatti, nel caso in cui il dominio dei carichi `e
rappresentato da un solo punto p, l’analisi di adattamento si riconduce semplicemente ad
una analisi limite in cui il carico cresce in modo proporzionale al fattore λ.
Le diﬀerenze con il processo standard di soluzione tipo path–following `e legata essenzialmente all’uso dello tensione aggiuntiva σ ∈ Es [λ]. Il metodo path–following durante
il processo di soluzione utilizza la tensione tatale σ ∈ E ed usa il moltiplicatore λ per
ampliﬁcare direttamente il carico ([26], [6]). In questo caso l’equilibrio `e espresso dalla
condizione:
r[u, λ] := s[u] − λp = 0
(3.5.5)
in cui il vettore s `e la risposta nodale corrispondente alla tensione σ associata allo spostamento u ed espressa come:
σ := σa [σ 0 + Eε, 0] ,
σ∈E
Di conseguenza, lo schema correttivo (3.2.2) diventa [6]:
u˙ j
−r[uj , λj ]
K e −p
,
=
−pT 0
0
λ˙ j
(3.5.6)
(3.5.7)
in cui, −p := ∂r/∂λ `e usato al posto di yj . Ovviamente, per ogni campo di deformazione
cinematicamente compatibile εp := Dup , si ha:
σ T εp dv = pT up > 0
(3.5.8)
B
che garantisce la convergenza per le stesse motivazioni riportate nei teoremi 3.3.1 e 3.4.1.
Per l’equivalenza che si ha tra i due schemi, ci si aspettano anche stesse caratteristiche di
convergenza. Ci`
o `e dimostrato dai risultati numerici. Si noti comunque, che lo schema
o non `e
(3.2.2) richiede la valutazione del vettore yj ad ogni ciclo di iterazione, mentre ci`
necessario nel caso dello schema standard (3.5.7). Per questo motivo lo schema proposto
non `e particolarmente conveniente in problemi di analisi limite in cui i processi di carico
che si analizzano sono relativi ad una singola condizione.
Capitolo 4
L’algoritmo di rientro sul dominio
ammissibile
Ingredienti principali per l’analisi a shekedown sono il dominio ridotto Es [λ] e l’algoritmo
di rientro nel dominio. La costruzione del dominio di adattamento deﬁnito dalla (2.1.2)
richiede che sia nota la funzione fs [σ, λ]. In generale, l’espressione analitica di tale funzione
non `e nota a priori ma `e possibile determinare solo per punti la funzione f [σ] da cui
essa deriva. Si rende necessario, quindi, formulare la costruzione del dominio ridotto di
shakedown per contesti del tutto generali.
4.1
Il dominio ridotto per una generica funzione di snervamento
In generale, la funzione f [σ] non `e nota in forma analitica ma si pu`
o determinare il suo
andamento mediante una costruzione per punti, ciascuno dei quali `e deﬁnito dalla coppia
{σ k , nk } tale che:
∂f ,
nk :=
(4.1.1)
f [σk ] = 0
∂σ σ k
La coppia {σ k , nk } deﬁnisce un piano tangente alla superﬁcie del dominio nel punto σ k .
Usando m coppie, `e possibile deﬁnire un poliedro di m-facce esterno alla funzione f [σ]
che approssima la superﬁcie del dominio elastico con un insieme di superﬁci piane f [σ] :=
{f1 [σ], f2[σ], · · · , fm [σ]}T che deﬁniscono il dominio nella forma:
fk [σ] ≡ {nTk σ − ck } ≤ 0 ,
ck = σ Tk nk
k = 1, · · · , m
(4.1.2)
La funzione di shakedown f s [σ, λ] := {fs1 [σ, λ], · · · , fsm [σ, λ]}T, introdotta con la (2.1.1),
si ottiene dalla (4.1.2) mediante la semplice espressione:
f s [σ, λ] := N T σ − c + λ¯b
in cui la matrice N raccoglie nelle colonne le normali ai piani:
∂fsm [σ, λ]
∂fs1 [σ, λ]
N = N [σ, λ] := n1 , · · · , nm =
, ··· ,
∂σ
∂σ
32
(4.1.3)
(4.1.4)
L’algoritmo di rientro sul dominio ammissibile
mentre i vettori c e b sono deﬁniti come:
c := {c1 , · · · , cm}T , b¯ := {¯b1 , · · · , ¯bm}T
33
con ¯bk := max {nTk σ e [t]}
σ e ∈ Se
(4.1.5)
Ricordando la (1.2.2), che deﬁnisce il dominio delle sollecitazioni elastiche Se , si ha:
min T
p
αi nk σ ei se nTk σ ei < 0
¯bk =
(4.1.6)
aki
con
aki :=
nTk σ ei se nTk σ ei ≥ 0
αmax
i
i=1
¯ `e valutato una sola volta all’inizio dell’analisi. La rappresentazione lineare a
Il vettore b
tratti del dominio Es [λ] consente di determinare, in ogni punto di Gauss p, il valore minimo
λp oltre il quale il dominio di shakedown Es [λ] risulta vuoto. Infatti, considerati due piani
di normali opposte nk e −nk , e indicando i piani corrispondenti con l’apice ±, si ha:
+
ck + c−
k
λp := min +
(4.1.7)
k
bk + b−
k
¯ deﬁnito dalla (2.1.3) si determina facilmente come il minimo
Il valore del moltiplicatore λ
dei valori di λp fatto su tutti i punti di Gauss:
¯ = min{λp}
(4.1.8)
λ
p
tale valore coincide con quello del moltiplicatore denominato, in letteratura, di adattamento
plastico.
La rappresentazione lineare a tratti del dominio elastico si presenta come una formulazione che pu`
o essere applicata in contesti del tutto generici. Come sar`
a dimostrato in
seguito dai risultati numerici ottenuti, il numero di piani che approssimano il dominio elastico inﬂuenza la soluzione ﬁnale. L’inﬁttimento dei piani aumenta l’accuratezza ma peggiora
l’eﬃcienza computazionale, in quanto l’algoritmo di rientro del dominio risulta poco eﬃciente quando il numero di piani diventa elevato. Si rende necessario, quindi, eﬀettuare una
scelta accurata del numero di piani e della loro distribuzione.
4.2
L’algoritmo di return–mapping in plasticit`
a multifalda
Partendo da un punto noto {λ0 , u0 , σ0 } e da un predittore elastico σ ∗ = σ 0 + Eε[u] non
necessariamente contenuto in Es [λ], la tensione ammissibile σ a [ε, λ] corrispondente a u e
¯ si ottiene dal seguente algoritmo di rientro nel dominio ammissibile:
λ ≤ λ,
σ := σ a[ε, λ] = σ ∗ − Eεp = σ ∗ + ∆σ
,
σ ∈ Es [λ]
(4.2.1)
che ci consente di determinare la tensione σ a contenuta in Es [λ], corrispondente alla dea multifalda,
formazione εp associata attraverso la legge di ﬂusso (2.1.8), che per plasticit`
diventa:
µk fk = µk f˙k = 0
p
p
T
∀k
(4.2.2)
,
ε := ε [µ, N ] = N µ , µ := {µ1 , · · · , µm }
µk ≥ 0
in cui N `e la matrice deﬁnita dalla (4.1.4) mentre il vettore µ raccoglie i moltiplicatori
plastici associati a ciascuno dei piani che individuano la superﬁcie del dominio di shakedown.
Ricordando la (4.1.3), la condizione di ammissibilit`
a della tensione σa si scrive semplicemente:
(4.2.3)
f s [σ, λ] := N T (σ∗ − Eεp ) − c + λ¯b ≤ 0
L’algoritmo di rientro sul dominio ammissibile
4.3
34
Il problema di rientro come problema di programmazione
quadratica QP
La condizione di Haar–Karman (2.4.3) che, in un contesto di plasticit`
a multifalda diventa:
⎧
⎨ min∆σ : ΠHK [∆σ] ≡ 1 ∆σT F ∆σ
2
con b = c − N Tσ ∗ − λ¯b (4.3.1)
QP :
⎩
subject to : f s [∆σ, λ] ≡ N T ∆σ − b ≤ 0
consente di ottenere l’algoritmo di rientro nel dominio ammissibile deﬁnito dalle (4.2.1)–
(4.2.3) nella forma di un problema di programmazione quadratica strettamente convesso.
Indichiamo con K = {1, 2, .., m} il set completo delle condizioni di vincolo e con QP (K)
il problema (4.3.1). Risolto il problema QP (K) deﬁniamo set attivo A l’insieme dei vincoli
per cui si veriﬁcno le condizioni:
fsk [∆σ, λ] = 0
,
µk > 0
∀k ∈ A
(4.3.2)
cio`e, in soluzione solo i vincoli attivi sono veriﬁcati con l’uguaglianza. Indicando con N a
la matrice che raccoglie le ma normali nk dei piani attivi, e con un pedice a le quantit`
a
relative ai soli piani attivi, in soluzione si ha la condizione ottimale:
N a µ[λ] = −E −1 ∆σ
(4.3.3)
I moltiplicatori plastici si determinano facilmente considerando che in soluzione sono valide
entrambi le condizioni:
σ = σ ∗ − EN aµ
(4.3.4)
N Ta σ − ca + λ¯ba = 0
da cui si ottiene:
µ[λ] = H −1 N Ta σ ∗ − ca + λ¯ba
con
H := N Ta EN a
(4.3.5)
La matrice tangente K t e i vettori yj e sj usati nel procedimento iterativo di soluzione si
determinano ricorrendo alle relazioni incrementali:
⎧
p
⎨ dσ[ε, λ] = E(dε − dε )
T
(4.3.6)
⎩ df [σ, λ] = ∂f dσ + ∂f dλ
∂σ
∂λ
Assumendo che le normali ai piani nk sono costanti in λ e σ, l’incremento di deformazione
plastica diventa semplicemente:
dεp := dεp [µ, N [σ, λ]] = N dµ
(4.3.7)
e di conseguenza, per l’incremento della funzione di shakedown si ottiene:
∂f
∂f T
dλ
E(dε − dεp ) +
∂σ
∂λ
= N Ta E(dε − N a dµ) + f ,λ dλ
df [σ, λ] =
(4.3.8)
L’algoritmo di rientro sul dominio ammissibile
35
Dovendo essere, nel caso di deformazione plastica, dµ > 0 e df = 0, si ricava facilmente:
(4.3.9)
dµ = H −1 N Ta Edε + f ,λ dλ
che sostituito nella (4.3.6) fornisce:
¯
dσ[ε, λ] = (E − EN aH −1 N Ta E)dε − EN a H −1 bdλ
(4.3.10)
Dalla relazione incrementale:
dσ[ε, λ] =
∂σ[ε, λ]T
∂σ[ε, λ]
dε +
dλ = E tdε + σ t dλ
∂ε
∂λ
(4.3.11)
si riconoscono le quantit`
a:
E t := E − EN a H −1 N Ta E
,
σ t := −EN a H −1 b¯
(4.3.12)
in cui E t `e la matrice tangente locale. La matrice tangente del problema algebrico, usata
in sostituzione della matrice secante K s non ancora nota a inizio passo, ed i vettori yj e sj
usati nel procedimento iterativo di soluzione possono essere valutati con una procedura di
assemblaggio delle quantit`a locali:
K t = A(DT E t[uj , λj ]D) ,
4.4
yj = A (DT σt [uj , λj ])
,
sj = A (DTσ[uj , λj ]) .
(4.3.13)
Cenni ai metodi di soluzione per problemi QP e superiorit`
a dei metodi di proiezione nel caso in esame
Gli algoritmi per risolvere il problema QP possono essere suddivisi in due categorie [30]:
varianti del metodo del simplesso e metodi di proiezione o Active set methods. I primi, come
l’algoritmo di Lemke, compiono operazioni di pivot su matrici di dimensione proporzionali
a (n × m) dove n `e la dimensione di σ e m `e il numero complessivo di vincoli del problema.
I secondi, invece, proiettano il problema in un sottoinsieme deﬁnito dai soli vincoli attivi
contenuti nel set attivo Ah ⊂ K e utilizzano operatori di dimensione al massimo (n × n).
In genere la dimensione dello spazio delle tensioni `e di poche unit`
a mentre il numero di
vincoli deve essere suﬃcientemente alto per consentire una accurata descrizione del dominio
elastico. Ci`o rende pi`
u eﬃciente l’uso di algoritmi del tipo Active Set methods per la
soluzione dei problemi di programmazione quadratica.
4.5
Il metodo di Golfbard e Idnani
Il problema di programmazione quadratica QP, cos`ı come riportato dalla (4.3.1), `e espresso
in forma primale, cio`e assumendo, come incognite del problema, le ∆σ dette appunto
variabili primali nel linguaggio dell’ottimizzazione. Analogamente, utilizzando la (4.3.3), lo
stesso problema pu`
o essere trasformato nella forma duale in cui i parametri incogniti sono
appunto le variabili duali µ.
Il metodo di Goldfarb e Idnani, detto Dual Active Set Method, sfrutta entrambe le propiet`a del problema considerato nella forma primale e nella forma duale. Il metodo `e di tipo
L’algoritmo di rientro sul dominio ammissibile
36
iterativo e genera una sequenza di vettori µh e ∆σh che soddisfano le condizioni di Kunh–
Tucker del problema originale, eccetto l’ammissibilit`
a delle variabili primali (f [∆σh , λ] ≤ 0)
per lo stesso problema. La soluzione ottenuta al ciclo h `e un punto ottimale, in termini di
variabili primali ∆σ, per il subproblema QP(A) deﬁnito a partire dal problema QP(K) considerando solo un sottoinsieme di vincoli linearmente indipendenti, al massimo (ma ≤ n),
contemporaneamente attivi. Ottimalit`
a della soluzione in variabili primali implica l’ammissibilit`
a in termini di variabili duali [30, 31], che sono calcolate usando il tensore di
Moore–Penrose N ∗a :
µ[λ] = −N ∗a g[∆t]
,
N ∗a = (N Ta EN a )−1 N Ta E
,
g[∆t] ≡ F ∆t
(4.5.1)
L’algoritmo iterativo, deﬁnisce una sequenza di S–coppie (∆σh , Ah) aggiungendo o rimuovendo un vincolo tra quelli violati nel problema da minimizzare, secondo la strategia descritta
in [32] e schematizzata nella tabella 4.6, mantenendo sempre l’ammissibilit`
a delle variabili
duali per i sottoproblemi intermedi ﬁnch`e non `e raggiunta l’ammissibilit`
a delle variabili
primarie per tutti i vincoli del problema. L’iterazione `e inizializzata con la S–coppia vuota
(∆σ = 0, A = ∅) ottenuta a partire dal predittore elastico e considerando nulli i vincoli.
Durante il processo iterativo `e mantenuta l’indipendenza lineare dei vincoli attivi. La
ˆ = 0 (vedi tab. 4.6) signiﬁca che, aggiungendo il corrispondente vincolo al set
condizione σ
attivo A la matrice N presenta colonne linearmente dipendenti. L’algoritmo, che procede
attivando alcuni dei vincoli violati, potrebbe richiede di eliminare alcuni dei vincoli attivi
per assicurare sempre l’ammissibilit`
a duale e l’ndipendenza dei vincoli.
Se all’inizio del Step 1 la soluzione corrente (∆σ(h) , A(h)) del subproblema QP(A(h) )
soddisfa tutte i vincoli K allora la soluzione `e anche soluzione del problema QP(K), ovvero
(∆σ, A) = (∆σ(h) , A(h)), altrimenti con un solo passo, completo o parziale, si riesce a
determinare la nuova soluzione (∆σ(h+1) , A(h+1)) oppure il problema non ha soluzione.
o ripetere. Inoltre, essendo
Poich`e ΠHR (∆σ (h+1)) > ΠHR (∆σ (h) ), la S–coppia non si pu`
le S–coppie in numero ﬁnito, l’algoritmo risolve il problema QP con un numero ﬁnito di
passi. Una descrizione pi`
u accurata del metodo pu`
o essere trovata in [32].
4.6
Cenno alla soluzione del problema senza linearizzare i
vincoli
La linearizzazione del dominio elastico E non `e un requisito essenziale richiesto dall’analisi
a shakedown e il dominio di shakedown Es [λ] si determina allo stesso modo attraverso
le equazioni (2.1.2). Nel caso di superﬁci non–lineari, l’algoritmo di return–mapping pu`
o
essere ottenuto mediante una procedura di ottimizzazione non–lineare, che in genere risolve
il problema come una sequenza di problemi di programmazione quadratica su vincoli lineari
QP [30]. Nonostante la natura non–lineare di questo approccio alternativo, la procedura di
return–mapping risulta molto eﬃciente [33].
Si noti, comunque, che l’approccio presentato, basato su una rappresentazione lineare
o essere adattato in contesti del tutto generici dove
del dominio di adattamento Es [λ], pu`
non si conosce l’espressione analitica della funzione del dominio elastico, come ad esempio,
il caso dei telai in cemento armato in cui il dominio elastico della sezione degli elementi si
pu`
o ricostruire soltanto per punti.
L’algoritmo di rientro sul dominio ammissibile
37
Tabella 4.6. Algoritmo active–set–method
Step 0 Inizializzazione La soluzione iniziale coincide con la soluzione del problema non vincolato:
ΠHK [0] ≡
1
∆σ T F ∆σ = 0;
2
∆σ (0) = 0;
A(0) = ∅;
mA = 0
Step 1 Controllo delle violazioni Si controllano i vincoli violati non contenuti nel set A(h) :
(h)
− bj > 0}
V = {j ∈ K \ A(h) : fj ≡ nT
j ∆t
IF (V = ∅)
ELSE
:
:
STOP. La soluzione corrente `e ammissibile e ottimale.
si sceglie p ∈ V e si pone
µ
, A+ = A ∪ {p}
n + = n p , µ+ =
0
Step 2 (Step Primale/Duale ) Valutazione della nuova S–coppia (∆σ (h) , A(h) ):
(a) Ricerca della direzione
IF mA > 0
ELSE
:
:
ricerca nello spazio primale
ricerca nello spazio duale
:
:
ˆ = E tn+
σ
ˆ = −N ∗ n +
µ
(b) Calcolo della lunghezza del passo s = min(s1 , s2 ) con:
(i) s1 il massimo passo nello spazio duale senza violare l’ammissibilit`
a delle variabili duali:
⎧
ˆ≤ 0 or mA = 0
⎪
⎨ ∞, if µ
+
µ
µ+
s1 =
j
k
, j∈A
=
min
⎪
⎩ µˆj >0 µ
ˆj
µ
ˆk
(ii) s2 il passo minimo nello spazio primale tale che il vincolo p risulta non violato
⎧
ˆ =0
⎨ ∞, if σ
(h)T +
s2 =
n
⎩ bp − ∆σ
ˆ Tn+
σ
(c) Si determina la nuova S–coppia e si compie il passo:
(i) Passo impossibile nello spazio primale e duale: IF (s = ∞) STOP : problema impossibile
(ii) Passo nello spazio duale: s2 = ∞ e s1 = ∞. Assegna:
ˆ
−µ
µ(h) = µ+ + s
;
1
elimina il vincolo k dal set attivo:
A(h) = A(h−1) \ {k};
mA = mA − 1;
aggiorna E t e N a e vai al Step 1.
(iii) Passo nello spazio duale:
ˆ ,
∆σ (h) = ∆σ (h−1) + sσ
ˆ
−µ
µ(h) = µ + + s
1
1
ˆ T n + ( s + µ+
ΠHR = ΠHR + s σ
mA +1 )
2
IF (s = s2 ) (passo completo): assegna µ(h) = µ+ e aggiungi il vincolo p al set attivo:
A(h) = A(h−1) ∪ {p}
,
mA = mA + 1
aggiorna E t e N e vai a Step 1.
IF (s = s1 ) (passo parziale): si elimina il vincolo k e si assegna:
A(h) = A(h−1) \ {k}
,
mA = mA − 1
update E t e N e vai a Step (2.a).
Capitolo 5
Esempio di implementazione: i
telai in c.a.
Il metodo iterativo per l’analisi a shakedown proposto in [10] e discusso nei capitoli precedenti, sar`
a esteso all’analisi di strutture intelaiate spaziali in cemento armato.
Il modello discreto del telaio spaziale sar`a ottenuto utilizzando l’elemento ﬁnito di asta
tridimensionale proposto da Petrolo–Casciaro. Questo elemento `e derivato dalla teoria del
De Saint Ven`
ant e tiene conto dell’interazione taglio–torsione che si genera nell’elemento
spaziale per eﬀetto dell’eccentricit`a dell’azione tagliante rispetto al centro di taglio. La matrice di rigidezza dell’elemento `e ottenuta a partire dalla soluzione numerica delle equazioni
diﬀerenziali che deﬁniscono il problema del solido spaziale di De Saint Ven`
ant.
Il comportamento elastoplastico delle sezioni sar`a modellato mediante una rappresentazione lineare a falde del dominio limite per pressoﬂessione deviata, ottenuta in modo
automatico a partire dalla geometria della sezione e dalle armature presenti [12]. Per ogni sezione si assumono 26 falde, ciascuna delle quali `e deﬁnita in funzione di un possibile
cinematismo plastico della sezione. L’algoritmo `e descritto in dettaglio e sono riportati e
discussi una serie di risultati numerici che mostrano l’eﬃcacia del metodo proposto come
concreto strumento di veriﬁca strutturale.
5.1
L’elemento finito asta 3d
La teoria del De Saint Ven`
ant [34] rappresenta una buona base teorica dalla quale derivare
l’elemento asta da usare nell’analisi ad elementi ﬁniti di telai spaziali, poich`e consente
una accurata descrizione dell’elemento monodimensionale in termini di azioni risultanti
applicate alle sezioni di estremit`a dell’elemento.
La costruzione della matrice di rigidezza dell’elemento richiede la valutazione delle propriet`
a elastiche della sezione trasversale, che rappresentano i fattori associati alle componenti di sforzo in funzione della geometria della sezione trasversale dell’elemento. Mentre
i fattori associati allo sforzo normale e al momento ﬂettente sono ottenuti facilmente a
partire dalla nota distribuzione delle tensioni normali, i fattori torsionali e di taglio, che
derivano dalla distribuzione delle tensioni tangenziali, richiedono la soluzione di un sistema di equazioni diﬀerenziali deﬁnite nel domino della sezione, in genere di non facile
determinazione.
38
Esempio di implementazione: i telai in c.a.
39
In genere, il contributo all’energia di deformazione associato al taglio pu`
o essere considerato piccolo. Inoltre, non `e necessaria una valutazione accurata delle tensioni di taglio.
Ci`o giustiﬁca l’uso di formule approssimate, come quelle dedotte da Bredt o dalla teoria di
Jourawski, per una variet`
a di geometrie della sezione, per cui la soluzione delle equazioni differenziali non `e necessaria. Comunque, questa approssimazione implica una specializzazione
in relazione al particolare problema che deve essere analizzato ed `e possibile solamente se si
hanno a priori le informazioni qualitative sulla distribuzione delle tensioni dovute al taglio.
Di conseguenza, questo modo di procedere `e conveniente in analisi “eseguite a mano”, ma
non `e adatto ad essere inserito in algoritmi numerici del tutto generici.
Nella presente formulazione si mostra come la teoria del De Saint Ven`ant rappresenta
un valido strumento per perfezionare una procedura automatica appropriata sia per deﬁnire
la matrice di rigidezza dell’elemento, sia per recuperare i valori delle tensioni tangenziali
interni alla sezione.
5.1.1
La teoria di De Saint Ven`
ant
Si consideri il solido cilindrico isotropico di Cauchy soggetto alle forze di superﬁcie applicate
alle sezioni di estremit`a, mostrato in ﬁgura 5.1. Il cilindro `e riferito ad un sistema locale
{x, y, z} di assi cartesiani baricentrico, orientato secondo le direzioni principali d’inerzia
della sezione trasversale contenuta nel piano {x, y}. Si indicano, inoltre, con Ω il dominio
ˆ la normale esterna al contorno di componenti nx , ny .
della sezione di contorno Γ e con n
Figura 5.1: Il solido cilindrico di De Saint Ven`ant
Le assunzioni di De Saint Ven`ant implicano σ xx = σ yy = σ xy = 0 (v. [34]). Combinando
con le equazioni di Cauchy (equilibrio)
σ,ij
j =0
,
i, j = x, y, z
(5.1.1)
e con le equazioni di Beltrami–Mitchell (compatibilit`
a cinematica)
σ,ij
kk +
1
σ,kk = 0
1 + ν ij
,
i, j, k = x, y, z
(5.1.2)
con ν fattore di Poisson, consentono di esprimere la tensine normale σ zz nella forma
σ zz = a + a1 x + a2 y − z(b + b1 x + b2 y)
(5.1.3)
Esempio di implementazione: i telai in c.a.
40
in cui a · · · b2 sono costanti di integrazione da determinare, mentre le tensioni tangenziali
σ zx e σ zy devono soddisfare il sistema di equazioni diﬀerenziali
σ zx ,x +σ zy ,y = b + b1 x + b2 y
ν
, ν¯ :=
(5.1.4)
1+ν
σ zy ,x −σ zx ,y = ν¯(b1 y − b2 x) + c
con c ulteriore costante di integrazione da determinare. La soluzione generale del sistema
(5.1.4) pu`
o essere espressa nella forma
σ zx
σ zx = w,x +¯
(5.1.5)
σ zy = w,y +¯
σ zy
in cui
1
σ
¯ zx := [bx + b1 (x2 − ν¯y 2 ) − cy] ,
2
1
σ
¯ zy := [by + b2 (y 2 − ν¯x2 ) + cx]
2
(5.1.6)
`e una soluzione particolare delle (5.1.4) e w := w[x, y] `e una funzione ausiliaria che soddisfa
l’equazione di Laplace
(5.1.7)
w,xx +wyy = 0 , {x, y} ∈ Ω
la cui condizione al contorno `e recuperata dall’equazione di equilibrio di Cauchy σ zx nx +
σ zy ny = 0 valida sulla superﬁcie laterale del cilindro:
σ zx nx + σ
¯ zy ny ) ,
w,n := w,x nx + w,y ny = −(¯
{x, y} ∈ Γ
(5.1.8)
Le equazioni (5.1.3), (5.1.5)–(5.1.8) deﬁniscono i campi di tensione a meno delle costanti
a · · · c che possono essere recuperate dalle condizioni di equilibrio globale
zz
zz
σ dΩ = N0 ,
σ ydΩ = Mx0 + zTy ,
σ zz xdΩ = −My0 + zTy
(5.1.9)
Ω
Ω
Ω
(σ zx y − σ zy x)dΩ = Mt
(5.1.10)
Ω
Si ottiene
a=
N
,
A
a1 = −
My0
Ty
Mx0
Tx
, a2 =
, b = 0 , b1 = −
, b2 = −
Jy
Jx
Jy
Jx
(5.1.11)
mentre, la costante c `e deﬁnita in forma implicita dalla (5.1.10). Le tensioni normali
σ zz sono espresse in forma esplicita dalla (5.1.3), mentre le tensioni tangenziali σ zx e σ zy
sono ottenute dalle equazioni (5.1.5) e richiedono la soluzione del problema diﬀerenziale
di Neumann (5.1.7), (5.1.8), di cui la soluzione esiste ed `e unica, a meno di una costante
inessenziale, sotto la condizione
w,n ds := − (¯
σ zx nx + σ
¯ zy ny )ds = 0
(5.1.12)
Γ
Γ
che nel nostro caso `e soddisfatta identicamente. Usando la notazione matriciale, la soluzione
delle equazioni di De Saint Ven`
ant, pu`
o essere scritta nella seguente forma compatta
⎧
My [z]
Mx [z]
N
⎨
−
x+
y
σ=
(5.1.13)
A
Jy
Jx
⎩
τ = ∇w − Tx b1 − Ty b2 − c b3
Esempio di implementazione: i telai in c.a.
41
in cui ∇ `e l’operatore gradiente,
σ := σ
zz
,
τ :=
σ zx
σ zy
(5.1.14)
N , Mx , My , Tx e Ty sono gli sforzi normale, ﬂettenti e di taglio che agiscono sulla sezione
trasversale, A, Jx , Jy sono l’area ed i momenti d’inerzia della sezione ed i vettori b1 , b2 e
b3 sono deﬁniti da
2
1
1
1
y
x − ν¯y 2
0
(5.1.15)
, b2 :=
, b3 :=
b1 :=
0
2Jy
2Jx y 2 − ν¯x2
2 −x
La funzione ausiliaria w espressa nella forma
w[x, y] := Txw1 [x, y] + Ty w2 [x, y] + cw3 [x, y]
ottenuta dalla soluzione dei tre problemi diﬀerenziali
{x, y} ∈ Ω
wj ,xx +wj ,yy = 0
ˆ
{x, y} ∈ Γ
wj ,n = bTj n
,
j = 1, 2, 3
consente di esplicitare la costante c deﬁnita dalla (5.1.10), ottenendo:
1
c = − (r1 Tx + r2 Ty − Mt ) , rj = 2 {bj − ∇wj }T b3 dΩ
r3
Ω
I campi di tensione possono essere riscritti in forma esplicita nella forma:
⎧
My [z]
Mx [z]
N
⎨
−
x+
y
σ=
A
Jy
Jx
⎩
τ = Dt
(5.1.16)
(5.1.17)
(5.1.18)
(5.1.19)
in cui t := {Tx, Ty , Mt}T e la matrice D := [d1 | d2 | d3 ] `e deﬁnita dai vettori colonna
d1 := ∇w1 − b1 − r1 d3
d2 := ∇w2 − b2 − r2 d3
d3 := (∇w3 − b3 )/r3
(5.1.20)
` stata ottenuta una soluzione del problema di De Saint Ven`
E
ant che consente di esprimere le
tensioni che si generano in ogni punto della sezione in forma esplicita, a meno delle funzioni
u approfondita del Principio
wj deﬁnite da un problema diﬀerenziale. Per una discussione pi`
di De Saint Ven`
ant si pu`
o consultare [35].
5.1.2
Generazione dell’elemento asta 3d
L’elemento ﬁnito di asta spaziale `e basato sulla teoria del De Saint Ven`
ant illustrata nella
sezione precedente. L’elemento `e descritto assegnando la lunghezza dell’asta, la geometria
della sezione e l’orientamento nello spazio. Quest’ultimo sar`
a assegnato deﬁnendo i versori
{ix , iy , iz } del sistema locale {x, y, z} nel sistema cartesiano globale {x1 , x2 , x3 }, mostrato
Esempio di implementazione: i telai in c.a.
42
Figura 5.2: Geometria dell’elemento
in ﬁgura 5.2. Gli estremi dell’elemento sono rappresentati dai nodi i e j, deﬁnendo con di
e dj la diﬀerenza di posizione tra il baricentro della sezione ed il nodo corrispondente.
La meccanica dell’elemento `e descritta dalle forze f i , f j e dai momenti mi , mj applicati
ai nodi estremi i e j dell’elemento e dagli spostamenti ui , uj e rotazioni φi , φj nodali
associati:
⎧
⎧
⎫
⎫
fi ⎪
ui ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎨
⎨
⎬
⎬
mi
φi
, ue :=
(5.1.21)
f e :=
f ⎪
u ⎪
⎪
⎪
⎪
⎪
⎩ j ⎪
⎩ j⎪
⎭
⎭
mj
φj
il lavoro esterno `e deﬁnito come:
f Te ue
(5.1.22)
mentre la matrice di rigidezza K e [12 × 12] dell’elemento, deﬁnita implicitamente attraverso
la relazione di Clapeiron
1
1
(5.1.23)
Φe = f Te ue = uTe K e ue
2
2
sar`
a ottenuta recuperando la relazione esplicita dell’energia di deformazione Φe .
Descrizione locale
Indicando con
f l = {fxi , fyi , fzi , mxi, myi , mzi , fxj , fyj , fzj , mxj , myj , mzj }T
ul = { uxi , uyi , uzi , ϕxi , ϕyi , ϕzi , uxj , uyj , uzj , ϕxj , ϕyj , ϕzj }T
(5.1.24)
i vettori che raccolgono, rispettivamente, forze e spostamenti nodali riferiti ad un sistema
locale baricentrico {x, y, z} orientato secondo le direzioni principali d’inerzia della sezione
dell’elemento, si ha:
(5.1.25)
f Tl ul = f Te ue
Esempio di implementazione: i telai in c.a.
43
Figura 5.3: Componenti locali delle azioni esterne
La relazione tra le componenti globali e le componenti locali `e ottenuta da semplici considerazioni geometriche (ved. ﬁg. 5.3). Indicando con il simbolo ∧ il prodotto vettoriale
esterno, si ha:
⎧
⎧
ϕxi = iTx φi
uxi = iTx ui + (φi ∧ di )T ix
⎪
⎪
⎪
⎪
⎪
⎪
T
T
⎪
⎪
uyi = iy ui + (φi ∧ di ) iy
ϕyi = iTy φi
⎪
⎪
⎪
⎪
⎨
⎨
T
T
uzi = iz ui + (φi ∧ di ) iz
ϕyi = iTz φi
,
(5.1.26)
uxj = iTx uj + (φj ∧ dj )Tix
ϕxj = iTx φj
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
u = iTy uj + (φj ∧ dj )T iy
ϕyj = iTy φj
⎪
⎪
⎩
⎩ yj
uzj = iTz uj + (φj ∧ dj )Tiz
ϕyj = iTz φj
La matrice di ﬂessibilit`
a della sezione
La matrice di ﬂessibilit`
a della sezione `e deﬁnita dal bilancio energetico
1
1
1 T
2
s Hs :=
σ dΩ +
τ T τ dΩ
2
2E Ω z
2G Ω
(5.1.27)
in cui E e G = E/2(1 + ν) sono, rispettivamente, il modulo elastico normale e tangenziale,
s := {N, Mx , My , Tx , Ty , Mt }T
(5.1.28)
`e il vettore che raccoglie le 6 componenti di sollecitazione che agiscono nella sezione. Usando
l’espressione ottenuta nella sezione 5.1.1, si ha:
"
!
Hf 0
(5.1.29)
H=
0 Ht
a ﬂessionale e tangenziale
in cui H f e H t sono, rispettivamente, la matrice di ﬂessibilit`
deﬁnite come:
⎡
⎤
1/A 0
0
1 ⎣
1
⎦
0 1/Jx 0
DT D dΩ
(5.1.30)
, H t :=
H f :=
E
G Ω
0
0 1/Jy
con la matrice D deﬁnit`
a dall’equazione (5.1.20).
Esempio di implementazione: i telai in c.a.
44
Componenti naturali di forze e spostamenti
L’equilibrio della trave implica
fyi + fyj = 0
fzi + fzj = 0
fxi + fxj = 0
mxi + mxj + fyi = 0 myi + myj − fxi = 0 mzi + mzj = 0
(5.1.31)
mentre, all’ascissa z, le sollecitazioni interne alla trave sono:
N [z] := fzj
Mx [z] := [mxi( − z) − mxj z]/
My [z] := [myi ( − z) − myj z]/
,
Tx[z] := (mxj + mxi )/
Ty [z] := (myj + myi )/
M t[z] :=
mzj
(5.1.32)
` conveniente esprimere il vettore delle sollecitazioni s[z] come una combinazione di 6 modi
E
naturali di tensione [36], tale che:
s[z] = S[z]T f n
in cui la matrice S ed il vettore delle tensioni naturali f n sono deﬁnite da:
⎧
⎡
⎤
fzj
1
·
·
·
· ·
⎪
⎪
⎪
⎪
⎢·
⎥
m
− mxj
1/2
·
·
·
·
⎪
xi
⎪
⎢
⎥
⎨
⎢·
⎥
m
·
1/2
·
·
·
yi − myj
⎥ , f n :=
S[z] := ⎢
⎢ · 1/2 − z/
⎥
mxi + mxj
·
· 1/ · ⎥
⎪
⎪
⎢
⎪
⎪
⎣·
m + myj
⎪
·
1/2 − z/ 1/ · · ⎦
⎪
⎩ yi
mzj
·
·
·
·
· 1
(5.1.33)
⎫
⎪
⎪
⎪
⎪
⎪
⎪
⎬
⎪
⎪
⎪
⎪
⎪
⎪
⎭
(5.1.34)
Dalla condizione di dualit`
a
Figura 5.4: Modi naturali di deformazione e tensioni
f Te ue = f Tn un
(5.1.35)
Esempio di implementazione: i telai in c.a.
si ottiene il vettore un delle deformazioni naturali associate a f n :
⎧
⎫
uzj − uzi
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
(ϕ
−
ϕ
)/2
⎪
⎪
xi
xj
⎪
⎪
⎨
⎬
(ϕyi − ϕyj )/2
un :=
(ϕxi + ϕxj )/2 + (uyj − uyi )/ ⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
(ϕyi + ϕyj )/2 + (uxi − uxj )/ ⎪
⎪
⎪
⎩
⎭
ϕzj − ϕzi
45
(5.1.36)
che, usando la (5.1.26), pu`
o essere espresso in termini di ue dalla relazione matriciale
un = T ne ue
in cui
(5.1.37)
⎡
T ne
⎤
−iTz
−(di ∧ iz )T
iTz
(dj ∧ iz )T
⎢ ·
⎥
iTx /2
·
−iTx /2
⎢
⎥
T
T
⎢ ·
⎥
/2
·
−i
/2
i
y
y
⎢
⎥
:= ⎢ T
T
T
T
T
T
⎥
/
i
/2
−
(d
∧
i
)
/
i
/
i
/2
+
(d
∧
i
)
/
−i
i
y
j
y
x
y
x
⎢ y
⎥
T
T
T
T
T
T
⎣ i / i /2 + (di ∧ ix) / −i / i /2 − (dj ∧ ix ) / ⎦
x
y
x
y
·
iTz
·
−iTz
(5.1.38)
Matrice di ﬂessibilit`
a e di rigidezza dell’elemento
Combinando le equazioni (5.1.34, 5.1.33) e (5.1.27, 5.1.29, 5.1.30), si ottiene l’espressione
dell’energia di deformazione dell’elemento:
1 T
1
s Hs dz = f Tn F n f n
(5.1.39)
Φe :=
2 0
2
in cui la matrice di ﬂessibilit`
a dell’elemento F n `e ottenuta, in forma esplicita, come:
⎤
⎡
h11 0
0
0
0
0
⎢
0
0
0
0 ⎥
h22 /4
⎥
⎢
⎢
0
0
0 ⎥
h33 /4
⎥
⎢
(5.1.40)
Fn = ⎢
h45 /
h56 ⎥
h22 /12 + h55 /
⎥
⎢
⎣
sym
h33 /12 + h44 / h46 ⎦
h66 Introducendo la matrice naturale di rigidezza K n deﬁnita da
K n := F −1
n
(5.1.41)
si ottiene, dalle equazioni (5.1.39) e (5.1.23), la relazione
1
1
Φe := uTe K e ue = uTn K n un
2
2
(5.1.42)
che consente di ricavare l’espressione della matrice di rigidezza K e dell’elemento in forma
esplicita:
(5.1.43)
K e = T Tne K n T ne
Esempio di implementazione: i telai in c.a.
5.1.3
46
Soluzione numerica delle equazione di De Saint Ven`
ant
La valutazione numerica della matrice di rigidezza dell’elemento tramite la (5.1.43), richiede
che siano note le componenti hij della matrice di ﬂessibilit`
a H della sezione trasversale
dell’asta, deﬁnita dall’equazione (5.1.29). Ci`
o implica la determinazione delle sub–matrici
H f e H t.
Il calcolo della matrice H f richiede la valutazione delle sole caratteristiche geometriche
2
dA , Jx :=
y dA , Jy :=
x2 dA
(5.1.44)
A :=
Ω
Ω
Ω
della sezione, che sono facilmente ottenibili, attraverso integrazione per parti, con la seguente
formula generale:
h
k
x y dΩ =
Ω
NL j=1
Γj
(xh+1 y k )ˆ
njx + (xh y k+1 )ˆ
njy
dΓ
h+k+2
(5.1.45)
ˆ jx e n
ˆ jx sono le componenti
in cui NL sono il numero di lati Γj del contorno della sezione, n
della normale esterna al lato Γj .
Il calcolo della matrice H t attraverso la (5.1.30) e (5.1.20), richiede la valutazione
numerica delle tre funzioni scalari wj [x, y] deﬁnite implicitamente dai problemi diﬀerenziali
(5.1.17). La soluzione dei problemi diﬀerenziali, per sezioni generiche, non pu`
o essere
ottenuta in forma chiusa, bens`ı soltanto attraverso una discretizzazione numerica.
La soluzione numerica `e stata ottenuta sia con il metodo degli elementi ﬁniti.
Soluzione agli elementi finiti
La ricerca della soluzione numerica con tecnica FEM richiede che il dominio della sezione
venga suddiviso in elementi [37], all’interno dei quali la funzione wj [x, y] `e interpolata in
funzione dei valori che la stessa assume ai nodi, nella forma
wj [x, y] = Ψ qe
,
q e := Ae q
(5.1.46)
in cui Ψ `e la matrice delle funzioni interpolanti e q e `e il vettore che raccoglie i valori nodali
della wj [x, y] che possono essere recuperate dal vettore globale q dei parametri nodali
attraverso l’operatore Ae .
La soluzione `e ottenuta minimizzando, nel dominio dei valori nodali ammissibili q, il
funzionale
1
2
2
ˆ
(wj ,x +wj ,y )dΩ − pj wj dΓ , pj := bTj n
(5.1.47)
Π[wj ] =
2 Ω
Γ
associato al problema diﬀerenziale (5.1.17). La stazionariet`a del funzionale (5.1.47) implica:
(qTe M e δqe − pe δqe ) = 0 , ∀ δqe
(5.1.48)
e
in cui la somma `e estesa a tutti gli elementi, mentre la matrice M e e il vettore pe sono
deﬁnite da:
T
T
( Ψ,x Ψ,x + Ψ,y Ψ,y ) dΩ , pe :=
pj Ψ dΓ
(5.1.49)
M e :=
Ωe
Γe
Esempio di implementazione: i telai in c.a.
47
Usando la procedura standard di assemblaggio, la (5.1.47) in forma discreta, fornisce il
seguente sistema lineare
ATe M e Ae , p :=
ATe pe
(5.1.50)
M q = p , M :=
e
e
la cui soluzione, per ognuno dei tre problemi diﬀerenziali deﬁniti dalla (5.1.17), permette
la valutazione numerica delle funzioni wj [x, y] necessarie al calcolo dei vettori colonna di
della matrice D. I coeﬃcienti htij della matrice di ﬂessibilit`
a H t sono ottenuti direttamente
dalla (5.1.30):
1
dT dj dΩ
(5.1.51)
htij =
G Ω i
in cui ciascuno dei vettori dk e deﬁnito dalla (5.1.20). L’integrale `e valutato mediante una
procedura numerica standard alla Gauss.
L’interesse a ricavare coeﬃcienti integrali da valutare su tutta la sezione, suggerisce
l’uso di elementi ﬁniti semplici tale da sempliﬁcare il pi`
u possibile il processo di soluzione
numerica. In particolare, sono stati utilizzati elementi triangolari standard a tre nodi ed
elementi triangolari a sei nodi. Questi ultimi, basati su una interpolazione quadratica della
funzione w[x, y], rappresentano un buon compromesso tra accuratezza del risultato ottenuto
ed eﬃcienza computazionale.
Elemento triangolare a tre nodi – T3
Si assume all’interno dell’elemento una interpolazione lineare per la funzione wj [x, y]:
Ψ=
in cui
1 α1 + β1 x + γ1 y
2Ae
α2 + β2 x + γ2 y
α3 + β3 x + γ3 y
⎧
⎧
⎧
⎨ α1 = x2 y3 − x3 y2 ⎨ α2 = x3 y1 − x1 y2 ⎨ α3 = x1 y2 − x2 y1
β = y2 − y3
β = y3 − y1
β = y1 − y2
⎩ 2
⎩ 3
⎩ 1
γ1 = x3 − x2
γ2 = x1 − x3
γ3 = x2 − x1
(5.1.52)
(5.1.53)
Ae `e l’area dell’elemento e {xi , yi} (i = 1, 2, 3) sono le coordinate dei vertici (v. ﬁg. 5.5).
Figura 5.5: Elemento T3 ed elemento T6
Esempio di implementazione: i telai in c.a.
Dall’equazione (5.1.49) si ottiene:
⎡ 2
β1 + γ12
Me = ⎣
sym.
48
⎤
β1 β3 + γ 1 γ 3
β2 β3 + γ 2 γ 3 ⎦
β32 + γ32
β1 β2 + γ 1 γ 2
β22 + γ22
(5.1.54)
mentre i vettori pe1 , pe2 e pe3 , associati rispettivamente alle tre condizioni al contorno w1,n,
w2,n e w3,n , deﬁnite dalle equazioni (5.1.17), sono ottenuti sommando i contributi dei lati
dell’elemento. Per il generico lato di nodi 1–2, si ha:
2
2
2
2
y2 − y1 2(x1 − ν¯y1 ) + (x1 + x2 ) − ν¯(y1 + y2 )
p1
=
(5.1.55)
pe1 =
p2
24Jy
2(x22 − ν¯y22 ) + (x1 + x2 )2 − ν¯(y1 + y2 )2
2
2
2
2
x1 − x2 2(y1 − ν¯x1 ) + (y1 + y2 ) − ν¯(x1 + x2 )
p1
=
(5.1.56)
pe2 =
p2
24Jx
2(y22 − ν¯x22 ) + (y1 + y2 )2 − ν¯(x1 + x2 )2
2
x2 + y22 + x1 x2 + y1 y2 − 2x21 − 2y12
1
p1
=
(5.1.57)
pe3 =
p2
12 −x21 − y12 − x1 x2 − y1 y2 + 2x22 + 2y22
I contributi dei lati di nodi 2–3 e 3–1 sono ottenuti in modo analogo, operando una permutazione ciclica degli indici. Si pu`
o osservare, che per la continuit`
a della funzione wj [x, y],
risultano diversi da zero solo i contributi che provengono dai lati che stanno sul contorno
della sezione.
Elemento triangolare a sei nodi – T6
Si assume una interpolazione quadratica per la funzione wj [x, y]:
ψ1 = 2ξ12 − ξ1
ψ 4 = 4 ξ2 ξ3
in cui
ξi =
ψ2 = 2ξ22 − ξ2
ψ 5 = 4 ξ3 ξ1
ψ3 = 2ξ32 − ξ3
ψ 6 = 4 ξ1 ξ2
(5.1.58)
i = 1, 2, 3.
(5.1.59)
1
(αi + βi x + γiy) ,
2Ae
con Ae , αi , βi e γi deﬁnite nella sezione precedente. Dall’equazione (5.1.49) si ottiene:
⎡
6k1 −k12 −k13
·
⎢
6k
−k
4k
2
23
23
⎢
1 ⎢
6k
4k
3
23
⎢
Me =
8k
12 ⎢
123
⎢
⎣ sym.
4k13
·
4k13
8k12
8k123
⎤
4k12
4k12 ⎥
⎥
· ⎥
⎥
8k13 ⎥
⎥
8k23 ⎦
8k123
(5.1.60)
con i coeﬃcienti kα deﬁniti da
k1 =
β12 + γ12
β 2 + γ22
β 2 + γ32
, k2 = 2
, k3 = 3
, k123 = k1 + k2 + k3
2Ae
2Ae
2Ae
k12 = k3 − k1 − k2 , k13 = k2 − k1 − k3 , k23 = k1 − k2 − k3 ,
(5.1.61)
Esempio di implementazione: i telai in c.a.
49
I vettori pe1 , pe2 e pe3 , associati alle tre condizioni al contorno dei problemi deﬁniti dalla
(5.1.17), per il generico lato di nodi 1–2, sono:
⎫
⎧ ⎫
⎧
⎨ p1 ⎬ y − y ⎨ 9x21 − x22 + 2x1 x2 + ν¯(y22 − 9y12 − 2y1 y2 ) ⎬
2
1
9x22 − x21 + 2x1 x2 + ν¯(y12 − 9y22 − 2y1 y2 )
pe1 = p2 =
(5.1.62)
⎭
⎩ ⎭
120Jy ⎩
p6
4[3x21 + 3x22 + 4x1 x2 − ν¯(3y12 + 3y22 + 4y1 y2 )]
⎫
⎧ ⎫
⎧
⎨ p1 ⎬ x − x ⎨ 9y12 − y22 + 2y1 y2 + ν¯(x22 − 9x21 − 2x1 x2 ) ⎬
1
2
(5.1.63)
9y22 − y12 + 2y1 y2 + ν¯(x21 − 9x22 − 2x1 x2 )
pe2 = p2 =
⎭
⎩ ⎭
120Jx ⎩
2
2
2
2
p6
4[3y1 + 3y2 + 4y1 y2 − ν¯(3x1 + 3x2 + 4x1 x2 )]
⎫
⎧ 2
⎧ ⎫
x1 + y12 − x1 x2 − y1 y2 ⎬
⎨
⎨ p1 ⎬
1
(5.1.64)
x2 + y 2 − x1 x2 − y1 y2
pe3 = p2 =
⎭
⎩ ⎭ 12 ⎩ 2 2 2 2
p6
2(x2 − x1 + y22 − y12 )
I contributi dei lati di nodi 2–3–4 e 3–1–5 sono ottenuti in modo analogo, operando una
permutazione ciclica degli indici.
Valutazione numerica della matrice di ﬂessibilit`
a Ht
Note le funzioni wj [x, y], i coeﬃcienti della matrice htij della matrice H t sono ottenuti
mediante integrazione numerica alla Gauss:
1
1 dTi dj dΩ =
di [xeg , yeg ]Tdj [xeg , yeg ]wg
(5.1.65)
htij =
G Ω
G e g
I vettori dk sono dati dall’equazione (5.1.20) e la somma `e estesa a tutti gli elementi e e
a tutti i punti di Gauss g interni all’elemento. L’integrazione della funzione, quadratica
per l’elemento T3 e quartica per l’elemento T6, `e realizzata mediante 3 e 9 punti di Gauss,
rispettivamente.
5.1.4
Analisi di sezioni in campo elastico
Sono stati eseguiti una serie di test numerici, determinando le caratteristiche elastiche di
sezioni generiche necessarie alla costruzione della matrice di rigidezza dell’alemento asta
3d. In seguito saranno presentati tre gruppi di risultati. Nel primo sono riportati dei test
adatti alla validazione del metodo proposto, in quanto la soluzione numerica ottenuta e
confrontata con analoga soluzione analitica o valori numerici di riferimento. Nel secondo
gruppo di risultati sar`
a mostrata la convergenza dei risultati numerici al variare della discretizzazione del modello, valutandone anche l’eﬃcienza computazionale delle strategie di
soluzione proposte. Inﬁne saranno presentati i risultati ottenuti per alcune sezioni generiche
di interesse tecnico.
Per rendere pi`
u agevole il confronto con valori di riferimento, la soluzione ottenuta `e
presentata sotto forma di coordinate del centro di taglio xc e yc , di orientamento del sistema
principale del taglio rappresentato dall’angolo αt che esso forma con il sistema principale
d’inerzia e in termini di fattori di taglio e torsione, χ1 , χ2 e χt , che intervengono nella
scrittura dell’energia di deformazione associata al taglio nella forma tradizionale:
2
¯ t2 1 T¯12
T¯22
M
τ
1
dΩ =
χ1 +
χ2 +
χt
(5.1.66)
2 Ω G
2 GA
GA
GJp
Esempio di implementazione: i telai in c.a.
50
in cui
T¯1 := Tx cos αt + Ty sin αt
,
T¯2 := −Tx sin αt + Ty cos αt
sono le componenti principali dello sforzo di taglio,
Jp := (x2 + y 2 )dΩ
(5.1.67)
(5.1.68)
Ω
rappresenta il momento di inerzia polare della sezione,
¯ t := Mt + Tx yc − Ty xc
M
(5.1.69)
`e il momento torsionale riferito al centro di taglio.
a tangenziale, le quantit`
a appena
Indicando con htij le componenti della matrice di ﬂessibilit`
introdotte sono date da:
xc = −
ht23
ht33
χ1 := GA(c1 + c2 + c3 ) ,
in cui
yc =
ht13
ht33
αt =
2 c3
1
arctan
2
c2 − c1
χ2 := GA(c1 + c2 − c3 ) ,
χt := GJp ht33
⎧
c1 :=(ht11 − 2ht13 yc + ht33 yc2 )/2
⎪
⎪
⎨
c2 :=(ht22 + 2ht23 xc + ht33 x2c )/2
)
⎪
⎪
⎩ c := (ht + ht x − ht y − ht x y )2 + (c − c )2
3
1
2
12
13 c
23 c
33 c c
(5.1.70)
(5.1.71)
(5.1.72)
Test di validazione
Sono stati analizzati i due test semplici mostrati in ﬁgura 5.6. I valori numerici ottenuti per
i parametri elastici della sezione sono riportati nelle tabelle 5.1 e 5.2. Tali valori coincidono
con i valori analitici esatti (sezione rettangolare: χ1 = χ2 = 5/6 per ν = 0 e χt ≈ 1.82204,
per ogni ν; sezione triangolare: χ2 ≈ 1.2444 per ν = 0.5) e con i valori numerici ottenuti
da Schramm et al. [38] per la sezione rettangolare.
Figura 5.6: Geometria delle sezioni analizzate
Esempio di implementazione: i telai in c.a.
ν = 0.0
1.20000
1.20000
χ1
χ2
χt
xc = 0 yc = 0
ν = 0.3
1.27479
1.20056
1.82204
ν = 0.5
1.35605
1.20118
Tabella 5.1: Sezione rettangolare
51
ν = 0.0
1.43949
1.24260
ν = 0.3 ν = 0.5
χ1
1.46738 1.49767
1.24348 1.24444
χ2
2.11447
χt
xc = 0 yc = −0.17962
Tabella 5.2: Sezione triangolare
Test di confronto
Sono presentati due test, relativi alla sezione trapezoidale e alla sezione a “C” mostrate in
ﬁgura 5.7 e ﬁgura 5.8, con l’obiettivo di evidenziare l’eﬃcienza dei metodi di soluzione proposti nel caso di sezioni tozze e di sezioni a pareti sottili. In entrambi i casi, i valori numerici
di riferimento dei parametri elastici sono ottenuti numericamente, usando discretizzazioni
molto ﬁtte. Per la sezione a “C”, la teoria approssimanta di Jourawski fornisce valori che
diﬀeriscono poco dalla soluzione numerica xc = 31.078 e χ2 = 2.3880.
ν = 0.0
χ1
1.34681
χ2
1.18408
61.30030
αT
xc = 1.63713; yc
ν = 0.3
ν = 0.5
1.37708
1.41004
1.18556
1.18709
62.15280
62.82870
= 1.38895; χt = 1.64811
Figura 5.7: Sezione tozza trapezoidale
ν = 0.0
ν = 0.3
ν = 0.5
χ1
2.86230 2.87395 2.88346
2.27017 2.27033 2.27024
χ2
xc = 30.3777; χt = 59.0865
Figura 5.8: Sezione a pareti sottile a “C”
Le ﬁgure 5.9 e 5.12 mostrano, nel caso di ν = 0.3, l’errore di discretizzazione in funzione
Esempio di implementazione: i telai in c.a.
52
Figura 5.9: Errore di discretizzazione e numero di variabili
Figura 5.10: Errore di discretizzazione e tempo di
Figura 5.11: Tempo relativo di analisi
analisi
del numero di variabili del modello discreto. Risulta evidente la migliore prestazione dell’elemento T6 rispetto all’elemento T3. Si noti, comunque, che valori accurati utilizzabili
in ambito tecnico (errore < 1%) si ottengono anche con discretizzazioni abbastanza rade,
specie nel caso della sezione compatta.
La prestazione degli elementi ﬁniti utilizzati `e valutata anche in termini di eﬃcienza
computazionale, facendo riferimento al tempo di analisi. Le ﬁgure 5.10 e 5.13, che mostrano
l’errore di discretizzazione in funzione del tempo impiegato per l’analisi, evidenziano la
prestazione migliore dell’elemento T6 rispetto all’elemento T3. Si noti, inoltre, che una
parte non trascurabile del tempo di analisi richiesto dalla soluzione FEM `e impiegato per
generare la mesh di dominio, come mostrato in ﬁgura 5.10 e ﬁgura 5.13. La percentuale
di tempo richiesto in questa fase preliminare dipende essenzialmente dal numero totale di
elementi e dalla complessit`a del dominio da tassellare. Ci`o penalizza fortemente l’elemento
T3 rispetto all’elemento T6, poich`e quest’ultimo richiede mesh pi`
u rade.
Esempio di implementazione: i telai in c.a.
53
Figura 5.12: Errore di discretizzazione e numero di variabili
Figura 5.13: Errore di discretizzazione e tempo di
analisi
Figura 5.14: Tempo relativo di analisi
Esempio di implementazione: i telai in c.a.
54
Alcuni test di interesse tecnico
In questa sezione sono riportati ulteriori test che mostrano come il metodo proposto si
presta bene a determinare le caratteristiche elastiche di sezioni generiche. In tutti i casi
analizzati `e stato assunto ν = 0.3. L’analisi `e stata eﬀettuata con l’elemento T6, usando
una mesh relativamente ﬁtta. Alcuni dei valori ottenuti sono stati confrontati con analoghi
di riferimento.
Yc
1.567
1.569
[39]
[40]
metodo
proposto
χ1
4.3290
1.569
χ2
χt
1.6686 4.4490
Figura 5.15: Sezione da ponte. Geometria e caratteristiche elastiche
Yc/b
a/b
0.00
0.25
0.50
0.75
esatto
0.5111
0.6327
0.8668
1.0922
[41]
0.5067
0.6287
0.8562
1.0624
metodo
proposto
0.5089
0.6301
0.8650
1.0912
χ1
1.30524
1.43957
1.82390
2.34377
χ2
1.16744
1.35179
1.70015
1.93492
Figura 5.16: Sezione curva a parete sottile. Geometria e caratteristiche elastiche
xc
yc
αT
χ1
χ2
χt
79.465
20.534
-45.00
2.12432
1.52215
5.11438
Figura 5.17: Geometria della sezione e caratteristiche elastiche
χt
1.6887
2.3233
5.2342
23.747
Esempio di implementazione: i telai in c.a.
5.2
55
Definizione del dominio elastico delle sezioni
La deﬁnizione del dominio elastico della sezione degli elementi pressoinﬂessi ricade nel caso
di funzioni generiche di cui non si conosce la funzione analitica che lo descrive. Si ricorre
quindi ad una sua ricostruzione per punti che ne consente una descrizione linearizzata del
dominio.
Si considera una sezione in cemento armato, di forma generica e distribuzione arbitraria
delle armature, deﬁnita dalle coordinate geometriche dei vertici e dalle posizioni e diametri
delle armature. Assunto un comportamento elastoplastico sia per il calcestruzzo che per
l’acciaio, il dominio elastico della sezione pu`
o essere deﬁnito come dominio di interazione
fra sforzo normale N e le due componenti di momento ﬂettente Mx ed My .
o essere determinato a parOgni punto {σ k , nk } della superﬁcie di snervamento pu`
tire da un possibile meccanismo plastico deﬁnito in termini di deformazione assiale ε0 e
deformazione ﬂessionale χx e χy , dal vettore
εpk := {χx , χy , ε0 }T
(5.2.1)
che consente di individuare in modo univoco la posizione dell’asse neutro ed il contorno della
sezione reagente Ωr . Nell’ipotesi di conservazione delle sezioni piane, la deformazione longitudinale ε[x] all’interno della sezione varia linearmente ed `e espressa tramite l’equazione
di Bernoulli–Navier:
ε[x] = ε0 + xχx − yχy
∀x := {x, y}T ∈ Ω
(5.2.2)
Le corrispondenti sollecitazioni limite σk := {Mxk , Myk , Nk }T possono essere calcolate dalle
equazioni di equilibrio
⎧
nb
⎪
⎪
⎪
σc [ε]ydΩ +
Asi σs [εsi ]yi
Mxk =
⎪
⎪
⎪
Ωr
⎪
i=1
⎪
⎪
nb
⎨
Myk = −
σc [ε]xdΩ −
Asi σs [εsi ]xi
(5.2.3)
⎪
Ωr
⎪
i=1
⎪
⎪
nb
⎪
⎪
⎪
⎪
σc [ε]dΩ +
Asi σs [εsi ]
⎪
⎩ Nk =
Ωr
i=1
dove nb `e il numero delle barre in acciaio, {xi , yi } ed Asi sono rispettivamente la posizione e
l’area della singola barra, σc [ε] `e la tensione nel calcestruzzo e σs [εsi ] la tensione nell’acciaio.
p
Scelti un numero m di possibili meccanismi plastici εk , rimane deﬁnito il poliedro di
p
m-facce, ciascuna delle quali di normale nk := εk e tangente al dominio elastico nel punto
σ k , che deﬁnisce il dominio nella forma (4.1.2):
fk [σ] ≡ {nTk σ − ck } ≤ 0 ,
ck = σ Tk nk
k = 1, · · · , m
(5.2.4)
dove σ := {Mx , My , N }T rappresenta la sollecitazione agente sulla sezione.
Ovviamente, l’accuratezza della rappresentazione dipende dal numero di facce del poliedro e dal modo in cui queste sono distribuite. Un compromesso eﬃcace tra accuratezza e
complessit`a computazionale pu`
o essere quello di scegliere, quali possibili cinematismi plastici
della sezione, le normali alle 26 facce del solido geometrico rappresentato in ﬁgura 5.18, le
Esempio di implementazione: i telai in c.a.
56
n5
n17
n25
n4
n9
n8
n3
n7
n1
e0
cx
cy
n6
εp
n1
n2
n3
n4
n5
n6
n7
n8
n9
n10
n11
n12
n13
χx
1
-1
0
0
0
0
ey
ey
-ey
-ey
ez
-ez
ez
χy
0
0
1
-1
0
0
ex
-ex
ex
-ex
0
0
0
ε0
0
0
0
0
1
-1
0
0
0
0
ex
ex
-ex
εp
n14
n15
n16
n17
n18
n19
n20
n21
n22
n23
n24
n25
n26
χx
-ez
0
0
0
0
eyz
eyz
-eyz
-eyz
-eyz
-eyz
eyz
eyz
χy
0
ez
ez
-ez
-ez
exz
exz
exz
exz
-exz
-exz
-exz
-exz
ε0
-ex
ey
-ey
ey
-ey
e
-e
e
-e
-e
e
e
-e
Figura 5.18: Cinematismi plastici
cui componenti sono corrette con i coeﬃcienti eαβ espressi in funzione delle sollecitazioni
limite relative alle direzioni principali di deformazione, nella forma:
ex :=| Mx2 − Mx1 | , ey :=| My4 − My3 | , ez :=| N6 − N4 |
ez
ez
ex + ey
, exz := ex , eyz := ey
e :=
2
ey
ex
(5.2.5)
Con l’introduzione dei coeﬃcienti eαβ si tiene conto del comportamento anisotropo della
sezione in termini di resistenze. Si osservi, infatti, che nel caso risultassero ex = ey = ez le
normali coinciderebbero con le normali ai piani del poliedro regolare di ﬁgura 5.18.
Si osservi che, operando in questo modo non si eﬀettua nessun controllo sulle deformazioni, ma ci si basa sulla sola resistenza a rottura della sezione. Un controllo sulle
deformazioni plastiche totali (veriﬁca a duttilit`a) `e realizzato, pi`
u coerentemente dal programma in sede di analisi, monitorando le sezioni e controllando le deformazioni plastiche
man mano accumulate, segnalando, se `e il caso, la formazione di deformazioni inaccettabili.
Il domino di resistenza di una sezione a “L”
Per illustrare i risultati ottenuti con l’approssimazione del dominio di interazione della
sezione, si far`a riferimento alla sezione a “L” di ﬁgura 5.19.
Il dominio di interazione e la sua approssimazione multifalda, per la sezione analizzata, sono rappresentati nelle ﬁgure 5.20–5.21. La rappresentazione riporta i valori delle
sollecitazioni limite adimensionalizzati e deﬁniti dalle relazioni:
mx =
Mxu
ex
,
my =
Myu
ey
,
n=
Nu
ez
(5.2.6)
in cui Mxu , Myu e Nu rappresentano le componenti della sollecitazione nella generica condizione di collasso della sezione. Si pu`
o osservare che la rappresentazione linearizzata si
presta bene a deﬁnire la superﬁcie del dominio limite per pressoﬂessione delle sezioni di
elementi in cemento armato.
Esempio di implementazione: i telai in c.a.
57
barra
1
2
3
4
5
6
7
8
9
10
x
3
16
27
44
57
57
57
44
27
27
y
3
3
3
3
3
15
27
27
27
42
barra
11
12
13
14
15
16
17
18
19
20
Calcestruzzo
Acciaio
ferri
x
27
27
27
15
3
3
3
3
3
3
y
58
72
87
87
87
72
58
42
27
15
C30/37
FeB44K
20 φ 20
Figura 5.19: Sezione a “L” analizzata
my
n
n
mx
my
mx
Figura 5.20: Proiezione del domino nei piani delle sollecitazioni principali
n
n
my
my
mx
mx
Figura 5.21: Rappresentazione tridimensionale del dominio di resistenza della sezione analizzata
Esempio di implementazione: i telai in c.a.
5.3
58
Risultati numerici
In questa sezione saranno presentati alcuni test numerici di analisi ad adattamento eseguiti
su telai in cemento armato. I primi due test si riferiscono a strutture semplici e riguardano
una trave continua appoggiata e un telaio piano. Con il terzo test, invece, sar`a analizzato
un telaio tridimensionale che corrisponde alla struttura di un ediﬁcio per civile abitazione.
I risultati numerici ottenuti sono rappresentati graﬁcamente sotto forma di domini di interazione che rappresentano le possibili condizioni di carico per cui si ha comportamento
elastico, adattamento e collasso della struttura.
5.3.1
La trave appoggiata
Il primo test riguarda l’analisi di shakedown eﬀettuta per la trave appoggiata di ﬁgura 5.22
sulla quale si considerano agenti i carichi p1 e p2 che variano nel tempo generando condizioni
di carico p[t] contenute nel dominio:
p1 = 2p2
0 ≤ α1 [t] ≤ 1
,
(5.3.1)
p[t] := α1 [t]p1 + α2 [t]p2 ,
0 ≤ α2 [t] ≤ 1
p2 = M0 /
in cui M0 `e il momento limite della sezione e `e rappresentata nella ﬁgura 5.22.
I risultati ottenuti sono sintetizzati nella ﬁgura 5.23, in cui `e stato rappresentato il dominio di shakedown, indicato con “S”, i cui punti corrispondono alle condizioni di carico
cui la struttura risponde elasticamente dopo un transitorio iniziale. Esso `e contenuto nel
dominio limite, il cui conﬁne identiﬁca le situazioni di collasso, ed `e esterno al domino elastico iniziale “E”, luogo delle condizioni di carico che non inducono nessuna deformazione
plastica se applicate sulla struttura nel suo stato naturale, indeformato e privo di sforzi.
Figura 5.22: Trave continua analizzata
Figura 5.23: Dominio di interazione
Esempio di implementazione: i telai in c.a.
5.3.2
59
Telaio piano
L’analisi `e stata eﬀettuata per il telaio piano di piccole dimensioni, la cui geometria e
caratteristiche meccaniche sono riportate in ﬁgura 5.24. Sono stati itilizzati un elemento
per ciascuna delle aste verticali e due elementi per l’asta orizzontale. Il dominio dei carici
`e deﬁnito come:
p[t] := α1 [t]p1 + α2 [t]p2
0 ≤ α1 [t] ≤ 1,
0 ≤ α2 [t] ≤ 2
(5.3.2)
con p1 = 60 KN e p2 = 45 KN .
Figura 5.24: Geometria del telaio piano
In ﬁgura 5.25 `e rappresentato il dominio elastico E, di shakedown S e di collasso C, al
variare del rapporto p2 /p1 . La rappresentazione `e adimensionalizzata rispetto al valore
p0 := 0.85bh2Rck 0.83/(γc) pari a 324.42KN , avendo assunto γc = 1.6, Rck la resistenza
caratteristica del calcestruzzo, b, h e i valori riportati nella ﬁgura 5.24.
p2
p1
Figura 5.25: Diagramma di interazione per il telaio piano
Esempio di implementazione: i telai in c.a.
60
Telaio spaziale
Il telaio spaziale riportato in ﬁgura 5.26 corrisponde alla struttura di un ediﬁcio per civile
abitazione. I pilastri e le travi sono state modellate con l’elemento asta 3d mentre ciascuno
dei campi di solaio `e stato modellato utilizzando due aste incernierate incrociate i cui nodi
coincidono con i vertici del campo di solaio.
Figura 5.26: Geometria del telaio spaziale
La struttura `e stata considerata sottoposta all’azione orizzontale p che simula l’azione
simica, le cui componenti {p1 , p2 } nel riferimento globale {x, y}, dipendono dalla direzione
β di incidenza del sima, nella forma:
p1 := p cosβ
,
p2 := p sinβ
(5.3.3)
Esempio di implementazione: i telai in c.a.
61
La forza p, la cui entit`
a variabile su ciascun livello `e riportata nella ﬁgura 5.26, `e applicata nel baricentro geometrico di ciascun livello. I carichi permanenti sono rappresentati
a dei carichi p3 e p4 , considerati
dall’azione p3 e i carichi accidentali dall’azione p4 . L’entit`
agenti uniformemente sulle travi, `e riportata nella tabella 5.3 per ciascun livello.
livello
1
2
p 3 [KN/m]
40.0
32.0
p4 [KN/m]
4.0
4.0
livello
3
4
p3 [KN/m]
23.0
18.0
p4 [KN/m]
4.0
4.0
Tabella 5.3: Carichi verticali uniformemente distribuiti sulle travi
Il dominio dei carichi `e stato deﬁnito nella forma:
p[t] :=
4
αi [t]pi
,
−1 ≤ α1 [t], α2 [t] ≤ 1,
1 ≤ α3 [t] ≤ 1,
0 ≤ α4 [t] ≤ 1
(5.3.4)
i=1
Anche in questo caso `e stato rappresentato, nella ﬁgura 5.27, il dominio dello stato elastico
E, di shakedown S e di collasso C, al variare del rapporto p2 /p1 in corrispondenza del
quarto livello. La rappresentazione `e riportata in forma adimensionale rispetto al valore
p0 := 0.85bh2 Rck 0.83/(γc) = 230.70KN .
Figura 5.27: Dominio di interazione del talaio spaziale
Capitolo 6
Esempio di implementazione: stati
piani di tensione e deformazione
Nel presente capitolo, l’analisi di shakedown discussa nei capitoli precedenti sar`
a estesa al
caso delle strutture bidimensionali piane caricate nel proprio piano e saranno analizzati i
casi di stato piano di tensione e stato piano di deformazione.
La particolarit`
a dell’analisi da eﬀettuare, molto simile all’analisi di tipo incrementale
utilizzata nell’analisi limite, ha richiesto la scelta di un elemento che fosse il pi`
u semplice possibile, adatto a discretizzazioni abbastanza ﬁtte, e privo di locking [42]. L’unica diﬀerenza
rispetto ad un classico elemento usato in analisi di plasticit`
a e che in questo caso l’elemento
deve essere tale da fornire una buona valutazione della soluzione elastica, necessaria per
¯
ottenere una stima accurata del valore di λ.
I risultati ottenuti mostreranno come l’elemento ﬁnito simplex proposto, che contiene
l’interpolazione delle sia tensioni sia degli spostamenti, si dimostra, rispetto ai requisiti
richiesti, un buon compromesso tra semplicit`
a e accuratezza.
6.1
L’elemento finito Simplex
Assumiamo un sistema di assi cartesiani (o, x, y) e indichiamo con u[x, y]T := {u[x, y], v[x, y]}
lo spostamento di un generico punto x := {x, y}T appartenente al continuo bidimensionale
di dominio B, con σ[x] := {σxx , σyy , σxy }T la tensione e con ε[x] = {εxx, εyy , εxy }T la deformazione. Indicando con b := {bx, by }T le forze di volume e con f := {fx , fy }T le forze
*
di superﬁcie agenti sulla parte Sf del contorno S = Sf Su del dominio B e assumendo
¯ in Su , la soluzione in forma discreta del problema elastico pu`
o essere
la condizione u = u
ottenuta mediante il classico principio di Hellingher–Reissner:
1 T −1
T
T
u bdV −
uT f dA = staz(u,σ ) (6.1.1)
ε σ − σ E σ dV −
ΠHR [σ, u] :=
2
B
B
Sf
La matrice elastica E per stati piani di tensione e stati piani di deformazione, rispettivamente, `e:
⎡
⎤
⎡
⎤
ν
1
0
1 ν 0
(1−ν)
E(1 − ν)
E ⎣
⎢ ν
1
0 ⎥
ν 1 0 ⎦ , E=
(6.1.2)
E=
⎣ (1−ν)
⎦
2
1−ν
(1 − 2ν)(1 + ν)
(1−2ν)
0 0 1−ν
0
0
2
2(1−ν)
62
Esempio di implementazione: stati piani di tensione e deformazione
63
Figura 6.1: L’elemento ﬁnito simplex
in cui `e stato indicato con E il modulo elastico di Young e con ν il coeﬃciente di Poisson.
Il dominio B del continuo bidimensionale `e suddiviso in elementi triangolari di area Ae
e spessore he ciascuno dei quali `e deﬁnito dai tre nodi posizionati ai vertici del triangolo.
All’interno dell’elemento si assume per gli spostamenti una interpolazionebilineare, mentre
nge
Agi , essendo
le tensioni si assumono costanti nell’area di inﬂuenza del nodo Ag := i=1
Agi := Aie /3 l’area nodale dell’elemento i–esimo e nge gli elementi collegati al nodo g
(v. ﬁg. 6.1).
Siano {1, 2, 3} i vertici nodali dell’elemento numerati in senso antiorario, indicheremo
con ui e σ i , rispettivamente, i valori che spostamento e tensione assumono al nodo i:
ui = u[xi , yi] ,
σ i = σ[xi , yi ]
i = 1, 2, 3
(6.1.3)
All’interno dell’elemento, invece, lo spostamento `e funzione dei valori che lo stesso assume
in corrispondenza dei nodi, attraverso la legge di interpolazione:
ue = N e [x, y]ue
(6.1.4)
in cui la matrice delle funzioni di forma N e [x, y] ed il vettore ue che raccoglie le componenti
nodali dello spostamento sono deﬁniti da:
N1 0 N2 0 N3 0
(6.1.5)
, ue := {u1 , u2 , u3 }T
N e [x, y] :=
0 N1 0 N2 0 N3
con
Ni [x, y] =
1
(ai + bi x + ciy) per
2Ae
i = 1, 2, 3
(6.1.6)
in cui i coeﬃcienti ai , bi, ci sono deﬁniti, attraverso una permutazione ciclica degli indici
{i, j, k} ≡ {1, 2, 3}, come
ai = xj yk − xk yj ,
bi = yj − yk ,
ci = xk − xj .
All’interno dell’elemento la deformazione compatibile εe risulta costante
⎧
⎫
⎡
b 0 b2 0 b3
⎨ u,x ⎬
1 ⎣ 1
v,y
0 c1 0 c2 0
= B e ue , B e [x, y] =
εe =
⎩
⎭
2Ae
u,y +v,x
c1 b 1 c2 b 2 c3
(6.1.7)
e deﬁnita come:
⎤
0
c3 ⎦
(6.1.8)
b3
Esempio di implementazione: stati piani di tensione e deformazione
64
La tensione σ g = σ[xg ], costante nell’area nodale Agi , `e espressa in termini di tensioni
naturali tg , che in questo caso facilitano l’applicazione dell’algoritmo rientro nel dominio
ammissibile e sono deﬁnte da:
⎧ ⎫
⎡
⎤
1 1 0
⎨ ts ⎬
1⎣
1 −1 0⎦ , tg := te
(6.1.9)
tg = Aσ g , A =
⎩ ⎭
2
tt
0 0 2
Indicando con q g il vettore che raccoglie gli spostamenti nodali ugi dei nodi che stanno ai
vertici degli elementi collegati al nodo g e con dg il vettore dei parametri delle interpolazioni
tensioni e spostamenti che interessano il nodo g
tg
T
, dg :=
(6.1.10)
q g := {ug1 , ug2 , · · ·ugnge }
qg
e introducendo la deformazione naturale γ e deﬁnita come:
γ e = A−1 εe = D e ue
,
De = A−1 B e
⎧
⎫
⎨ u,x +v,y ⎬
γ e := u,x −v,y
⎩
⎭
u,y +v,x
,
(6.1.11)
`e possibile trasformare la (6.1.1) in forma discreta attraverso l’equivalenza:
σ T εdV =
V
Ng
tTg
g=1
nge
A−1 εgi Agi hgi =
i=1
V
tTg γ g =
i=1
Ng
σ T E −1 σdV =
Ng
Ng
tTg Qg q g
(6.1.12)
i=1
Ng
σTg E −1 σ g Vg =
g=1
tTg F g tg
(6.1.13)
g=1
in cui si `e posto:
γ g := Qg qg
,
Vg :=
nge
Vgi
,
Vgi := Agi hgi
(6.1.14)
F g = A−1 E −1 A−T Vg
(6.1.15)
i=1
con:
Qg = {Dg1 Vg1 , Dg2 Vg2 , · · · Dgni Vgnge }
Utilizzando le quantit`
a appena introdotte, la matrice nodale del problema algebrico espresso
sia in termini di tensioni che di spostamenti diventa:
⎤
⎡
−F g De1 Ve1 De2 Ve2 . . . DeN VeN
⎥ ⎢ DT Ve1
0
0
...
0
⎥
⎢ e1
−F g Qg
⎥
⎢ D T Ve2
0
0
.
.
.
0
e2
=
(6.1.16)
Kg = ⎢
⎥
QTg 0
⎥
⎢
..
..
..
..
..
⎦
⎣
.
.
.
.
.
T
0
0
...
0
D eN VeN
con F g matrice diagonale. La soluzione del sistema lineare `e ottenuta condensando i
parametri tensione:
tg = F −1
g Qg q g = E g
ne i
i=1
Dgi ugi Vgi
,
E g = F −1
g
(6.1.17)
Esempio di implementazione: stati piani di tensione e deformazione
ottenendo, quindi, la matrice pseudo-compatibile
⎡
⎤
2
Vg1
DTg1 E g Dg1
Vg1 Vg2 D Tg1 E g Dg2 . . . Vg1 Vneg DTg1 E g Dneg
2
⎢ Vg1 Vg2 DTg2 E g D g1
Ve2
DTg2 E g Dg2
. . . Vg2 Vneg DTg2 E g Dneg ⎥
⎢
⎥
K cg ≡ ⎢
⎥
..
..
..
..
⎣
⎦
.
.
.
.
T
T
T
2
Vg1 Vneg Dneg E g Dg1 Vg2 Vneg D neg E g Dg2 . . . Vgneg Deni E g Deneg
65
(6.1.18)
Se indichiamo con d = {t, u}T il vettore che raccoglie i parametri di spostamento nel sistema
globale u = {u1 , u2 , · · · , uN }T e il vettore delle tensioni t = {t1 , t2 , · · · , tN }T del modello
discreto, attraverso un processo di assemblaggio esteso a tutti i punti di Gauss si ottiene:
A(K cg ) , s :=
A(sg ) =
A QTg tg
(6.1.19)
K e :=
g
g
g
o
in cui K e `e la matrice del sistema algebrico e s `e il vettore della risposta strutturale. Si pu`
osservare che la risposta strutturale pu`
o essere ottenuta facilmente mediante una procedura
standard di assemblaggio estesa a tutti gli elementi, nella forma:
s :=
,
se =
e
6.2
Ae he T De
tg
3
3
A(se )
(6.1.20)
g=1
La discretizzazione del dominio nel caso di stati piani di
tensione e deformazione
Si assume che il dominio elastico E sia deﬁnito dalla funzione di snervamento di Mises che,
per stati piani di tensione e stati piani di deformazione, rispettivamente, si scrive:
+
)
σy
σy
1 2
, f [t] = t2e + t2t − √
(6.2.1)
ts + t2e + t2t − √
f [t] =
3
3
3
in cui σy `e la tensione normale di snervamento.
Il dominio ammissibile di adattamento Es [λ] `e ottenuto mediante una appropriata linearizzazione a tratti della (6.2.1), usando un poliedro esterno regolare di m–facce, tale che il
k–esimo piano sia tangente al punto tk e la normale nk sia la normale alla funzione f [t] nel
punto tk . Ovviamente, nel caso di stati piani di deformazione il dominio elestico deﬁnito
dalla (6.2.1) `e rappresentato da un cerchio nel piano {te , tt}. La linearizzazione a tratti
dell’equazione del dominio (6.2.1) consente di descrivere la superﬁcie di snervamento come
un insieme di superﬁci piane nella forma f [t] := {f1 [t], f2 [t], · · · , fm[t]}T ed ottenere per il
dominio di shakedown la rappresentazione (4.1.3), nella forma:
, ¯bk := max {nTk te [t]} , ck = tTk nk k = 1, · · · , m
fsk [t, λ] = nTk t − ck + λ¯bk
te∈ Se
(6.2.2)
Ricordando la (1.2.2), che deﬁnisce il dominio delle sollecitazioni elastiche Se , si ha:
bk =
p
i=1
aki
con
aki :=
αmin
nTk tei se nTk tei < 0
i
max
αi nTk tei se nTk tei ≥ 0
(6.2.3)
Esempio di implementazione: stati piani di tensione e deformazione
66
Figura 6.2: Dominio elastico linearizzato (20 × 10) e dominio originale per stati piani di tensione
in cui bk `e valutato una sola volta all’inizio dell’analisi. Prima del processo iterativo `e
¯ oltre il quale esiste almeno un punto di Gauss p in cui il
valutato anche il moltiplicatore λ
dominio di shakedown Es [λ] `e vuoto.
Nei risultati numerici sar`
a mostrato l’inﬂuenza che la linearizzazione del dominio elastico
ha sul moltiplicatore di collasso e sul moliplicatore di shakedown.
Per gli stati piani di deformazione `e stata usata una discretizzazione regolare a 64 tratti
lineari. Per gli stati piani di tensione, invece, il dominio elastico `e stato discretizzato in
nr × nl piani, con nr la discretizzazione nel piano {te , tt} e con nl l’approssimazione nei
piani {te , ts } e {tt , ts } (v. ﬁg. 6.2).
6.3
Risultati numerici
In questa sezione saranno presentati alcuni test numerici che mostrano l’eﬃcienza del metodo proposto per la valutazione del moltiplicatore di shakedown di strutture bidimensionali
piane.
La sperimentazione numerica `e stata sviluppata in due parti. Nella prima l’attenzione
`e concentrata sugli aspetti che riguardano l’eﬃcienza della formulazione proposta rispetto
al tipo di elemento ﬁnito utilizzato, alla linearizzazione del dominio elastico e la rappresentazione della geometria della lastra. Nella seconda parte, sono stati ottenuti e confrontati
i risultati numerici di problemi classici di adattamento, presentando, inﬁne, una serie di
nuovi test.
C’`e da osservare che nel caso di plasticit`a alternata, il moltiplicatore di adattamento
¯ che pu`
coincide il moltiplicatore λ,
o essere valutato in modo veloce a valle dell’analisi
elastica e richiede un impegno computazionale simile alla valutazione del limite elastico.
¯ per`
o, non pu`
o essere stabilita a priori, per cui `e sempre necessario
La coincidenza λa = λ,
eﬀettuare l’intera analisi a shakedown, che `e la parte pi`
u onerosa. Ci`
o giustiﬁca l’impiego
di un elemento ﬁnito semplice capace di rendere pi`
u eﬃciente l’intera analisi.
¯
Inoltre, si capisce come sia importante una valutazione corretta del moltiplicatore λ.
` stato necessario selezionati test che non presentano singolarit`a nella soluzione elastica,
E
sostituendo gli spigoli con un contorno circolare, opportunamente discretizzato in elemento
Esempio di implementazione: stati piani di tensione e deformazione
67
rettilinei, in modo da evitare l’insorgere di tensioni spurie che avrebbero portato ad una
¯
sottostima del moltiplicatore λ.
L’analisi a shakedown richiede una soluzione elastica accurata capace di valutare correttamente le tensioni elastiche te . Inoltre, il processo incrementale richiede una discretizzazione ad elementi ﬁniti abbastanza ﬁtta per poter cogliere i meccanismi di deformazione
locali. In eﬀetti, il metodo proposto, consentirebbe di lavorare su mesh ﬁtte e ottenere la
soluzione elastica con un elemento ﬁnito accurato [43] ed utilizzare le tensioni cos`ı ottenute
su un elemento ﬁnito pi`
u semplice, impiegato all’interno del processo iterativo.
6.3.1
Lastra quadrata con foro centrale circolare
La lastra in esame `e un problema classico di analisi a shekedown (v. ﬁg. 6.3). Presenta
rapporto tra diametro del foro e lato pari a 0.2 ed `e soggetta ad un carico biassiale uniforme
p1 e p2 . L’analisi `e stata eﬀettuata nei casi di incremento proporzionale dei carichi p1 e p2
(analisi limite), variazione indipendente di p1 e p2 (0 ≤ α1 ≤ 1 e 0 ≤ α2 ≤ 1) e variazione
proporzionale dei carichi (0 ≤ α1 = α2 ≤ 1).
mesh1
mesh2
mesh3
mesh4
nodi
42
169
650
2496
elementi
72
288
1200
4896
Figura 6.3: Geometria e discretizzazioni della lastra quadrata con foro
Questo problema `e stato analizzato da diversi autori utilizzando tecniche diverse. Per
una visione completa dei risultati si pu`
o consultare il lavoro di Groß–Weedge [44], oppure
[45], [46].
L’analisi `e stata eﬀettuata per le discretizzazioni riportate nella ﬁgura 6.3. In particolare, indicando con il primo termine la discretizzazione del foro circolare e con il secondo la
Esempio di implementazione: stati piani di tensione e deformazione
68
discretizzazione del lato opposto, sono state usate le discretizzazioni (6 × 3) = 72 elementi,
(12 × 6) = 288 elementi, (24 × 12) = 1200 elementi e (48 × 24) = 4896 elementi per mesh
1, mesh 2, mesh 3, mesh 4 rispettivamente.
Nella tabella 6.1 sono riportati i risultati numerici ottenuti per il caso di stati piani di
tensione, utilizzando la mesh 4 e una discretizzazione del dominio elastico pari a 20 × 10.
I valori cos`ı ottenuti sono confrontati con quanto ottenuto da altri autori. Per ulteriori
approfondimenti riguardo le tecniche usate dagli altri autori `e utile consultare [44] e [46].
(p1 , p2 )
Belytschko [47]
Nguyen et al. [48]
Corradi et al. [2]
Genna [49]
Stein and Zhang [50]
Zhang et al. [51]
Gross-Weege [44]
Silveira [45]
Present (mesh 4)
(1,1)
0.704
0.767
0.893
0.882
0.894
0.902
λc
(1,0.5)
0.907
0.891
0.911
0.912
(1,0)
0.780
0.564
0.691
0.793
0.802
0.789
0.782
0.803
0.806
(1,1)
0.431
0.431
0.504
0.478
0.453
0.477
0.446
0.429
0.438
λa
(1,0.5)
0.501
0.514
0.579
0.566
0.539
0.549
0.524
0.500
0.508
(1,0)
0.571
0.557
0.654
0.653
0.624
0.647
0.614
0.594
0.604
Tabella 6.1: Lastra quadrata: valori di λc e λa a confronto
In ﬁgura 6.4 sono rappresentati il dominio elastico (E), il dominio di shakedown (S) e
di collasso (C) al variare del rapporto p1 /p2 . In questa caso, si noti che il moltiplicatore di
¯ Siamo in presenza di collasso per plasticit`
adattamento coincide sempre con il valore di λ.
a
alternata e l’accuratezza del moltiplicatore di shakedown dipende dall’accuratezza della
soluzione elastica.
Figura 6.4: Diagramma di interazione per la lastra forata nel caso di Figura 6.5: Funzione
trazioni p1 e p2 indipendenti e proporzionali
fn in corrispondenza dell’adattamento
Nel caso di variazione indipendente dei due carichi e per il rapporto p1 /p2 = 1, `e stata analiz-
Esempio di implementazione: stati piani di tensione e deformazione
69
zata l’inﬂuenza che si ha sul moltiplicatore di adattamento per eﬀetto della discretizzazione
ad elementi ﬁniti e per la diversa discretizzazione del domio elastico.
I risultati riportati nella tabella 6.2, confrontati con quelli ottenuti usando l’elemento
ﬁnito a deformazione costante (CT) e l’elemento a deformazione lineare (LT), evidenziano
una buona prestazione dell’elemento simplex soprattutto per mesh piuttosto rade. Si noti,
per`
o, che la soluzione elastica ottenuta con l’elemento simplex converge al valore esatto
pi`
u lentamente rispetto all’elemento LT e all’elemento CT. C’`e, per`
o, da considerare che a
parit`
a di mesh l’elemento LT presenta un numero maggiore di variabili.
Gli stessi risultati evidenziano che il moltiplicatore di collasso λc non `e praticamente
eﬀetto dell’errore di discretizzazione della geometria del foro, infatti, il valore ottenuto con
la mesh 1 presenta un errore del 3.6% rispetto a quello ottenuto con la mesh 4. Il valore
del moltiplicatore di shakedown λa, invece, che in questo caso `e strettamente legato alla
¯ risente maggiormente dell’errore
soluzione elastica poich`e concidente con il valore di λ,
introdotto con la discretizzazione del foro in elementi lineari.
dominio
(20 × 10)
mesh 1
mesh 2
mesh 3
mesh 4
CT
λc
0.8415 (0.8362)
0.8260 (0.8149)
0.8168 (0.8080)
0.8121 (0.8043)
¯
λa = λ
0.5216
0.4566
0.4394
0.4350
LT
λc
0.8179 (0.8092)
0.8119 (0.8042)
0.8083 (0.8017)
0.8065 (0.8007)
¯
λa = λ
0.4609
0.4396
0.4343
0.4329
Simplex
λc
0.8223 (0.8132)
0.8113 (0.8037)
0.8079 (0.8012)
0.8060 (0.8004)
¯
λa = λ
0.4536
0.4458
0.4428
0.4383
Tabella 6.2: Lastra quadrata: errore di discretizzazione nel caso di stati piani di tensione. In parentesi `e
riportato il valore ottenuto con la funzione di Mises esatta (λc `e calcolato per p1 /p2 = 0, λa per p1 /p2 = 1).
Riguardo all’errore introdotto con la discretizzazione del dominio elastico, si pu`
o osservare (v. tab. 6.3) che una discretizzazione (20 × 10) del dominio fornisce un risultato
abbastanza accurato, essendo, per ciascuna delle mesh, l’errore su λa inferiore al 2% rispetto valore ottenuto con una approssimazione del dominio (40 × 20). In generale, l’errore
contenuto nel moltiplicatore di collasso `e inferiore al 1.2% rispetto al valore che si ha senza
approssimazione del dominio elastico. In ﬁgura 6.5 `e riportata la mappa della funzione
di snervamento normalizzata fn = maxk {fsk } + 1, nello stato di adattamento, per il caso
p1 = 1 e p2 = 0.
mesh
mesh
mesh
mesh
1
2
3
4
(72 elem)
(288 elem)
(1200 elem)
(4896 elem)
Esatto E
λc
0.8132
0.8037
0.8012
0.8004
10 × 5
¯
λc
λa = λ
0.8640 0.4823
0.8589 0.4707
0.8579 0.4697
0.8578 0.4651
20 × 10
¯
λc
λa = λ
0.8223
0.4536
0.8113
0.4458
0.8079
0.4428
0.8060
0.4383
40 × 20
¯
λc
λa = λ
0.8160
0.4536
0.8060
0.4433
0.8030
0.4406
0.8020
0.4362
Tabella 6.3: Lastra quadrata: variazione di λc e λa con la discretizzazione del dominio elastico (λc `e
calcolato per p1 /p2 = 0, λa per p1 /p2 = 1).
La stessa analisi `e stata fatta nell’ipotesi di stati piani di deformazione, usando una
approssimazione regolare a 64 tratti del dominio elastico. I risultati sono riassunti nella
tabella 6.4. Si pu`
o notare, che in questo caso, l’eﬀetto del locking volumetrico distorce il
valore del moltiplicatore di collasso λc ottenuto con l’elemento a deformazione costante CT,
Esempio di implementazione: stati piani di tensione e deformazione
70
come ci si aspettava, mentre con l’elemento simplex proposto si ottengono buoni risultati
anche con mesh rade. La ﬁgura 6.6, che rapprentata l’evoluzione del processo di carico,
evidenzia il fenomeno del locking presente nell’elemento a deformazione costante. Il valore
¯ non risente del problema del locking.
di λa, poich`e in questo caso coincide con il valore di λ,
domain
64 elements
mesh 1
mesh 2
mesh 3
mesh 4
CT
λ∗c
1.0903 (1.0894)
1.0854 (1.0846)
1.0740 (1.0737)
1.0424 (1.0425)
¯
λa = λ
0.6082
0.5298
0.5073
0.5010
LT
λc
0.9498 (0.9491)
0.9357 (0.9351)
0.9282 (0.9275)
0.9267 (0.9257)
¯
λa = λ
0.5275
0.5037
0.4987
0.4973
Simplex
λc
0.9402 (0.9396)
0.9212 (0.9207)
0.9224 (0.9219)
0.9234 (0.9230)
¯
λa = λ
0.5282
0.5093
0.5065
0.5019
Tabella 6.4: Lastra quadrata: errore di discretizzazione per l’elemento simplex nel caso di stati piani di
deformazione. Il valore tra parentesi `e ottenuto utilizzando la funzione esatta di Mises. (*) I valori con
l’elemento CT sono ottenuti per uno spostamento massimo va = 0.05.
Figura 6.6: Lastra quadrata: evoluzione del carico e meccanismo di collasso per il caso di stati piani
deformazione
6.3.2
Ulteriori test con concentrazione di tensioni
Nel seguito sono presentati una serie di test eseguiti su lastre quadrate che presentano
fori di varia geometria e diverse condizioni di carico. La geometria `e la discretizzazione
in elementi ﬁniti, per ciascuno dei test, sono rappresentate nella ﬁgura 6.7. Tutti i test
sono stati analizzati sia nell’ipotesi di stati piani di tensione, sia nell’ipotesi di stati piani
di deformazione.
Per il caso di stati piani di tensione e per ogni test `e stata usata la discretizzazione 20×10
del dominio elastico, poich`e nel test precedente si `e visto che fornisce buoni risultati. Per
stati piani di deformazione, invece, `e stata utilizzata una rappresentazione del dominio
regolare a 64 tratti. La geometria della lastra `e stata descritta sostituendo agli spigoli
Esempio di implementazione: stati piani di tensione e deformazione
test 1.
test 2.
test 3.
test 4.
71
Figura 6.7: Geometria e discretizzazione in elementi ﬁniti
interni un arco di cerchio approssimato con 32 elementi lineari nel caso del test 1, test 2 e
test 4, e con 12 elementi lineari per il test 3.
Come per il caso precedente, per ciascun test, `e stato analizzato sia il caso di variazione indipendente di p1 e p2 (0 ≤ α1 ≤ 1 and 0 ≤ α2 ≤ 1), sia il caso di variazione
proporzionale (0 ≤ α1 = α2 ≤ 1). Nella tabella 6.7 sono riportati i valori del moltiplicatore
di shakedown e del moltiplicatore di collasso ottenuti per stati piani di tensione e per stati
piani di deformazione. In ciascuno dei casi, per la presenza di concentrazioni di tensioni, il
¯ L’accuratezza dei risultati, quindi,
moltiplicatore di shakedown coincide con il valore di λ.
dipende dall’accuratezza della soluzione elastica.
6.3.3
Alcuni test con collasso incrementale
Nei casi analizzati in precedenza `e stata veriﬁcata la coicidenza del moltiplicatore di shake¯ Ci`
down con il valore di λ.
o si veriﬁca per la presenza di un escursione alta e locale delle
Esempio di implementazione: stati piani di tensione e deformazione
test
test 1
test 2
test 3
test 4
nodi
1522
972
821
1522
elementi
2877
1824
1510
2877
E
σy
=
=
200000
10
72
ν
t
=
=
0.3
1
Tabella 6.6: Propriet`a del materiale
(t=spessore della lastra)
Tabella 6.5: Discretizzatione
(p1 , p2 )
test 1
test 2
test 3
test 4
λc
(1,1)
5.795
5.644
6.941
5.770
Stati piani di tensione
λa,indipendente
(1,0)
(1,1)
(1,0)
5.046
1.631 1.786
5.526
2.931 3.669
5.979
1.821 2.684
5.820
1.826 1.986
(∗)
λa
(1,1)
1.974
3.550
4.625
2.176
λc
(1,1)
10.016
7.925
11.671
9.970
Stati piani di deformazione
λa,indipendente
(1,0)
(1,1)
(1,0)
5.794
1.882 2.061
6.121
3.348 4.236
6.851
2.088 3.094
6.515
2.107 2.292
(∗)
λa
(1,1)
2.279
4.084
5.299
2.513
Tabella 6.7: Valori dei moltiplicatori λc e λa per stati piani di tensione e stati piani di deformazione. (*)
Carichi combinati
tensioni, fenomeno tipico della soluzione elastica dei problemi analizzati, che porta al collasso per plasticit`
a alternata. Comunque, come dimostrato dai test che seguono, si pu`
o
¯ e λc , per cui il
veriﬁcare che il moltiplicatore di adattamento `e diverso dal valore di λ
collasso si ha per incremento di deformazione plastica (ratcheting).
Analisi di un telaio semplice
Il primo test riguarda il telaio semplice riportato in ﬁgura 6.8 analizzato per due diﬀerenti
condizioni di vincolo indicati con caso a (spostamento verticale del lato e-f impedito) e caso
b (spostamento verticale ed orizzontale impediti sul lato e-f ). In entrambi i casi l’analisi `e
stata condotta per il dominio dei carichi riportato nella tabella 6.8.
La soluzione elastica analitica per il test in esame presenterebbe tensioni inﬁnite in
corrispondenza della discontinuit`
a della normale al bordo. Tale circostanza, frequente per
le strutture bidimensionali piane, pu`
o essere eliminata solo con opportuni raccordi degli
spigoli. Il fenomeno non `e presente nella soluzione ad elementi ﬁniti che presenta un valore
¯ = 0, anche per discretizzazioni del dominio abbastanza ﬁtte, come quella analizzata.
di λ
Con il presente test si intende mostrare come sia possibile determinare il meccanismo di
adattamento anche quando questo `e diverso dal meccanismo di collasso.
p1
p2
=
=
3.0
1.0
0.4
0.4
≤
≤
α1 [t]
α2 [t]
≤
≤
1
1
Tabella 6.8: Dominio dei carichi per il telaio semplice
Per entrambi i casi analizzati, il moltiplicatore di shakedown ottenuto `e diverso dal
¯ e anche da λc , i quali, per una serie di combinazini di carico, sono riportati nella
valore di λ
tabella 6.10. Inﬁne, in ﬁgura 6.9 `e rappresentata la curva del processo incrementale ed il
meccanismo di ratcheting
Esempio di implementazione: stati piani di tensione e deformazione
73
Figura 6.8: Telaio semplice
Figura 6.9: Processo evolutivo del carico e meccanismo di ratcheting per caso a e caso b
Analisi al collasso
caso a.
λc [q2 ] = 2.831
λc [q3 ] = 2.975
λc [q4 ] = 2.645
caso b.
λc [q2 ] = 4.207
λc [q3 ] = 7.804
λc [q4 ] = 3.949
Analisi di shakedown
caso
λe =
λa =
¯ =
λ
a.
1.203
2.473
2.940
caso
λe =
λa =
¯ =
λ
b.
1.355
3.925
4.518
Figura 6.10: Dominio dei carichi e risultati dell’analisi per il telaio semplice
Esempio di implementazione: stati piani di tensione e deformazione
74
Trave continua simmetrica
Anche in questo ultimo test, la cui geometria, propriet`
a dei materiali e discretizzazione
sono riportati nella ﬁgura 6.11, il collasso avviene per ratcheting. Il dominio dei carichi `e
riportato nella tabella di ﬁgura 6.12.
¯ = 5.304, mentre i
L’analisi ha fornito il moltiplicatore di shakedown λa = 3.244 e λ
valori di λc , per una serie di combinazioni di carico, sono riportati nella tabella 6.12. In
ﬁgura 6.12 `e rappresentata la curva di carico in termini di moltiplicatore λa e spostamento
verticale del punto a.
Figura 6.11: Geometria e discretizzazione ad elementi ﬁniti della trave continua
Dominio dei carichi
p1 = 2
0.6 ≤ α1 ≤ 1.0
0.0 ≤ α2 ≤ 1.0
p2 = 1
Moltiplicatori di collasso
(p1 , p2 )
λc
(2.0, 0.0)
3.280
(0.0, 1.0)
8.718
(1.2, 1.0)
5.467
(2.0, 1.0)
3.280
Figura 6.12: Meccanismo di ratcheting, dominio dei carichi e fattore di sicurezza a shakedown
Considerazioni conlcusive
Nel presente lavoro di tesi `e stato presentato il metodo numerico appropriato per l’analisi
di shakedown proposto da Casciaro–Garcea [10]. Il metodo, di tipo iterativo–incrementale,
`e capace di determinare la soluzione di adattamento in contesti FEM del tutto generici e si
presenta molto pi`
u eﬃciente di altri metodi proposti in letteratura ([52], [53]), specie per
problemi a grandi dimensioni.
Il metodo, le cui caratteristiche di convergenza sono ampiamente discusse in [10], dal
punto di vista implementativo, presenta lievi diﬀerenze rispetto all’algoritmo standard tipo
path–following usato attualmente in analisi elasto–plastiche per valutare il percorso di equilibrio di una struttura. Questa caratteristica lo rende facilmente perfezionabile in contesti
FEM del tutto generici, come `e stato dimostrato dalle applicazioni fatte per l’analisi di telai
spaziali in cemento armato e per l’analisi di strutture bidimensionali piane.
L’analisi di adattamento di telai spaziali `e stata condotta su un modello numerico generato modellando gli elementi travi e pilastri della struttura in cemento armato mediante
elementi ﬁniti monodimensionali. L’elemento ﬁnito utilizzato `e stato derivato dalla teoria
del De Saint Ven`
ant `e tiene conto dell’iterazione taglio–torsione che si genera nell’elemento spaziale per eﬀetto dell’eccentricit`a dell’azione tagliante rispetto al centro di taglio. Il
comportamento elastoplastico della sezione dell’elemento `e stato deﬁnito mediante una rappresentazione lineare a falde del dominio limite per pressoﬂessione deviata, ottenuta in modo
automatico a partire dalla geometria della sezione e dalle armature presenti. La scelta di
utilizzare 26 piani per approssimare il dominio elastico della sezione `e risultata un buon
compromesso tra errore introdotto e impegno computazionale richiesto.
L’analisi delle strutture bidimensionali `e stata condotta utilizzando l’elemento ﬁnito
simplex. Tale elemento, di tipo misto, `e molto semplice e si presta bene ad essere usato in
contesti di plasticit`
a dove `e necessaria una discretizzazione abbastanza ﬁtta del dominio.
Inoltre, la necessit`
a di avere anche una soluzione elastica accurata, rende l’elemento ﬁnito
simplex da preferire rispetto a elementi standard triangolari a tre e a sei nodi.
I risultati numerici ottenuti hanno dimostrato come il metodo si presta bene a determinare il moltiplicatore di adattamento di strutture generiche a dimensione reale. Inoltre,
hanno confermato che la valutazione della sicurezza della struttura basata sulla sola analisi
limite, nel caso di strutture su cui agiscono carichi che variano nel tempo, fornisce fattori di sicurezza superiori a quelli previsti dall’adattamento in campo elastico. Ci`o rende
necessaria, per le strutture reali, la veriﬁca della sicurezza basata sul moltiplicatore di adattamento elastico. In questo contesto, il metodo esposto potrebbe essere considerato uno
strumento appropriato per la veriﬁca di strutture in ambito tecnico–professionale.
75
Bibliografia
[1] G. Ceradini, C. Gavarini. Limit analysis and linear programming. Giornale del Genio Civile, Gen/Feb
(1965).
[2] L. Corradi, A.Zavelani. A linear programming approach to shakedown analysis of structures. Comp.
Meth. Appl. Mech. Engng, 3, 37-53 (1974).
[3] G. Maier. A Matrix Structural Theory of Piecewise Linear Elastoplasticity. Meccanica, 5, 54-66 (1970).
[4] M. Capurso. A displacement bounding principle in shakedown of structures subjected to cyclic loads.
Int. J. Solids Structures, 10, 77-92 (1974).
[5] R. Casciaro, L. Cascini. A mixed formulation and mixed ﬁnite elements for limit analysis. Int. J.
Numer. Meth. Engrg, 18, 211-243 (1982).
[6] R. Casciaro, M. Mancuso. Un approccio numerico al problema dell’elastoplasticit`
a incrementale.
Omaggio a Giulio Ceradini, Dipartimento di Ingegneria Strutturale e Geotecnica, Univ. di Roma
La Sapienza, 1988
[7] R. Casciaro, L. Cascini. Limit analysis by incremental iterative procedure. Proceeding of the IUTAM
Conference on Defrmation and failure of Granualr Materials, Delft, 1982.
[8] E. Melan. Zur Plastizit¨
at des ra¨
umlichen continuum. Ing. Arch, 9, 116-126, (1938)
[9] W.T. Koiter. General theorems for elastic–plastic solids. Progress in solids mechanics, ed I. N. Sneddon
and R Hill, 165-221, North–Holland, Amesterdam, 1960.
[10] R. Casciaro, G. Garcea. An iterative method for shakedown analysis. Computer Methods in Applied
Mechanics and Engineering, 191, 5761-5792 (2002)
[11] A.S. Petrolo, R. Casciaro. 3D beam element based on De Saint Venant’s rod theory. Computers and
Structures, 82, 2471-2481 (2004)
[12] M. Malena, A.S. Petrolo, R. Casciaro. Caratterizzazione di domini di interazione e suo uso nell’analisi
non lineare di telai 3D in c. a. Report 29, LabMec, Dipartimento di Strutture, UNICAL, 2003.
[13] A. K¨
onig. Shakedown of elastic–plastic structures. Fundamental Studies in Engineering 7, Elsevier,
1987.
[14] R.T. Rockafellar. Convex analysis. Princeton univ. Press, Princeton, 1972
[15] G. Ceradini. Sull’adattamento dei corpi elasto plastici soggetti ad azioni dinamiche. Gionale del Genio
Civile, 106, n.4/5, 239-250.
[16] G. Ceradini. C. Gavarini. Applicazione della programmazione lineare in problemi di adattamento
plastico dinamico. Giornale del Genio Civile, 106, n.8, 239-250 (1969)
[17] C. Polizzotto, G. Borino, S. Caddemi, P. Fuschi. Theorems of restricted dynamic shakedown. Int. J.
Mech. Sci., vol. 35 n◦ 9, 787-801 (1993).
76
77
[18] G. Borino, C. Polizzotto. Dinamic shakedown of structures under repeated seismic loads. J. Engng.
Mechanics, 1306-1314 (1995).
[19] H. Bleich. Ueber die Bemessung statisch umbestimmter Stanhltragwerke unter Ber¨
uchsichtigung des
elastisch–plastischen Verhaltendes Baustoﬀes. Bauingenieur (1932).
[20] W.T. Koiter. A new general theorem on shakedown of elastic–plastic structures. Proc. Koninkl. Ned.
Akad. Wet., B 59, 24–34 (1956)
[21] M. Capurso. Variational theorems and related methods of solution for the rate problem in elastoplastic
solids with singular yield surface. J. Struct. Mech., 3:1 (1974).
[22] Simo, Kennedy, Govindjee. Non–smooth multisurface plasticity and viscoplasticity. Loading/Unloading
conditions and numerical algoritms. Int. J. Num. Meth.Engrg., 26, 2161-2185 (1998).
[23] L. Corradi. On compatible ﬁnite element models for elastic plastic analysis. Meccanica, 13, 133-150
(1978)
[24] L.Corradi. On stress computation in displacement ﬁnite element models. Comp. Meth. Appl. Mech.
Engng, 54, 325-339 (1986)
[25] E. Riks. An incremental approach to the solution of snapping and buckling problems. Int. J. Solids
Structures, 15, 529-551 (1979).
[26] E. Riks. Buckling analysis of elastic structures; a computational approach. Advances in Applied
Mechanics , E. van der Giessen, T.Y Wu eds., Academic Press, 1997
[27] G.Garcea, G.A.Trunﬁo, R. Casciaro. Mixed formulation and locking in path following nonlinear
analysis. Comput. Meth. Appl. Mech.Engrg., 165 1-4, 247-272 (1998).
[28] G. Maier. Upper bounds on dynamic deformations of elastoplastic continua. Meccanica , 9, 30-42
(1974).
[29] M. Capurso, L. Corradi and G. Maier. Bounds on deﬂections and displacemnts in shakedown theory.
Mat´eriaux et Structures sous Chargment Cyclique, 231-244, Ass. Amicale des Ingenieurs Ancient El´eves
de E.N.P.C., Paris, 1979.
[30] R. Fletcher. Practical methods of optimization, J. Wiley & S., 1969.
[31] E. Gill, W. Murray. Numerically stable methods for quadratic programming. Mathematical
Programming, 14, 349-372 (1978).
[32] D. Goldfarb, A. Idnani. A numerically stable dual method for solving strictly convex quadratic
programs. Mathematical Programming, 27, 1-33 (1983).
[33] G. Armentano. Il principio di Haar-Karman e la programmazione matematica in plasticit`
a multifalda.
PhD thesis Computational Mechanics, UNICAL, Cosenza, 2004.
[34] EAH. Love. A treatise on the mathematical theory of elasticity. New York, Dover, 1944.
[35] G. Fichera. Remarks on Saint Ven`
ant principle. I.N. Vekua 70th anniversary volume, Moskow, 1977.
[36] JH. Argyris et al. Finite Element Method - the Natural Approach. Comp. Meth. Appl. Mech. Engng.
17/18, 1-106 (1979).
[37] OC. Zienkiewic. The ﬁnite element method in engineering science. New York, McGraw-Hill, 1984.
[38] U. Schramm, V. Rubenchik, W. Pilkey. Beam stiﬀness matrix based on the elasticity equations.
International Journal for Numerical Methods in Engineering, 40, 211-232 (1997).
78
[39] E.J. Sapountzakis. Solution of non-uniform torsion of bars by an integral equation method. Computers
and Structures , 77, 659-667 (2000).
[40] F. Gruttmann, W. Wagner, R. Sauer. Zur Berechnung von Woolbfunckion und Torsionskennwerten
beliebiger Stabquerschnitte mit der Methode der Finiten Elemente. Bauingenieur, 73(3), 138-143
(1998).
[41] S.I. Chou. Determination of centers of ﬂexure using the boundary element method. Engineering
Analysis with Boundary Elements, 12, 321-324 (1993).
[42] T. Belytschko, W.K.. Liu, B. Moran. Nonlinear ﬁnite elements for continua and structures. John Wiley
& Sons, LTD, 2000.
[43] C.Felippa. A study of optimal membrane triangles with drilling freedoms. Comput. Meth. Appl. Mech.
Engrg., 193, 2125-2168 (2003).
[44] J. Groß-Wedge. On the numerical assessment of the safety factor of elastic-plastic structures under
variable loading. Int. J. Mech. Sci., 39(4), 417-433 (1997).
[45] N. Zouanin, L. Borges, J.L. Silveira. An algorithm for shakedown analysis with nonlinear yeld function.
Comp. Meth. Appl. Mech. Engng., 191, 2463-2481 (2002).
[46] X. Zhang, Y. Liu, Z. Cen. Boundary element methods for lower bound limit and shakedown analysis.
Eng. Anal.Boundary Methods, 28, 905-917 (2004).
[47] T. Belytschko. Plane stress shakedown analysis by ﬁnite element. Int. J. Mech. Sci., 14, 619-625
(1972).
[48] H. Nguyen-Dang, L. Palgen. Shakedown analysis by displacement method and equilibrium ﬁnite
element. Proc. SMIRT 5, Berlin paper L3/3, 1979.
[49] F. Genna. A nonlinear inequality, ﬁnite element approach to the direct computation of shakedown load
safety factor. Int. J. Mech. Sci., 30, 769-789 (1988).
[50] E. Stein, G. Zhang. Shakedown with nonlinear strain–hardening including structural computation
using ﬁnite element method. Int. J.Plasticity, 8, 1-31 (1992).
[51] Yaung-GAo Zhang. An iteration algorithm for kinematic shakedown analysis. Comp. Meth. Appl.
Mech. Eng., 127, 217-226 (1995).
[52] G. Maier, V. Carvelli, G. Cocchetti. On direct methods for shakedown and limit analysis. Plenary
Lecture, 4th Euromech Solid Mechanics Conference, Metz, June 2000, European Journal of Mechanics
A/Solids, 19, Special Issue, S79-S100, 2000
[53] A.R.S. Ponter, K.F. Karter. Shakedown state simulation tecniques based on linear elastic solutions.
Comp. Methods Appl. Mech. Engng., 140, 259-279 (1997).

Download Report