MSS2 10-11 Elementi Finiti

Principi
p ed applicazioni
pp
del
metodo degli elementi finiti
Formulazione base con approccio agli
spostamenti
p
PRINCIPIO DEI LAVORI VIRTUALI
Data una certa statica: sforzi σij, forze di volume Fi e forze di superficie fj;
si dimostra che imporre la seguente equazione:
∫ σijεijdv = ∫ Fsj jdv + ∫ fjs jds + ∫ rj sjds
V
V
SF
SV
∀ cinematica congruente, cioè deformazioni εij e spostamenti si, tali che:
⎛ ∂si ∂s j ⎞
1
εij =
+
in V
2 ⎜⎜ ∂X ∂X ⎟⎟
i ⎠
⎝ j
equivale
i l ad
d iimporre lle equazioni
i i di equilibrio:
ilib i
[ si = si
su S V ]
METODO DEGLI ELEMENTI FINITI PER UN PROBLEMA
ELASTICO-LINEARE 2D
Si consideri un problema piano,
piano il cui dominio sia quello rappresentato sotto e
si voglia determinare lo stato tenso-deformativo del solido indotto dai carichi e
dai vincoli presenti.
PROCEDURA
1
1.
Suddivido il mio dominio (bidimensionale) in un certo numero di parti,
parti ad
esempio triangolari, dette elementi finiti (mesh) e numero i nodi della griglia;
2
4
1
5
3
9
8
10
6
7
11
12
14
13
Lo scopo è quello di trasformare un problema differenziale avente come incognite
dei campi (funzioni),
(funzioni) in un problema algebrico avente come incognite degli scalari.
scalari
Più precisamente si passa dall’avere come incognite primarie le funzioni sx(x,y) e
sy(x,y) ad avere come incognite gli spostamenti Ux e Uy dei nodi della griglia.
2.
Approssimo
pp
linearmente il campo
p di spostamenti
p
su ciascuno di q
questi
elementi: (si noti l’introduzione delle numerazione locale dei nodi a
livello di elemento)
⎧⎪s x ( x, y ) = a1 + b1x + c1y
⎨
⎪⎩s y ( x, y ) = a2 + b2 x + c 2 y
U3y
3
U2y
Cambio di parametri
d’interpolazione
d
interpolazione
2
U2x
U3x
U1y
1
U1x
⎧⎪s x ( x,y ) = N1 ( x,y ) U1x + N2 ( x,y ) U2x + N3 ( x,y ) U3x = Ni ( x,y ) Uix
⎨
⎪⎩s y ( x,y ) = N1 ( x,y ) U1y + N2 ( x,y ) U2y + N3 ( x,y ) U3y = Ni ( x,y ) Uiy
dove : Ni ( x, y ) =
x j ym − xm y j + ( y j − ym ) x + ( x j − xm ) y
2Δ
vale che: Ni ( x j , y j ) = δij
⇒
i
⎧s x ( x j , y j ) = Ujx
⎪
⎨
⎪⎩s y ( x j , y j ) = Ujy
i
1 xi
⎛
⎜ 2Δ = det 1 x
j
⎜
⎜
1 xm
⎝
⎞
⎟
yj
⎟
y m ⎟⎠
yi
Le tra funzioni di forma
dell’elemento triangolare
Dalla combinazione delle tre funzioni di forma si ha l’approsimazione (lineare) del campo
di spostamenti sul singolo elemento finito
sx
U1x
U3x
sx
U2x
1
3
2
In forma vettoriale ho:
⎡ s x ⎤ ⎡N1 0 N2 0 N3
⎢s ⎥ = ⎢
⎣ y ⎦ ⎣ 0 N1 0 N2 0
se
Ne
⎡ U1x ⎤
⎢U ⎥
⎢ 1y ⎥
0 ⎤ ⎢U2x ⎥
⎢
⎥
N3 ⎥⎦ ⎢U2y ⎥
⎢U3x ⎥
⎢
⎥
⎢⎣U3y ⎥⎦
Ue
Approssimazione del
campo di spostamenti
sul singolo elemento
finito funzione dei soli
spostamenti nodali
Matrice delle
funzioni di
forma
Vettore delle
componenti di
spostamento nodali
dell’elemento
Applicando l’operatore differenziale di congruenza ho:
⎡ ε x ⎤ ⎡N1,x 0 N2,x
0 N3,x
⎢ ⎥ ⎢
0 N2,y
0
⎢ ε y ⎥ = ⎢ 0 N1,y
⎢⎣ γ xy ⎥⎦ ⎢⎣N1,y N1,x N2,y N2,x N3,y
εe
Be
⎡ U1x ⎤
⎢U ⎥
1y ⎥
0 ⎤⎢
⎥ ⎢U2x ⎥
N3,y ⎥ ⎢
⎥
U
2y
⎥
N3,x ⎥⎦ ⎢
⎢U3x ⎥
⎢
⎥
⎢⎣U3y ⎥⎦
Ue
caso piano
i
⎡∂
⎢
⎧ ε x ⎫ ⎢ ∂x
⎪ ⎪
ε = ⎨ε y ⎬ = ⎢ 0
⎢
⎪γ ⎪ ⎢
⎩ xy ⎭
∂
⎢
⎣⎢ ∂y
⎤
0⎥
⎥
∂ ⎥ ⎧sx ⎫
⎨ ⎬
∂y ⎥ ⎩s y ⎭
∂⎥
⎥
∂x ⎥⎦
3.
Fase di assemblaggio: sostituisco ora nell’equazione dei lavori virtuali scritta
in forma vettoriale:
T
T
T
F = forza di volume
∫ ε σdv = ∫ s Fdv + ∫ s fds
f = forza di superficie
V
V
S
F
∑ ∫ εTe σedv = ∑ ∫ seTFedv + ∑ ∫
e
e
e
Ve
Ve
seT feds
SFe
∑ UTe ∫ BeTσedv = ∑ UeT ∫ NeTFedv + ∑ UeT ∫ NeT feds
e
e
e
Ve
Ve
SFe
Ma i gradi di libertà locali Ue di ciascun elemento li posso esprimere in funzione
di quelli globali U, attraverso le cosiddette matrice booleane di connettività:
Ue = LeU
⎡ U1x ⎤
⎢U ⎥
⎢ 1y ⎥
⎢U2x ⎥
⎢
⎥
U
⎢ 2y ⎥
⎢U ⎥
⎢ 3x ⎥
⎢⎣U3y ⎥⎦
Ue
⎡ U1x ⎤
⎢U ⎥
⎢ 1y ⎥
⎢ . ⎥
⎢
⎥
.
⎢
⎥
⎢ . ⎥
⎢
⎥
U
⎢ Nx ⎥
⎢U ⎥
⎣ Ny ⎦
U
Le è una matrice di zeri e uno
ESEMPIO:
2
1
1 3
3 2
2
5
1
• 3x2 = 6 gradi di libertà ∀ elemento
3
2
1
3
• 5x2
2 = 10 gradi di libertà
à globali
2 3 1
4
⎡ U1x ⎤
⎢U ⎥
⎢ 1y ⎥
⎢U2x ⎥
U2 = ⎢ ⎥ =
⎢U2y ⎥
⎢U ⎥
⎢ 3x ⎥
⎢⎣U3y ⎥⎦
6x1
⎡0
⎢0
⎢
⎢1
⎢
⎢0
⎢0
⎢
⎣0
0
0
0
1
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
L 2 6x10
0
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
U
⎡ U1x ⎤
0⎤ ⎢ ⎥
U1y
⎥
⎢
⎥
0
⎥ ⎢ . ⎥
0⎥ ⎢ ⎥
⎥⋅⎢ . ⎥
0⎥ ⎢ ⎥
.
0⎥ ⎢ ⎥
⎥ ⎢U5 x ⎥
0⎦ ⎢ ⎥
⎣U5 y ⎦
10x1
Segue che:
UT ∑ LTe ∫ BeTσ e dv = UT ∑ LTe ∫ NeTFe dv + UT ∑ LTe
e
[U1X
Ve
e
T
∀UT
SFe
U1X ... UNX UNX ]1x2N ∑ LTe ∫ BeTσ edv =
e
[U1X
e
Ve
∫ Ne feds
Ve
⎡i ⎤
⎢i ⎥
⎢⎥
≡ R int
⎢i ⎥
⎢⎥
⎢i ⎥
⎢⎣ i ⎥⎦ 2Nx1
U1X ... UNX UNX ]1x2N ∑ LTe ∫ NeTFedv + [U1X U1X ... UNX UNX ]1x2N ∑ LTe
e
Ve
⎡i ⎤
⎢i ⎥
⎢⎥
⎢i ⎥ ≡RFest
⎢⎥
⎢i ⎥
⎢i ⎥
⎣ ⎦2Nx1
e
∫ Ne feds
T
SFe
⎡i ⎤
⎢i ⎥
⎢⎥
f
≡ R est
⎢i ⎥
⎢⎥
⎢i ⎥
⎢⎣i ⎦⎥ 2Nx1
∀ [U1X U1X ... UNX UNX ]1x2N
Dovendo valere ∀UT congruente allora deve valere anche per la seguente scelta di UT
UT = U1T = [1 0 ... 0 0]1x2N
Da cui:
[1
0 ... 0 0]1x2N ∑ LTe ∫ BeTσ edv =
e
[1
Ve
⎡i ⎤
⎢i ⎥
⎢⎥
≡ R int
⎢i ⎥
⎢⎥
⎢i ⎥
⎢⎣i ⎥⎦ 2Nx1
0 ... 0 0]1x2N ∑ LTe ∫ NeTFedv + [1 0 ... 0 0]1x2N ∑ LTe
e
Ve
⎡i ⎤
⎢i ⎥
⎢⎥
⎢ i ⎥ ≡RFest
⎢⎥
⎢i ⎥
⎢i ⎥
⎣ ⎦2Nx1
Da cui, la prima equazione (scalare):
f
R int (1) = R Fest (1) + R est
(1)
e
∫ Ne feds
T
SFe
⎡i ⎤
⎢i ⎥
⎢⎥
f
≡ R est
⎢i ⎥
⎢⎥
⎢i ⎥
⎢⎣i ⎥⎦ 2Nx1
Adottando le altre scelte per UT segue che:
UT = UiT
f
R int ( i ) = R Fest ( i ) + R est
(i)
Alla fine tutte le equazioni
q
che p
posso scrivere p
possono essere espresse
p
in termini
vettoriali:
f
R int = R Fest + R est
O equivalentemente:
∑ LTe ∫ BeTσedv = ∑ LTe ∫ NeTFedv + ∑ LTe ∫ NeT feds
e
Ve
e
Ve
e
SFe
Equilibrio in forma debole e discretizzato per elementi finiti
Applichiamo ora le equazioni del legame costitutivo e di congruenza in forma vettoriale:
σ e = Dε e = DBeUe = DeBeL eU
S tit i
Sostituisco
nelle
ll equazioni
i i di equilibrio
ilib i d
debole
b l iin fforma di
discretizzata:
ti
t
∑ LTe ∫ BeTDeBeLeUdv = ∑ LTe ∫ NeTFedv + ∑ LTe ∫ NeT feds
e
Ve
e
Ve
e
SFe
e trovo:
∑ LTe ∫ BeTDeBedv LeU = ∑ LTe ∫ NeTFedv + ∑ LTe ∫ NeT feds
e
e
e
Ve
Ve
Ke
SFe
Fev
Fes
∪ ∫ BTeDeBedv ⋅ U = ∪ ∫ NeTFedv + ∫ NeT feds
e
e
Ve
Ke = Matrice di
rigidezza dell’elemento
finito di dimensione 6x6
per un elemento piano a
tre nodi
Ve
K
SFe
F
P
K ⋅U = P
SISTEMA RISOLVENTE
K = matrice di rigidezza
g
globale del solido di
dimensione dipendente
dal numero totale di nodi
della discretizzazione
discreti a ione
NB. Il prodotto con le matrici booleane non viene mai eseguito (sarebbe troppo
oneroso dal punto di vista computazionale); semplicemente le matrici di rigidezza e
i vettori dei carichi di ciascun elemento vengono assemblati nella matrice di
rigidezza e nel vettore dei carichi globali.
ESEMPIO:
2
1
1 3
3 2
2
5
1
• 3x2 = 6 gradi di libertà ∀ elemento
3
2
1
3
• 5x2
2 = 10 gradi di libertà
à globali
2 3 1
4
⎡ U1x ⎤
⎢U ⎥
⎢ 1y ⎥
⎢U2x ⎥
U2 = ⎢ ⎥ =
⎢U2y ⎥
⎢U ⎥
⎢ 3x ⎥
⎢⎣U3y ⎥⎦
6x1
⎡0
⎢0
⎢
⎢1
⎢
⎢0
⎢0
⎢
⎣0
0
0
0
1
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
L 2 6x10
0
0
0
0
1
0
0
0
0
0
0
1
0
0
0
0
0
0
U
⎡ U1x ⎤
0⎤ ⎢ ⎥
U1y
⎥
⎢
⎥
0
⎥ ⎢ . ⎥
0⎥ ⎢ ⎥
⎥⋅⎢ . ⎥
0⎥ ⎢ ⎥
.
0⎥ ⎢ ⎥
⎥ ⎢U5 x ⎥
0⎦ ⎢ ⎥
⎣U5 y ⎦
10x1
2
2
2
2
⎡K11
K12
K13
K14
⎢
2
2
2
K
K
K
22
23
24
⎢
2
2
⎢
K 33 K 34
K2 = ⎢
K 244
⎢
⎢
Sym
⎢
⎢⎣
K
K
K
K
K
2
15
2
25
2
35
2
45
2
55
2
2
⎡K 33
K 34
⎢
2
Tabella delle
K
44
⎢
incidenze per
⎢
l’elemento 2
⎢
⎢
⎢
T
K = L 1K 1L 1 + ⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣
K
K
K
K
K
K
2
16
2
26
2
36
2
46
2
56
2
66
⎤
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎥⎦
2
K 31
K 241
2
K11
Loc.
1x ( = 1)
1y ( = 2 )
Glob.
2X ( = 3 )
2Y ( = 4 )
2x ( = 3 )
1X ( = 1)
2y ( = 4 ) 1Y ( = 2 )
3x ( = 5 ) 4X ( = 7 )
3y ( = 6 ) 4Y ( = 8 )
2
K 32
K 242
2
K12
K 222
0
0
0
0
0
0
0
0
0
0
0
Sym
LT2 K 2 L 2
2
K 35
K 245
2
K15
2
K 25
0
0
2
K 55
2
K 36
K 246
2
K16
2
K 26
0
0
2
K 56
K 266
0
0
0
0
0
0
0
0
0
0⎤
⎥
0⎥
0⎥
⎥
0⎥
0⎥ T
⎥ + L 3K 3L 3
0⎥
0⎥
⎥
0⎥
0⎥
⎥
0 ⎥⎦
4.
Imposizione delle condizioni al contorno
Vi sono due tipi di condizioni al contorno: carichi applicati (dette naturali), i quali
compaiono nel vettore dei termini noti P; spostamenti imposti (essenziali), i quali
vengono imposti direttamente sui nodi interessati e il sistema lineare viene ridotto
alle
ll sole
l righe
i h colonne
l
non iinteressate
t
t d
da ttalili condizioni
di i i all contorno.
t
NB: Se uno spostamento è assegnato la relativa forza esterna non può essere
prescritta e rimane incognita
incognita.
ESEMPIO:
2
1
1 3
3 2
y
2
5
1
1
2 3 1
x
3
2
3
δ
4
In questo esempio: Ux1=0; Uy3= -δ; Uy4=0; Ux5=0; e il sistema risolvente diventa:
L’imposizione dei vincoli si effettua nel sistema assemblato, per prima cosa
eliminando le equazioni relative ai gradi di libertà vincolati.
⎡K11 K12
⎢
K 22
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎢⎣
K13
K 23
K 33
Sym
K14
K 24
K 34
K 44
K15
K 25
K 35
K 45
K 55
K16
K 26
K 36
K 46
K 56
K 66
K17
K 27
K 37
K 47
K 57
K 67
K 77
K18
K 28
K 38
K 48
K 58
K 68
K 78
K 88
K19
K 29
K 39
K 49
K 59
K 69
K 79
K 89
K 99
K110 ⎤ ⎡ 0 ⎤ ⎡ R1 ⎤
⎥⎢
⎥ ⎢ ⎥
K 210 ⎥ ⎢ U1y ⎥ ⎢ F2 ⎥
K 310 ⎥ ⎢U2x ⎥ ⎢ F3 ⎥
⎥⎢
⎥ ⎢ ⎥
K 410 ⎥ ⎢U2y ⎥ ⎢ F4 ⎥
K 510 ⎥ ⎢U3x ⎥ ⎢ F5 ⎥
⎥⎢
⎥=⎢ ⎥
K 610 ⎥ ⎢ −δ ⎥ ⎢R6 ⎥
K 710 ⎥ ⎢U4x ⎥ ⎢ F7 ⎥
⎥⎢
⎥ ⎢ ⎥
0
K 810 ⎥ ⎢
⎥ ⎢R 8 ⎥
K 910 ⎥ ⎢ 0 ⎥ ⎢R9 ⎥
⎥⎢
⎥ ⎢ ⎥
K1010 ⎦⎥ ⎢⎣U5y ⎥⎦ ⎢⎣F10 ⎥⎦
⎡ K12 K 22 K 23 K 24 K 25 K 26 K 27 K 28 K 29
⎢
⎢ K 31 K 32 K 33 K 34 K 35 K 36 K 37 K 38 K 39
⎢ K 41 K 42 K 43 K 44 K 45 K 46 K 47 K 48 K 49
⎢
⎢ K 51 K 52 K 53 K 54 K 55 K 56 K 57 K 58 K 59
⎢ K 71 K 72 K 73 K 74 K 75 K 76 K 77 K 78 K 79
⎢
⎢⎣K101 K102 K103 K104 K105 K106 K107 K108 K109
⎡ 0 ⎤
⎢U ⎥
⎢ 1y ⎥
K 210 ⎤ ⎢U2x ⎥ ⎡ F2 ⎤
⎥ ⎢U ⎥ ⎢ ⎥
K 310 ⎥ ⎢ 2y ⎥ ⎢ F3 ⎥
K 410 ⎥ ⎢U3x ⎥ ⎢ F4 ⎥
⎥=⎢ ⎥
⎥⎢
K 510 ⎥ ⎢ −δ ⎥ ⎢ F5 ⎥
K 710 ⎥ ⎢⎢U4x ⎥⎥ ⎢ F7 ⎥
⎥
⎢ ⎥
K1010 ⎥⎦ ⎢ 0 ⎥ ⎣F10 ⎦
⎢
⎥
0
⎢
⎥
⎢U5y ⎥
⎣
⎦
Quindi bisogna portare al secondo membro i termini noti legati al valore imposto
dello spostamento.
⎡K 22 K 23 K 24
⎢
K 33 K 34
⎢
⎢
K 44
⎢
⎢
⎢
⎢
⎢⎣
K 25
K 35
K 45
K 55
K 27
K 37
K 47
K 57
K 77
⎡ K 26 ⎤
K 210 ⎤ ⎡ U1y ⎤ ⎡ F2 ⎤
⎥⎢
⎥ ⎢ ⎥
⎢
⎥
K 310 ⎥ ⎢U2x ⎥ ⎢ F3 ⎥
K
⎢ 36 ⎥
⎢ K 46 ⎥
K 410 ⎥ ⎢U2y ⎥ ⎢ F4 ⎥
⎥⎢
⎥ = ⎢ ⎥ +δ⎢
⎥
K 510 ⎥ ⎢U3x ⎥ ⎢ F5 ⎥
K
⎢ 56 ⎥
⎢ K 76 ⎥
K 710 ⎥ ⎢U4x ⎥ ⎢ F7 ⎥
⎥⎢
⎥ ⎢ ⎥
⎢
⎥
U
F
K1010 ⎥⎦ ⎢⎣ 5y ⎥⎦ ⎢⎣ 10 ⎥⎦
⎢⎣K106 ⎦⎥
SISTEMA RISOLVENTE
con le CONDIZIONI AL CONTORNO imposte
Una volta risolto il sistema posso determinare le reazioni vincolari nel seguente modo
((ovvero usando le equazioni
q
cancellate p
prima):
)
⎡ R1 ⎤ ⎡K11
⎢R ⎥ ⎢K
⎢ 6 ⎥ = ⎢ 61
⎢R8 ⎥ ⎢K 81
⎢ ⎥ ⎢
⎣R9 ⎦ ⎣K 91
K12
K 62
K 82
K 92
K13 K14 K15 K16 K17 K18 K19
K 63 K 64 K 65 K 66 K 67 K 68 K 69
K 83 K 84 K 85 K 86 K 87 K 88 K 89
K 93 K 94 K 95 K 96 K 97 K 98 K 99
⎡ 0 ⎤
⎢U ⎥
⎢ 1y ⎥
⎢U2x ⎥
⎢
⎥
K110 ⎤ ⎢U2y ⎥
K 610 ⎥ ⎢U3x ⎥
⎥⎢
⎥
K 810 ⎥ ⎢ −δ ⎥
⎥
K1010 ⎦ ⎢U4x ⎥
⎥
⎢
⎢ 0 ⎥
⎥
⎢
0
⎢
⎥
⎢U5y ⎥
⎣
⎦
5.
Risoluzione del sistema lineare
Metodi diretti:
• metodo di eliminazione di Gauss;
•…
Metodi iterativi (calcolano la soluzione come limite di una successione di vettori):
• metodi di Richardson stazionari e non;
•…
6.
Ricostruzione della soluzione: una volta calcolata la soluzione in termini di
spostamenti nodali,
nodali posso calcolare il campo deformativo e di sforzi locali con le
seguenti relazioni già viste prima:
ε e ( x, y ) = Be ( x, y ) Ue = Be ( x, y ) LeU
σ e ( x, y ) = D ( x, y ) ε e ( x, y ) = D ( x, y ) Be ( x, y ) Ue = De ( x, y ) Be ( x, y ) LeU
INTEGRAZIONE NUMERICA:
⎛ G
⎞
f (ξi )Hi ⎜⎜ = ∑ fiHi ⎟⎟
∫−1 f (ξ)dξ ≅ ∑
i=1
⎝ i=1
⎠
1
G
Tra le tante tecniche di integrazione numerica, una delle più efficienti è quella di
Gauss, in cui la posizione dei punti di Gauss ξi e dei relativi pesi Hi sono
determinati in modo tale che con n p
punti di Gauss, un p
polinomio di ordine 2n-1
possa essere integrato esattamente.
INTEGRAZIONE GAUSSIANA:
Gli integrali, che definiscono la matrice di rigidezza e il vettore dei carichi nodali
dell’elemento, non sono mai calcolati analiticamente, ma numericamente e
solitamente con la tecnica di Gauss:
G
K e = ∫ B DeBe dv ≅ ∑ BTeiDeiBeiHi
T
e
i=1
Ve
G
F = ∫ N F dv ≅ ∑ NTeiFeiHi
V
e
T
e e
Ve
Fes =
i=1
G
T
T
N
f
ds
≅
N
∫ e e ∑ ei feiHi
SFe
i=1
L’indice i è un indice che va da 1 al numero di punti di Gauss utilizzati per calcolare
numericamente l’integrale. In genere per un elemento triangolare si usano G=1 o
G=3 punti di Gauss. I valori di Hi sono i corrispondenti pesi.
L’ELEMENTO a QUATTRO NODI
Altri esempi di ELEMENTI FINITI PIANI e delle loro TECNICHE d’INTEGRAZIONE: