Claudio Borri
Luca Salvatori
MECCANICA COMPUTAZIONALE
Capitolo 6
Introduzione all’analisi non-lineare
Rev. 12/06/2006
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 1/25
Argomenti trattati nel capitolo 6
Non-linearità fisica e geometrica
Metodi di analisi
Procedura incrementale
Metodi iterativi
Criteri di convergenza
Esempio di non-linearità geometrica: biella in
formulazione lagrangiana totale
Esempio di non-linearità fisica: biella con legame
elastico non-lineare
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 2/25
1
Introduzione
P
Struttura
V
INPUT
Sistema
fisico
OUTPUT
(carichi esterni)
V = V ( P,"t" )
(spostamenti,
deformazioni,
sforzi, …)
sistema di equazioni
(in generale non-lineare)
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 3/25
Caso lineare
V = V ( P,"t" ) = K -1eP
sistema di equazioni lineari
☺ Vantaggi:
soluzione con metodi standard
computazionalmente poco oneroso
vale il “principio” di sovrapposizione degli effetti
(combinazioni di carico, sovrapposizione modale, …)
L’approssimazione lineare non è sempre applicabile!
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 4/25
2
Possibili fonti di non-linearità
Gruppi di equazioni fondamentali
STATICHE
(equilibrio)
COSTITUTIVE
(legame)
CINEMATICHE
(congruenza)
non-linearità fisica
non-linearità
geometrica
(plasticità dell’acciaio,
fessurazione del cls,…)
σ ij
(equilibrio della
struttura deformata:
= σ ij (ε," t ")
l’eq. in conf. indef. è
un controsenso
fisico!!)
1 ⎛ ∂ui ∂u j ∂uk ∂uk ⎞
εij = ⎜ +
+
⎟
2 ⎜⎝ ∂x j ∂xi ∂xi ∂x j ⎟⎠
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 5/25
Condizione di equilibrio di un sistema discreto
e linearizzazione
Immaginiamo di poter scrivere l’equilibrio in forma vettoriale come:
FFint
(V) = P
int(V) = P
condizione di equilibrio
Fint vettore delle forze interne (dipendente dalla struttura)
V spostamenti nodali (incogniti)
P carichi esterni (dati)
Caso lineare: Fint(V)=Klin V
Caso non-lineare: linearizzazione attorno a V0 in equilibrio (Fint(V0)=P0 ):
Fint ( V + Δ V ) = Fint ( V0 ) +
∂ Fint
∂V
Δ V + ... = P0 + Δ P = P
V0
K T ⋅ ΔV P − Fint
K T ( V0 ) =
∂Fint
∂V
equilibrio linearizzato
matrice di rigidezza tangente
V0
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 6/25
3
Obbiettivo dell’analisi non-lineare e
la procedura incrementale-iterativa
Obbiettivo: tracciare la risposta carico-spostamento della struttura.
Si assume il carico dipendente da un solo parametro scalare λ
(moltiplicatore di carico):
P = λ Pref
Come il FEM discretizza la struttura nello spazio, la procedura
incrementale effettua la discretizzazione del percorso caricospostamento: dato uno stato in equilibrio definito dalla coppia (Vi,λi)
la procedura incrementale stabilisce come si debba ricercare il
successivo stato (Vi+1,λi+1).
Poiché le equazioni di equilibrio sono non-lineari, ogni coppia in
equilibrio è ricercata con l’ausilio di un metodo iterativo.
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 7/25
“Ingredienti”
Procedura incrementale (criterio con cui vengono ricercate le coppie (V,λ)):
• Controllo di carico
• Controllo di spostamento
• “A lunghezza d’arco”
Metodo iterativo (metodo con cui si risolvono le equazioni di equilibrio):
• Rigidezza tangente aggiornata (Newton-Raphson o “N-R”)
• Rigidezza tangente non aggiornata (metodi “modificati”)
• Rigidezza elastica (o “rigidezza costante” o “rigidezza lineare”)
• “Rigidezza” secante (quasi-Newton)
Criterio di convergenza (grandezza con cui si misura l’errore e che determina la
fine delle iterazioni):
• Norma dei residui
• Norma degli spostamenti
• Norma energetica
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 8/25
4
Esempio: contr. di carico + N-R + err. ass. residui
Si parte da una posizione di equilibrio nota
(λ0=0,V0=0)
(e.g. supponendo che la struttura scarica
sia in configurazione indeformata).
λi = λi-1+Δλ
All’inizio di ogni passo incrementale si
V*i,0 = Vi-1
assume come prima ipotesi per gli
spostamenti il valore al passo precedente.
Fint = Fint(V*i,k-1)
Con esso si valutano le forze interne e la
KT = KT(V*i,k-1)
rigidezza tangente (metodo di Newton).
Si calcola il vettore dei residui.
Ri,k = λiPref-Fint
Si calcola la correzione degli spostamenti
ciclo
e il relativo valore aggiornato.
ΔV*i,k = KT-1 Ri,k
iterativo
Si verifica la bontà dell’approssimazione
V*i,k = V*i,k-1+ΔV*i,k
(k=1,2,…)
(in questo esempio con l’errore assoluto
misurato sui residui): se non c’è
||Ri,k||<δtoll
convergenza si continua ad iterare,
no
altrimenti si accetta il valore V*i,k come in
sì
equilibrio con il moltiplicatore λi.
Vi = V*i
La procedura ha termine quando si
raggiunge un valore massimo prestabilito
λi ≥ λmax
no
di λ. Il risultato sono le coppie (Vi,λi).
sì
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
ciclo
incrementale
(i=1,2,…)
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 9/25
FINE
Procedura incrementale: controllo di carico
λPref
punto limite (snap-through)
(Vi,λiPref)
(V2,λ2Pref)
(V1,λ1Pref)
V
☺ corrisponde ad un modo “naturale” di applicare il carico per passi
impossibile seguire il percorso di equilibrio oltre il punti limite
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 10/25
5
Procedura incrementale: controllo di spostamento
P
punto limite (snap-through)
(λiVref,Pi)
(λ2Vref, P2)
punto di ritorno
(snap-back)
(λ1Vref,P1)
λVref
Il moltiplicatore si applica agli spostamenti e si ricercano i carichi corrispondenti:
tutti gli spostamenti assegnati devono essere fra loro proporzionali, ciò accade
solo in pochi modelli fisici (e.g. carichi concentrati in prove sperimentali).
☺ consente di superare i punti limite
ma non i punti di ritorno
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 11/25
Procedura incrementale: “a lunghezza d’arco”
λPref
punto limite (snap-through)
A = Δλ 2 + ( cΔV )
(V2,λ2Pref)
(Vi,λiPref)
2
punto di ritorno
(snap-back)
(V1,λ1Pref)
V
non si ha il controllo né sul carico né sugli spostamenti
☺ è possibile seguire ogni percorso di equilibrio
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 12/25
6
Metodi iterativi: Newton-Raphson
λPref
ΚT(V*i,2)
ΚT(V*i,1)
(Vi+1,λi+1Pref)
λi+1Pref
ΚT(Vi)
i,2,Fint
(V*
Ri,1
(V*
i,2))
(V*i,1,Fint(V*i,1))
λiPref
☺
☺ Convergenza
Convergenza
quadratica
quadratica(poche
(poche
iterazioni).
iterazioni).
KKT viene
calcolata ed
T viene calcolata ed
invertita
invertitaad
adogni
ogni
passo
(onere
passo (onere
computazionale).
computazionale).
(Vi,λi)
V
ΔV*i,1 ΔV*i,2 ΔV*i,3
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Ri,2
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 13/25
Metodi iterativi: Newton-Raphson modificato
λPref
ΚT(Vi) ΚT(Vi)
λi+1Pref
ΚT(Vi)
(V*i,2,Fint(V*i,2))
*
i,1,Fint(V i,1))
(V*
λiPref
(Vi,λi)
ΔV*i,1 ΔV*i,2 ΔV*i,3
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
(Vi+1,λi+1Pref)
Ri,1
Ri,2
☺
☺ KKTTviene
vienecalcolata
calcolataed
ed
invertita
una
invertita unasola
solavolta:
volta:
ciascuna
ciascunaiterazione
iterazioneèèpiù
più
rapida.
rapida.
La
Laconvergenza
convergenzadiviene
diviene
lineare
lineareeedunque
dunquesono
sono
necessarie
necessariepiù
piùiterazioni.
iterazioni.
V
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 14/25
7
Metodi iterativi: quasi-Newton (della secante)
λPref
ΚS(V*i,1)
ΚS(V*i,2)
λi+1Pref
ΚT(Vi)
i,2,Fint
(V*
(V*
i,1,Fint
λiPref
(V*
(Vi+1,λi+1Pref)
i,2))
(V*
i,1))
(Vi,λi)
Ri,1
•• KKT viene
calcolata ed invertita solo
T viene calcolata ed invertita solo
alalprimo
passo,
primo passo,aiaipassi
passisuccessivi
successivi
sisiutilizza
la
matrice
utilizza la matricedi
dirigidezza
rigidezza
secante
secante(più
(piùrapida
rapidada
dacalcolare).
calcolare).
•• La
Lavelocità
velocitàdi
diconvergenza
convergenzaèè
intermedia
fra
intermedia fraililmetodo
metododi
diN-R
N-R
standard
e
quello
modificato
standard e quello modificato
V
ΔV*i,1 ΔV*i,2 ΔV*i,3
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Ri,2
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 15/25
Metodi iterativi: rigidezza iniziale
Il metodo della rigidezza iniziale (o elastica o costante) utilizza
sempre la matrice di rigidezza elastica.
Il metodo richiede in genere molte iterazioni e presenta la
convergenza più lenta.
In alcuni casi tuttavia rende l’analisi più stabile.
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 16/25
8
Metodi iterativi applicati
alla procedura “a lunghezza d’arco”
iterazioni lungo arco di
circonferenza
iterazioni lungo normale fissa
Esistono diverse varianti della procedura
“a lunghezza d’arco” che, assieme al tipo di
iterazione con cui vengono implementati,
prendono nomi specifici (e.g. il metodo di
Riks-Wemper-Wessels è una procedura a
controllo di lunghezza d’arco che proietta
lungo la normale, e itera alla Newton
aggiornando le matrici).
iterazioni lungo normale
aggiornata (R-W-W)
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 17/25
Criteri di convergenza
La procedura iterativa termina quando:
errore relativo o assoluto
dell’iterazione corrente
(calcolato dal programma)
δk =
ek
< δ toll
eref
precisione ricercata
(fissata dall’utente)
Le possibili varianti sono:
a) nella scelta della grandezza da normare:
1) norma dei residui (Ri,k)
2) norma delle correzioni degli spostamenti (ΔV*i,k)
3) norma energetica (prodotto di Ri,k per ΔV*i,k)
b) nella scelta del tipo di norma (euclidea, max-abs, sum-abs, …).
c) nella scelta della quantità di riferimento (se si considera l’errore assoluto
è eref = 1.0, se si considera l’errore relativo è eref = e1).
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 18/25
9
Scelta dell’algoritmo
Non è possibile definire a priori l’algoritmo “migliore”,
ma va scelto di volta in volta in base a:
Non-linearità fisica o geometrica
Hardening o softening
Dimensione del problema
Carichi applicati in singoli punti o su intere
superfici
Non-linearità improvvise o graduali
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 19/25
Esempio di elemento geometricamente N-L:
biella in formulazione lagrangiana totale
y
v1
x
coordinate locali
v4
v2
v3
E
Interpolazione degli spostamenti:
derivate degli spostamenti:
A A
⎡1 − ξ
⎡ux ⎤ ⎢ 2
⎢u ⎥ = ⎢
⎣ y⎦ ⎢ 0
⎣⎢
∂u x −u1 + u3
=
∂x
A
⎤ ⎡ v1 ⎤
0 ⎥⎢ ⎥
v2
⎥⎢ ⎥
1 + ξ ⎥ ⎢ v3 ⎥
0
⎢ ⎥
2 ⎦⎥ ⎣⎢ v4 ⎦⎥
∂u y −u2 + u4
=
∂x
A
1+ ξ
2
0
1− ξ
2
,
trascurabile rispetto al termine lineare
∂u x 1 ⎛ ∂u x ⎞ 1 ⎛ ∂u y ⎞
+ ⎜
⎟
⎟ + ⎜
∂x 2 ⎝ ∂x ⎠ 2 ⎝ ∂x ⎠
2
Deformazioni (congruenza non-lineare): ε xx =
Sforzi (legame lineare):
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
2
σ xx = Eε xx
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 20/25
10
Rigidezza elastica e rigidezza geometrica
2
Energia interna:
2
1
1 A
EAA ⎡ ∂u x 1 ⎛ ∂u y ⎞ ⎤
⎢
⎥ =
+ ⎜
U = ∫ σ xx ε xx dV = ∫ EAε xx2 dx =
⎟
2 V
2 0
2 ⎢ ∂x 2 ⎝ ∂x ⎠ ⎥
⎣
⎦
trascurabile rispetto al
2
4
2
termine del terzo ordine
EAA ⎡⎛ ∂u x ⎞ ∂u x ⎛ ∂u y ⎞ 1 ⎛ ∂u y ⎞ ⎤
⎢⎜
=
+
⎜
⎟ + ⎜
⎟ ⎥
⎟
∂x ⎝ ∂x ⎠ 4 ⎝ ∂x ⎠ ⎥
2 ⎢⎝ ∂x ⎠
⎣
⎦
2
2
AE
EA ( v3 − v1 )
N ( v4 − v2 )
( v3 − v1 )
dove: N =
+
A
A
A
2
2
Primo teorema di Castigliano:
∂ 2U
∂vi ∂v j
kij =
⎡1
⎢
AE ⎢ 0
kT = kE + kG =
A ⎢ −1
⎢
rigidezza elastica (si ottiene
⎣0
considerando la congurenza lineare).
Matrice di rigidezza tangente
in coordinate locali:
kE
k G rigidezza geometrica (dipende dallo
0 −1 0 ⎤
⎡0 0
⎥
0 0 0 ⎥ N ⎢⎢0 1
+
0 1 0⎥ A ⎢0 0
⎥
⎢
0 0 0⎦
⎣0 −1
kE
0⎤
0 −1⎥⎥
0 0⎥
⎥
0 1⎦
0
kG
stato di sollecitazione).
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 21/25
Trasformazione di coordinate e commento
Trasformazione di coordinate:
k T glob = Q T k TQ
⎡ c s 0 0⎤
⎢−s c 0 0⎥
⎥
Q=⎢
⎢ 0 0 c s⎥
⎢
⎥
⎣ 0 0 −s c ⎦
con:
Y
⎧c = cos α
⎨
⎩ s = sin α
y
α
x
X
coordinate globali
L’approccio classico appena visto, è formulato in modo tale che l’espressione
delle forze interne abbia la forma approssimata:
fint glob ( v glob ) = k T ( v glob ) v glob
la cui derivata non è kT!
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 22/25
11
Esempio di elemento fisicamente N-L:
legge elastica bi-lineare a trazione
σ
Legame:
se ε ≤ ε 0
⎧ Eε
⎩ Eε 0 + E1 ( ε − ε 0 ) se ε > ε 0
σ =⎨
1
E
E1
1
ε0
Legame
tangente:
D=
ε
∂σ ⎧ E se ε ≤ ε 0
=⎨
∂ε ⎩ E1 se ε > ε 0
σ
Attenzione: in scarico e ricarico segue sempre lo stesso
percorso, NON è un comportamento plastico nel quale i
percorsi di scarico hanno pendenza parallela al tratto iniziale:
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
ε
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 23/25
Biella elastica bi-lineare av trazione
Congruenza (lineare):
Legame (non-lineare):
v −v
ε= 3 1
A
se ε ≤ ε 0
⎧ Eε
σ =⎨
⎩ Eε 0 + E1 ( ε − ε 0 ) se ε > ε 0
Vettore delle forze
interne in coordinate
locali (equilibrio):
⎡− N ⎤
⎡ −1⎤
⎢ 0 ⎥
⎢0⎥
fint = ⎢ ⎥ = Aσ ⎢ ⎥
⎢N ⎥
⎢1⎥
⎢ ⎥
⎢ ⎥
⎣ 0 ⎦
⎣0⎦
Matrice di rigidezza
tangente in coordinate
locali:
⎡1
⎢
∂f
∂f ∂σ ∂ε AD ⎢ 0
=
k T = int = int
∂v
∂σ ∂ε ∂v
A ⎢ −1
⎢
⎣0
2
v1
v4
v3
Linearizzazione del
legame:
D=
∂σ ⎧ E se ε ≤ ε 0
=⎨
∂ε ⎩ E1 se ε > ε 0
0 −1 0 ⎤
0 0 0 ⎥⎥
0 1 0⎥
⎥
0 0 0⎦
Segue la consueta trasformazione in coordinate globali.
Naturalmente è possibile formulare elementi con non-linearità fisiche e
geometriche contemporaneamente!
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 24/25
12
Nel prossimo capitolo
Cenni di stabilità computazionale
Università degli Studi di Firenze
Dipartimento di Ingegneria Civile
Corso di Meccanica Computazionale
Claudio Borri
Luca Salvatori
(rev. 12/06/2006) Capitolo 6: 25/25
13