Elementi finiti - Parte III - POLITECNICO DI TORINO

getto
didattica in ret
progetto
didattica in rete
Elementi finiti
Parte III
A. Gugliotta
Politecnico di Torino, maggio 2002
Dipartimento di Meccanica
otto editore
ELEMENTI FINITI
Parte III
A. GUGLIOTTA
P OLITECNICO DI TORINO
WWW. POLITO . IT
INDICE – III
5.
ELEMENTI ISOPARAMETRICI ........................................135
5.1 SISTEMA DI RIFERIMENTO LOCALE E FUNZIONI DI FORMA .... 135
Elementi monodimensionali ...........................................................136
Elementi triangolari ........................................................................139
Elementi quadrangolari ...................................................................142
Elementi solidi tetraedri ..................................................................147
Elementi solidi parallelepipedi ........................................................ 149
5.2 CALCOLO DELLA MATRICE DI RIGIDEZZA ............................... 152
5.3 CARICHI NODALI EQUIVALENTI ............................................. 160
Carichi di linea ...............................................................................160
Carichi di volume ...........................................................................164
Effetti termici.................................................................................. 164
5.4 CALCOLO DELLE TENSIONI .................................................... 165
5.5 PROBLEMI RELATIVI AGLI ELEMENTI ISOPARAMETRICI ........... 169
Scelta dell'ordine di integrazione .....................................................169
Distorsione degli elementi............................................................... 171
Integrazione selettiva e modi incompatibili .....................................174
Autovalori della matrice di rigidezza ................................................180
Patch test ........................................................................................187
6.
SOLIDI ASSIALSIMMETRICI ............................................189
6.1 INTRODUZIONE ..................................................................... 189
6.2 CARICO ASSIALSIMMETRICO ................................................... 190
6.3 CARICO NON ASSIALSIMMETRICO .......................................... 194
i
7.
PIASTRE INFLESSE ............................................................201
7.1 RICHIAMI DELLA TEORIA DELLE PIASTRE INFLESSE ................. 201
Sistemi di riferimento ......................................................................201
Tensioni e carichi per unità di lunghezza .........................................202
7.2 IPOTESI DI KIRCHOFF ............................................................. 203
Spostamenti e deformazioni ............................................................203
Relazioni tensioni-deformazioni ......................................................204
7.3 IPOTESI DI MINDLIN .............................................................. 205
Spostamenti e deformazioni ............................................................205
Relazioni tensioni-deformazioni ...................................................... 206
7.4 FORMULAZIONE DI RIGIDEZZA NELL’IPOTESI DI KIRCHOFF .. 207
7.5 LA TRAVE DI TIMOSHENKO .................................................... 209
Elemento lineare .............................................................................210
Elemento parabolico .......................................................................213
7.6 MODELLI DISCRETI
(LEGAME DISCRETO TRA FLESSIONE E TAGLIO) ...................... 225
7.7 L’ELEMENTO PIASTRA DI MINDLIN ........................................ 229
Matrice di rigidezza......................................................................... 229
Carico nodale equivalente a carico superficiale distribuito ............... 231
ii
5. ELEMENTI ISOPARAMETRICI
5.1 SISTEMA DI RIFERIMENTO LOCALE E FUNZIONI DI FORMA
Per semplificare la caratterizzazione di elementi di forma geometrica complessa si
introduce il concetto di elemento di riferimento, di forma geometrica semplice,
definito in uno spazio naturale e adimensionale.
La formulazione si basa sul fatto di creare una corrispondenza biunivoca tra un
elemento di forma qualsiasi nel sistema cartesiano (x , y , z ) e un elemento di
forma semplice nel sistema naturale, o locale, ( x, , z ).
Un generico elemento nel sistema fisico (x, y, z) viene cioè pensato come la trasformazione di un elemento di forma semplice, e sempre lo stesso, nel sistema
locale (x, , z ); nei casi mono-, bi- e tri-dimensionale si hanno rispettivamente
le seguenti relazioni:
n
n
n
x =
x =
 ni (x )xi
i=1
 ni (x ,h)xi
i=1
n
y =
x =
y =
 ni (x ,h)yi
i=1
 ni (x ,h ,z )xi
i=1
n
 ni (x ,h ,z )yi
5.1
i=1
n
z =
 ni (x ,h ,z ) zi
i=1
dove xi , yi , zi (i = 1, n) sono le coordinate cartesiane degli n nodi dell'elemento,
ni sono le funzioni di trasformazione o funzioni di forma definite nel sistema di
coordinate naturale (x, , z ) dell'elemento.
Le funzioni di forma ni sono in numero pari ai nodi dell'elemento e variano tra 0
e 1; esse assumono valore 1 in corrispondenza del nodo i di coordinate (xi , yi , zi )
e valore nullo in corrispondenza degli altri nodi dell'elemento.
135
ELEMENTI ISOPARAMETRICI
5.1.1 Elementi monodimensionali
Si consideri un elemento a due nodi di lunghezza l (fig. 5.1):
Fig. 5.1 – Elemento a 2 nodi: sistema di riferimento naturale.
la posizione di un punto P sul segmento può essere individuata dalla sua coordinata fisica xP oppure dalle misure dei segmenti l1 ed l2 o dalle loro misure adimensionali:
l
L 1 = ---1
l
l
L 2 = ---2
l
5.2
essendo naturalmente:
L1 + L2 = 1
5.3
In corrispondenza del nodo 1 si ha L1 = 1, L2 = 0, mentre nel nodo 2 è L1 = 0,
L2 = 1; si può allora descrivere la coordinata x di un punto nella forma:
x = L1 x1 + L2 x2
5.4
Fig. 5.2 – Elemento a 2 nodi: sistema di riferimento naturale.
Volendo esprimere la coordinata fisica x in funzione di una sola coordinata naturale x, variabile tra 0 e 1 (fig. 5.2) si effettua il cambio di variabile:
L1 = 1 – x
136
L2 = x
5.5
ELEMENTI ISOPARAMETRICI
e quindi:
x = ( 1 – x )x 1 + xx 2
5.6
È usuale riferire l'origine del sistema naturale nel baricentro geometrico dell'elemento (fig. 5.3), con la coordinata naturale x variabile tra –1 e 1. Si avrà allora:
1–x
L 1 = ----------2
1+x
L 2 = -----------2
5.7
Fig. 5.3 – Elemento a 2 nodi: sistema di riferimento naturale.
1+x
1–x
x = -----------x 1 + ------------ x 2
2
2
5.8
x = n1 x1 + n2 x2
5.9
avendo definito le funzioni di forma n1 e n 2:
1–x
n 1 = ----------2
1+x
n 2 = -----------2
5.10
La figura 5.4 mostra l'andamento delle funzioni di forma n1 e n2.
Fig. 5.4 – Funzioni di forma per l'elemento monodimensionale a 2 nodi.
Per un elemento a tre nodi (fig. 5.5) la coordinata x è data da:
x = n1 x1 + n2 x2 + n3 x3
5.11
137
ELEMENTI ISOPARAMETRICI
Fig. 5.5 – Elemento monodimensionale a 3 nodi.
dove:
n 1 = L 1 ( 2L 1 – 1 )
n 2 = L 2 ( 2L 2 – 1 )
n 3 = 4L 1 L 2
5.12
e ponendo l'origine del sistema di riferimento naturale nel baricentro geometrico
dell'elemento si ha:
1
n 1 = – --- x ( 1 – x )
2
1
n 2 = --- x ( 1 + x )
2
n3 = 1 – x 2
5.13
Fig. 5.6 – Funzioni di forma per l'elemento monodimensionale a 3 nodi.
La figura 5.6 mostra l'andamento delle funzioni di forma. Le tabelle 5.1 e 5.2
riassumono le funzioni di forma e le loro derivate rispettivamente per gli elementi monodimensionali a 2 nodi e a 3 nodi.
Tab. 5.1 – Funzioni di forma e loro derivate per l'elemento monodimensionale a 2 nodi
138
NODO
1
2
n
1
--- ( 1 – x )
2
1
--- ( 1 + x )
2
∂n § ∂x
1
– --2
1
--2
ELEMENTI ISOPARAMETRICI
Tab. 5.2 – Funzioni di forma e loro derivate per l'elemento monodimensionale a 3 nodi
NODO
1
2
3
n
1
– --- x ( 1 – x )
2
1
--- x ( 1 + x )
2
( 1 – x )2
∂n § ∂x
1
– --- + x
2
1
--- + x
2
– 2x
5.1.2 Elementi triangolari
Nel caso di elementi triangolari le funzioni interpolatrici sono definite in un
sistema di riferimento naturale basato sulle coordinate d'area.
La figura 5.7 illustra un generico elemento triangolare individuato dai tre nodi 1,
2, 3; i lati del triangolo sono individuati da numeri corrispondenti ai vertici
opposti.
Un punto P all'interno dell'elemento lo suddivide in tre triangoli di aree A1, A2,
A3 e può quindi essere individuato, oltre che dalle coordinate cartesiane (x, y),
anche dalle misure delle aree A1, A2, A3 o dai loro rapporti rispetto all'area totale
del triangolo, A.
A
L 1 = -----1A
A
L 2 = -----2A
A
L 3 = -----3A
5.14
Fig. 5.7 – Elemento triangolare a 3 nodi.
I tre rapporti L1, L2, L3 rappresentano le coordinate naturali o d'area per un
triangolo (fig. 5.8). Chiaramente deve essere:
L1 + L2 + L3 = 1
5.15
139
ELEMENTI ISOPARAMETRICI
Fig. 5.8 – Elemento triangolare a 3 nodi: sistema di riferimento naturale.
Come già visto per il caso monodimensionale, le coordinate (x, y) di un punto P
potranno essere descritte da:
x = L1 x1 + L2 x2 + L3 x3
5.16
y = L1 y1 + L2 y2 + L3 y3
Effettuando un cambio di variabile:
L1 = 1 – x – h
L2 = x
L3 = h
5.17
le coordinate x, y divengono:
x = ( 1 – x – h )x 1 + x x 2 + h x 3
5.18
y = ( 1 – x – h )y 1 + x y 2 + h y 3
x = n1 x1 + n2 x2 + n3 x3
5.19
y = n1 y1 + n2 y2 + n3 y3
avendo definito le funzioni di forma n1, n2 e n3:
n1 = 1 – x – h
n2 = x
n3 = h
La figura 5.9 mostra l'andamento della funzione di forma n1.
140
5.20
ELEMENTI ISOPARAMETRICI
Fig. 5.9 – Elemento triangolare a 3 nodi: funzione di forma n1.
Per un elemento triangolare a 6 nodi (fig. 5.10) si ricavano le seguenti funzioni:
n 1 = L 1 ( 2L 1 – 1 ) = 1 – 3 ( x + h ) + 2 ( x + h ) 2
n 4 = 4L 1 L 2 = 4x ( 1 – x – h )
n 2 = L 2 ( 2L 2 – 1 ) = x ( 2x – 1 )
n 5 = 4L 2 L 3 = 4xh
n 3 = L 3 ( 2L 3 – 1 ) = h ( 2h – 1 )
n 6 = 4L 1 L 3 = 4h ( 1 – x – h )
5.21
Fig. 5.10 – Elemento triangolare a 6 nodi.
La figura 5.11 mostra l'andamento delle funzioni di forma relative rispettivamente ad
un nodo d'angolo e ad un nodo intermedio per un elemento a sei nodi.
Fig. 5.11 – Elemento triangolare a 6 nodi: funzioni di forma n1e n5.
141
ELEMENTI ISOPARAMETRICI
Le tabelle 5.3 e 5.4 illustrano le funzioni di forma e le loro derivate rispettivamente per gli elementi triangolari a 3 nodi e a 6 nodi.
Tab. 5.3 – Funzioni di forma e loro derivate per l'elemento triangolare a 3 nodi
NODO
1
2
3
n
1–x–h
x
h
∂n § ∂x
–1
1
0
∂n § ∂h
–1
0
1
Tab. 5.4 – funzioni di forma e loro derivate per l'elemento triangolare a 6 nodi
NODO
1
2
3
4
5
6
n
l ( 2l – 1 )
x ( 2x – 1 )
h ( 2h – 1 )
4xl
4xh
4hl
∂n § ∂x
1 – 4l
4x – 1
0
4(l – x)
4h
– 4h
∂n § ∂h
1 – 4l
0
4h – 1
– 4x
4x
4(l – h)
l = 1–x–h
5.1.3 Elementi quadrangolari
La figura 5.12 mostra il sistema di coordinate naturali per un elemento quadrangolare. L'origine del sistema di riferimento naturale è nel baricentro geometrico
dell'elemento, mentre gli assi locali x ed intersecano i lati del quadrilatero nella
loro mezzeria.
Fig. 5.12 – Elemento quadrangolare a 4 nodi.
142
ELEMENTI ISOPARAMETRICI
Le funzioni di forma per un elemento quadrangolare a quattro nodi sono date da:
1
n 1 = --- ( 1 – x ) ( 1 – h )
4
1
n 3 = --- ( 1 + x ) ( 1 + h )
4
1
n 2 = --- ( 1 + x ) ( 1 – h )
4
1
n 4 = --- ( 1 – x ) ( 1 + h )
4
5.22
Esse sono il risultato di una interpolazione bilineare, che si può descrivere come
segue. Sul lato 1-2 la trasformazione lineare che fornisce l'ascissa x o l'ordinata y
(o qualsiasi altra funzione f variabile linearmente al variare di x) è data dalla 5.8:
– x 1 + x Ï f1 ¸
f 12 = 1---------------------- Ì ý
2
2 Ó f2 þ
5.23
Sul lato 3-4, allo stesso modo:
+ x 1 – x Ï f3 ¸
f 34 = 1----------- ----------- Ì ý
2 Ó f4 þ
2
5.24
Considerando allora una data ascissa c , e interpolando linearmente la funzione
f secondo l'ordinata h :
– h 1 + h Ï f 12 ¸
f = 1----------ý
- ------------- Ì
2 Ó f 34 þ
2
5.25
Sostituendo le 5.23 e 5.24 nella 5.25 e definendo un unico vettore dei valori
nodali fi , si ha:
1–x 1+x
----------- ------------ 0
0
2
2
–h 1+h
f = 1----------- ------------2
2
1+x 1–x
0
0 ------------ ----------2
2
Ï
Ô
Ô
Ì
Ô
Ô
Ó
f1 ¸
Ô
f2 Ô
ý
f3 Ô
Ô
f4 þ
Ï f1 ¸
Ô Ô
Ô f2 Ô
(
1
+
x
)
(
1
+
h
)
(
1
–
x
)
(
1
+
h
)
(
1
–
x
)
(
1
–
h
)
(
1
+
x
)
(
1
–
h
)
f = ---------------------------------- ----------------------------------- ----------------------------------- ----------------------------------- Ì ý
4
4
4
4
Ô f3 Ô
Ôf Ô
Ó 4þ
5.26
5.27
La figura 5.13 mostra l'andamento della funzione di forma n1.
143
ELEMENTI ISOPARAMETRICI
Fig. 5.13 – Elemento quadrangolare: funzione di forma n 1 .
Fig. 5.14 – Elementi quadrangolari parabolici.
La figura 5.14 illustra due elementi quadrangolari a lati curvi, entrambi con tre
nodi per lato; il primo elemento è detto elemento Serendipity; il secondo, avente
un nodo anche nel baricentro geometrico, è detto elemento Lagrange.
Nel secondo caso si utilizza un polinomio completo di secondo grado, mentre
nel primo viene a mancare il termine quadratico misto x 2 y 2.
Le funzioni di forma dei due elementi possono essere ricavate mediante sovrapposizione delle singole componenti, come illustrato in figura 5.15.
La funzione di forma relativa al nodo 5 (fig. 5.15b) può essere ottenuta come
prodotto di una parte quadratica (1 – x2) nella direzione del lato e di una parte
lineare (1 – )/2 nella direzione perpendicolare al lato:
1
n 5 = --- ( 1 – x 2 ) ( 1 – h )
2
5.28
In modo analogo si ricava la funzione di forma n8 relativa al nodo 8 (fig. 5.15c):
1
n 6 = --- ( 1 – x ) ( 1 – h 2 )
2
144
5.29
ELEMENTI ISOPARAMETRICI
Fig. 5.15 – Elementi quadrangolari parabolici: funzioni di forma.
La funzione di forma relativa al nodo 1 deve essere tale da avere valore nullo,
oltre che nei punti 2, 3, 4, 6, 7, anche nei punti 5 e 8. Essa può quindi essere
ottenuta sottraendo alla funzione bilineare (fig. 5.15a), che fornisce valori pari a
0.5 in corrispondenza dei nodi 5 e 8, metà delle funzioni relative ai nodi 5 e 8,
avendosi:
1
1
1
1
n 1 = --- ( 1 – x ) ( 1 – h ) – --- n 5 – --- n 8 = --- [ ( 1 – x ) ( 1 – h ) ( – x – h – 1 ) ]
4
2
2
4
5.30
La figura 5.15d mostra la funzione di forma relativa al nodo 1 per un elemento a
8 nodi di tipo Serendipity. In modo analogo si ottengono le altre funzioni di
forma relative ai nodi d'angolo. Questo modo di procedere suggerisce inoltre la
possibilità di ottenere elementi con numero di nodi variabile per lato, come ad
esempio indicato in figura 5.16.
Fig. 5.16 – Elemento quadrangolare con numero di nodi variabile.
145
ELEMENTI ISOPARAMETRICI
La tabella 5.5 riassume le funzioni di forma per elementi con numero di nodi
variabile tra 4 e 8. Le tabelle 5.6 e 5.7 riassumono le derivate delle funzioni di
forma rispettivamente per gli elementi quadrangolari a 4 nodi e a 8 nodi.
Tab. 5.5 – Funzioni di forma per elementi quadrangolari con numero di nodi variabile
tra 4 e 8
N
FUNZIONE
+ (SE C'È 5)
1
1
--- ( 1 – x ) ( 1 – h )
4
1
– --- n 5
2
2
1
--- ( 1 + x ) ( 1 – h )
4
1
– --- n 5
2
3
1
--- ( 1 + x ) ( 1 + h )
4
4
1
--- ( 1 – x ) ( 1 + h )
4
+ (SE C'È 6)
+ (SE C'È 7)
+ (SE C'È 8)
1
– --- n 8
2
1
– --- n 6
2
1
– --- n 7
2
1
– --- n 6
2
1
– --- n 7
2
1
– --- n 8
2
1
--- ( 1 – x 2 ) ( 1 – h )
4
5
1
--- ( 1 + x ) ( 1 – h 2 )
4
6
1
--- ( 1 – x 2 ) ( 1 + h )
4
7
1
--- ( 1 – x ) ( 1 – h 2 )
4
8
Tab. 5.6 – Derivate delle funzioni di forma per l'elemento quadrangolare a 4 nodi
NODO
1
2
3
4
∂n § ∂x
(1 – h)
----------------4
(1 – h)
----------------4
(1 + h)
----------------4
(1 + h)
----------------4
∂n § ∂h
(1 – x)
---------------4
(1 + x)
----------------4
(1 + x)
----------------4
(1 – x)
---------------4
146
ELEMENTI ISOPARAMETRICI
Tab. 5.7 – Derivate delle funzioni di forma per l'elemento quadrangolare a 8 nodi
NODO
1
2
3
4
∂n § ∂x
( 1 – h ) ( 2x + h )
--------------------------------------4
( 1 – h ) ( 2x – h )
-------------------------------------4
( 1 + h ) ( 2x + h )
--------------------------------------4
( 1 + h ) ( 2x – h )
--------------------------------------4
∂n § ∂h
( 1 – x ) ( 2h + x )
-------------------------------------4
( 1 + x ) ( 2h – x )
-------------------------------------4
( 1 + x ) ( 2h + x )
-------------------------------------4
( 1 – x ) ( 2h – x )
------------------------------------4
NODO
5
6
7
8
∂n § ∂x
– x (1 – h)
1 – h2
--------------2
–x ( 1 + h )
1 – h2
– --------------2
∂n § ∂h
1 – x2
– -------------2
–( 1 – x ) h
1 – x2
-------------2
–( 1 – x ) h
5.1.4 Elementi solidi tetraedri
Il sistema di riferimento locale è descritto in termini di coordinate di volume L1,
L2, L3, L4, in modo analogo al caso degli elementi triangolari piani:
V
L 1 = -----1V
V
L 2 = -----2V
V
L 3 = -----3V
V
L 4 = -----4V
5.31
Fig. 5.17 – Elemento tetraedro lineare.
Le funzioni di forma e le loro derivate per l'elemento tetraedro lineare a 4 nodi
(fig. 5.17) sono riassunte in tabella 5.8.
147
ELEMENTI ISOPARAMETRICI
Tab. 5.8 – Funzioni di forma e loro derivate per l'elemento tetraedro a 4 nodi
NODO
1
2
3
4
n
1–x–h–z
x
h
z
∂n § ∂x
–1
1
0
0
( ∂n ) § ( ∂h )
–1
0
1
0
∂n § ∂z
–1
0
0
1
La tabella 5.9 riassume le funzioni di forma e le loro derivate per l'elemento
tetraedro quadratico a 10 nodi (fig. 5.18).
Tab. 5.9 – Funzioni di forma e loro derivate per l'elemento triangolare a 10 nodi
NODO
n
∂n § ∂x
∂n § ∂h
∂h § ∂z
1
l ( 2l – 1 )
1 – 4l
1 – 4l
1 – 4l
2
x ( 2x – 1 )
4x – 1
0
0
3
h ( 2h – 1 )
0
4h – 1
0
4
z ( 2z – 1 )
0
0
4z – 1
5
4xl
4(l – x)
– 4x
– 4x
6
4xh
4h
4x
0
7
4hl
– 4h
4(l – h)
– 4h
8
4zl
– 4z
– 4z
4(l – z)
9
4xz
4z
0
4x
10
4hl
0
4z
4h
l = 1–x–h–z
148
ELEMENTI ISOPARAMETRICI
Fig. 5.18 – Elemento tetraedro parabolico.
5.1.5 Elementi solidi parallelepipedi
Le funzioni di forma per gli elementi solidi parallelepipedi possono essere ricavate a partire da quelle per gli elementi quadrangolari piani aggiungendo la terza
coordinata .
Fig. 5.19 – Elemento solido (brick).
In modo compatto le funzioni di forma e le loro derivate per l'elemento lineare a
8 nodi (fig. 5.19) possono essere scritte come:
149
ELEMENTI ISOPARAMETRICI
1
n i = --- ( 1 + xx i ) ( 1 + hh i ) ( 1 + zz i )
8
∂n i
1
= --- x i ( 1 + hh i ) ( 1 + zz i )
8
∂x
5.32
∂n i
1
= --- h ( 1 + xx i ) ( 1 + zz i )
8 i
∂h
∂n i
1
= --- z i ( 1 + xx i ) ( 1 + zz i )
8
∂z
dove i valori di x i ,h i ,z i sono riportati in tabella (tab. 5.10).
Tab. 5.10 – Valori di x i , h i , z i per i nodi d'angolo
NODO
1
2
3
4
5
6
7
8
x
–1
1
1
–1
–1
1
1
–1
h
–1
–1
1
1
–1
–1
1
1
z
–1
–1
–1
–1
1
1
1
1
Fig. 5.20 – Elemento solido a 20 nodi.
Per l'elemento quadratico a 20 nodi (fig. 5.20) le funzioni di forma e le loro derivate sono:
– per i nodi sui vertici dell'elemento (1, 2, 3, 4, 5, 6, 7, 8), (tab. 5.10):
150
ELEMENTI ISOPARAMETRICI
1
n i = --- ( 1 + xx i ) ( 1 + hh i ) ( 1 + zz i ) ( – 2xx i + hh i + zz i )
8
∂n i
1
= --- x i ( 1 + hh i ) ( 1 + zz i ) ( – 1 + 2xx i + hh i + zz i )
8
∂x
∂n i
1
= --- h ( 1 + xx i ) ( 1 + zz i ) ( – 1 + xx i + 2hh i + zz i )
8 i
∂h
5.33
∂n i
1
= --- z i ( 1 + xx i ) ( 1 + zz i ) ( – 1 + xx i + hh i + 2zz i )
8
∂z
– per i nodi con x = 0 (9, 11, 13, 15):
NODO
9
11
13
15
h
–1
1
–1
1
z
–1
–1
1
1
1
n i = --- ( 1 – x 2 ) ( 1 + hh i ) ( 1 + zz i )
4
∂n i
1
= – --- x ( 1 + hh i ) ( 1 + zz i )
2
∂x
5.34
∂n i
1
= --- h i ( 1 – x 2 ) ( 1 + zz i )
4
∂h
∂n i
1
= --- z i ( 1 – x 2 ) ( 1 + hh i )
4
∂z
– per i nodi con = 0 (10, 12, 14, 16):
NODO
10
12
14
16
x
1
–1
1
–1
z
–1
–1
1
1
151
ELEMENTI ISOPARAMETRICI
1
n i = --- ( 1 + xx i ) ( 1 – h 2 ) ( 1 + zz i )
4
∂n i
1
= --- x i ( 1 – h 2 ) ( 1 + zz i )
4
∂x
5.35
∂n i
1
= – --- h ( 1 + xx i ) ( 1 + zz i )
∂h
2
∂n i
1
= --- z i ( 1 + xx i ) ( 1 – h 2 )
4
∂z
– per i nodi con = 0 (17, 18, 19, 20):
NODO
1
18
19
20
x
1
1
–1
–1
h
–1
1
1
–1
1
n i = --- ( 1 + xx i ) ( 1 + hh i ) ( 1 – z 2 )
4
∂n i
1
= --- x i ( 1 + hh i ) ( 1 – z 2 )
4
∂x
∂n i
1
= --- h i ( 1 + xx i ) ( 1 – z 2 )
4
∂h
5.36
∂n i
1
= – --- z ( 1 + xx i ) ( 1 + hh i )
∂z
2
5.2 CALCOLO DELLA MATRICE DI RIGIDEZZA
Le leggi dello spostamento sono definite, così come per le coordinate, nel sistema
di riferimento naturale; nei casi mono-, bi- e tri-dimensionale si ha, rispettivamente:
152
ELEMENTI ISOPARAMETRICI
m
m
u =
m
u =
 ni (x )ui
i=1
 ni (x ,h )ui
i=1
m
v =
u =
i=1
m
v =
 ni (x ,h )vi
i=1
 ni (x ,h ,z )ui
 ni (x ,h ,z )vi
5.37
i=1
m
w =
 ni (x ,h ,z )wi
i=1
dove ui , vi , wi (i = 1,m) sono gli spostamenti nodali dell'elemento, ni sono le
funzioni di forma definite nel sistema di coordinate naturali (, , ) dell'elemento.
È da notare che le leggi per definire gli spostamenti possono essere diverse da
quelle utilizzate per definire la geometria dell'elemento; nelle espressioni degli
spostamenti la sommatoria è estesa all'indice m, mentre nelle espressioni della
geometria la sommatoria è estesa all'indice n.
È tuttavia utilizzata la formulazione isoparametrica, nella quale è utilizzato lo
stesso numero di parametri (n = m) per definire sia la geometria sia gli spostamenti.
Altri tipi di formulazione danno luogo a elementi subparametrici e superparametrici; nei primi la legge per definire la geometria ha un grado minore di quella
utilizzata per definire gli spostamenti (n < m), nei secondi ha un grado maggiore
(n > m).
Riferendosi per semplicità al caso bidimensionale si avrà quindi:
m
u =
 ni (x ,h )ui
i=1
m
v =
5.38
 ni (x ,h )vi
i=1
Per il calcolo della matrice di rigidezza [k] dell'elemento:
[k] =
ÚV [ b ]T [ E ] [ b ] dV = ÚA [ b ]T [ E ] [ b ]h dA
5.39
è necessario calcolare la matrice [b] che lega le deformazioni agli spostamenti:
{e} = [b]{s}
5.40
avendo definito il vettore {e} come:
{ e } T = { e xx e yy g xy }
5.41
153
ELEMENTI ISOPARAMETRICI
con:
e xx =
∂u
∂x
e yy =
∂v
∂y
g xy =
∂u ∂v
+
∂y ∂x
5.42
Sostituendo agli spostamenti u e v le espressioni delle funzioni di forma si ha:
n
e xx =
i=1
n
e yy =
∂n i
 ∂ x ui
∂n i
 ∂ y vi
i=1
n
g xy =
∂n i
5.43
∂n i
 ÊË ∂ y ui + ∂ x viˆ¯
i=1
e in forma matriciale, esplicitando il solo termine i-esimo:
º
º
∂n i
∂x
{e} = º
º
0
∂n i
∂y
º
º
º
∂n i
∂y
∂n i
∂x
º
0
º
Ï
Ô
º Ô
Ô
Ô
º Ì
Ô
Ô
Ô
º Ô
Ó
¸
Ô
Ô
Ô
Ô
ý
vi Ô
Ô
º Ô
º Ôþ
º
º
ui
5.44
Sono quindi da calcolare le derivate delle funzioni di forma ni rispetto alle variabili x, y, mentre dette funzioni dipendono esplicitamente da x, . Per la generica
funzione si può scrivere:
∂n i ∂x ∂n i ∂y
∂n i
=
+
∂x ∂x ∂y ∂x
∂x
∂n i ∂x ∂n i ∂y
∂n i
=
+
∂x ∂h ∂y ∂h
∂h
e in forma matriciale:
154
5.45
ELEMENTI ISOPARAMETRICI
Ï
Ô
Ô
Ì
Ô
Ô
Ó
∂n i ¸
∂x ∂y
Ô
∂x Ô
∂x ∂x
ý =
∂n i Ô
∂x ∂y
Ô
∂h ∂h
∂h þ
Ï
Ô
Ô
Ì
Ô
Ô
Ó
Ï
∂n i ¸Ô
Ô
Ô
∂x Ô
ý = [ J ]Ì
Ô
∂n i Ô
Ô
Ô
∂h þ
Ó
Ï
Ô
Ô
Ì
Ô
Ô
Ó
∂n i ¸
Ô
∂x Ô
ý
∂n i Ô
Ô
∂y þ
5.46
∂n i ¸Ô
∂x Ô
ý
∂n i Ô
Ô
∂y þ
5.47
dove [ J ] è detta matrice jacobiana e mette in relazione le derivate nelle coordinate naturali a quelle nel piano fisico:
∂x ∂y
[ J ] = ∂x ∂x
∂x ∂y
∂h ∂h
5.48
e i coefficienti della matrice sono dati da:
∂x
=
∂x
∂x
=
∂h
n
Â
∂n i
x
∂x i
Â
∂n i
x
∂h i
i=1
n
i=1
∂y
=
∂x
∂y
=
∂h
n
∂n i
 ∂x
i=1
n
Â
i=1
yi
∂n i
y
∂h i
5.49
Ammesso che l'inversa della matrice jacobiana esista, cioè che si abbia una relazione biunivoca tra le coordinate naturali e quelle fisiche dell'elemento, si ha:
Ï
Ô
Ô
Ì
Ô
Ô
Ó
Ï
∂n i ¸
Ô
Ô
Ô
∂x Ô
ý = [ J ] –1 Ì
Ô
∂n i Ô
Ô
Ô
∂y þ
Ó
∂n i ¸
Ô
∂x Ô
ý
∂n i Ô
Ô
∂h þ
5.50
155
ELEMENTI ISOPARAMETRICI
Ï
Ô
Ô
Ì
Ô
Ô
Ó
∂n i ¸
Ô
∂x Ô
1
ý = --------------det [ J ]
∂n i Ô
Ô
∂y þ
Ï
∂y
∂y Ô ∂n i
– -----∂ h ∂x Ô ∂ x
Ì
∂x ∂x Ô ∂n i
– -----Ô
∂h ∂ x Ó ∂ h
¸
Ô
Ô
ý
Ô
Ô
þ
5.51
È così possibile valutare i vari coefficienti della matrice [b] e procedere quindi al
calcolo della matrice di rigidezza [k], che non può, in generale, essere fatto utilizzando una integrazione analitica. Si deve quindi ricorrere ad uno schema di integrazione numerica.
Prima però di eseguire l'integrazione bisogna esprimere anche il termine dA in
funzione delle variabili naturali x, ; la figura 5.21 mostra un elemento di area
dA, compreso tra le linee di coordinate x e x + dx e e + d. I lati del parallelogramma così individuato valgono:
Ï
Ô
Ô
{a} = Ì
Ô
Ô
Ó
¸
∂x Ô
∂x Ô
ýd x
∂y Ô
∂x Ô
þ
Ï ∂x ¸
Ô
Ô
Ô ∂h Ô
{b} = Ì
ý dh
∂y
Ô
Ô
Ô ∂h Ô
Ó
þ
5.52
Fig. 5.21 – Elemento di area dA.
L'area dA è data dal prodotto vettoriale di {a} per {b}:
∂x ∂y ∂y ∂x ˆ
dA = { a } ∑ { b } = Ê
–
d xd h = det [ J ] d x d h 5.53
Ë ∂ x ∂ h ∂ x ∂ h¯
e la matrice di rigidezza [k]:
156
ELEMENTI ISOPARAMETRICI
[k] =
ÚA
[ b ] T [ E ] [ b ]h dA =
1
1
Ú–1 Ú–1 [ b ]T [ E ] [ b ]det [ J ] d x d h
5.54
Uno degli schemi di integrazione numerica più utilizzati è il metodo di Gauss,
secondo il quale si ha, supponendo costante lo spessore h dell'elemento:
m
[k] =
m
  wi wj [ b ]ij [ E ] [ b ]ij det [ J ]ij h
T
5.55
i = 1j = 1
dove m è l'ordine di integrazione, wi, wj sono i pesi di integrazione, [b]ij e
det [ J ]ij sono calcolati nei punti di integrazione xi, j.
ESEMPIO 5.1
Calcolare le matrici di rigidezza nel caso di elementi monodimensionali di tipo
lineare e parabolico.
Elemento asta lineare
Nel caso di un elemento asta isoparametrico a due nodi le espressioni della
geometria e degli spostamenti sono date da:
1+x
1–x
x = -----------x 1 + ------------ x 2
2
2
1+x
1–x
u = -----------u 1 + ------------ u 2
2
2
5.56
La deformazione nell'elemento è data da:
e xx =
dn 1 dn 2 Ï u 1 ¸
du
=
Ì ý = [b]{s}
dx
d x d x Ó u2 þ
5.57
e:
dn i
dx
=
dx
dx
– 1 dn
i
dx
= [ J ] –1
dn i
dx
5.58
dove [ J ] è la matrice jacobiana della trasformazione. Le derivate delle funzioni di forma rispetto a x sono:
dn 1
1
= – --dx
2
dn 2
1
= --dx
2
5.59
essendo inoltre l = x 2 – x 1 la lunghezza dell'elemento si ha:
dn 2
dn 1
1
1
dx
l
=
x +
x = – --- x 1 + --- x 2 = --dx 1 dx 2
2
2
dx
2
5.60
157
ELEMENTI ISOPARAMETRICI
per cui:
2
[ J ] –1 = --l
l
[ J ] = --2
l
det [ J ] = --2
5.61
e la matrice [ b ] :
[ b ] = – 1--l
1
--l
5.62
La matrice di rigidezza, considerando costanti l'area A ed il modulo elastico
E è:
l
1
l
[ k ] = EA Ú [ b ] T [ b ] dx = EA Ú [ b ] T [ b ] --- d x
2
0
–1
5.63
EA 1 – 1
[ k ] = ------l –1 1
5.64
Elemento asta parabolico
Lo spostamento u è espresso da:
Ï u1 ¸
Ô Ô
1
1
2
u = – --- x ( 1 – x ) --- x ( 1 – x ) ( 1 – x ) Ì u 2 ý
2
2
Ôu Ô
Ó 3þ
5.65
Le derivate delle funzioni di forma rispetto a x sono:
dn 2
1
= --- + x
dx
2
dn 1
1
= – --- + x
2
dx
dn 3
= – 2x
dx
5.66
per cui la matrice jacobiana è:
[ J ] = Ê – 1--- + xˆ
Ë 2
¯
Ï x1 ¸
Ô Ô
Ê 1--- + xˆ – 2x Ì x ý
2
Ë2
¯
Ô Ô
x
Ó 3þ
5.67
funzione lineare di c . Se il nodo centrale 3 è in mezzeria dell'elemento la
matrice jacobiana è una costante e si ha:
l
[ J ] = --2
2
[ J ] –1 = --l
l
det [ J ] = --2
5.68
per cui la matrice [ b ] è:
2
[ b ] = --- Ê – 1--- + xˆ
l Ë 2
¯
158
Ê 1--- + xˆ – 2x
Ë2
¯
5.69
ELEMENTI ISOPARAMETRICI
Supponendo costanti l'area A ed il modulo elastico E , si ha:
2
Ê x – 1---ˆ
Ë
¯
2
1
x 2 – --4
1
– 2x Ê x – ---ˆ
Ë
2¯
1
x 2 – --4
2
Ê x + 1---ˆ
Ë
2¯
1
– 2x Ê x + ---ˆ
Ë
2¯
1
– 2x Ê x – ---ˆ
Ë
2¯
1
– 2x Ê x + ---ˆ
Ë
2¯
4x 2
1
2EA
[ k ] = ---------- Ú
l
–1
dx
5.70
che integrata analiticamente da:
7
--6
2EA 1
[ k ] = ---------- --l
6
4
– --3
1
4
--- – --6
3
7
4
--- – --6
3
4
8
--– --3
3
5.71
ESEMPIO 5.2
Calcolare la matrice di rigidezza per l'elemento piano a tre nodi in stato di
tensione piano o deformazione piano.
Gli spostamenti u e v sono dati da:
u = ( 1 – x – h )u 1 + xu 2 + hu 3
5.72
v = ( 1 – x – h )v 1 + xu 2 + hv 3
le derivate delle funzioni di forma rispetto alle coordinate locali sono:
∂n 1
= –1
∂x
∂n 2
= 1
∂x
∂n 3
= 0
∂x
∂n 1
= –1
∂h
∂n 2
= 0
∂h
∂n 3
= 1
∂h
5.73
e gli elementi della matrice jacobiana:
∂x
=
∂x
∂x
=
∂h
[J] =
3
∂n i
 ∂ x xi
i=1
3
Â
i=1
= x2 – x1
∂n i
x = x3 – x1
∂h i
( x2 – x1 ) ( y2 – y1 )
( x3 – x1 ) ( y3 – y1 )
∂y
=
∂x
∂y
=
∂h
3
∂n i
 ∂ x yi =
i=1
3
Â
i=1
y2 – y1
∂n i
y = y3 – y1
∂h i
5.74
5.75
159
ELEMENTI ISOPARAMETRICI
e l'inversa:
( y3 – y1 ) –( y2 – y1 )
1
[ J ] –1 = --------------det [ J ] – ( x – x ) ( x – x )
3
1
2
1
5.76
dove:
det [ J ] = ( x 2 – x 1 ) ( y 3 – y 1 ) – ( x 3 – x 1 ) ( y 2 – y 1 ) = 2A
5.77
essendo A l'area dell'elemento triangolare. Le derivate delle funzioni di forma
rispetto alle coordinate fisiche sono quindi:
Ï
Ô
Ô
Ì
Ô
Ô
Ó
∂n i ¸
Ô
∂x Ô
1 Ï yi + 1 – yi + 2 ¸
ý = -----Ì
ý
2A
∂n i Ô
Ó xi + 2 – xi + 1 þ
Ô
∂y þ
i = 1, 2, 3
5.78
e infine la matrice [ b ] :
[b] =
1
-----2A
5.79
( y2 – y3 )
0
( y3 – y1 )
0
( y1 – y2 )
0
0
( x3 – x2 )
0
( x1 – x3 )
0
( x2 – x1 )
( x3 – x2 )
( y2 – y3 )
( x1 – x3 )
( y3 – y1 )
( x2 – x1 )
( y1 – y2 )
simile a quella ottenuta nell'esempio 4.6 mediante le coordinate generalizzate, tenendo conto del fatto che qui si è ordinato il vettore degli spostamenti nodali { s } nel seguente modo:
{ s } T = { u 1 v 1 u 2 v 2 u 3 v 3}
5.80
5.3 CARICHI NODALI EQUIVALENTI
5.3.1 Carichi di linea
È qui ricavata l'espressione per il carico nodale equivalente a carichi distribuiti
lungo un contorno nel caso di elementi piani isoparametrici. Si considerano lati
aventi due o tre nodi, tipici di elementi a 3 o 4 nodi e 6 o 8 nodi rispettivamente.
La figura 5.22 illustra i due tipi di contorno; i nodi 1 e 2 nel caso lineare e 1, 3, 2
nel caso parabolico individuano il contorno su cui agiscono i carichi per unità di
lunghezza e stabiliscono il verso della coordinata naturale x sulla linea.
In entrambi i casi illustrati in figura 5.22 sono definiti un carico distribuito tan160
ELEMENTI ISOPARAMETRICI
genziale qt , positivo se diretto nel verso individuato dai nodi 1 e 2, e un carico
distribuito normale qn , positivo se agisce nel senso che si ottiene ruotando in
senso antiorario il vettore tangente alla linea.
Fig. 5.22 – Carichi di linea per elementi lineari e parabolici.
Le coordinate di un punto qualsiasi sulla linea sono definite da:
n
x =
n
 ni xi
y =
i=1
 ni yi
5.81
i=1
dove n è il numero di nodi lungo la linea e le ni sono le funzioni di forma definite che per un contorno definito da due nodi sono:
1–x
n 1 = ----------2
1+x
n 2 = -----------2
5.82
mentre per un lato definito da tre nodi:
1
n 1 = – --- x ( 1 – x )
2
1
n 2 = --- x ( 1 + x )
2
n3 = 1 – x 2
5.83
I carichi distribuiti tangenziale e normale al contorno possono essere espressi in
funzione dei valori che tali carichi assumono in corrispondenza dei nodi, in un
modo del tutto analogo alla definizione della geometria o della cinematica:
n
qt =
 ni qt
i=1
n
i
qn =
 ni qn
i
5.84
i=1
Si nota che per lati definiti da due nodi è possibile considerare al massimo carichi
distribuiti variabili linearmente, mentre per contorni definiti da tre nodi si potrà
avere al massimo una legge parabolica.
Il vettore del carico nodale equivalente è:
161
ELEMENTI ISOPARAMETRICI
{ fe }q =
Úl [ n ]T { q } dl = Úl [ n ]T { dq }
5.85
dove {dq}, vettore delle forze elementari agenti secondo gli assi coordinati x e y,
è dato da:
{ dq } T = { dq x
dq y }
5.86
Le componenti dqx, dqy possono essere ricavate considerando un tratto dl di
linea (fig. 5.23):
dq x = qt dl cos a – qn dl sin a
dq y = qt dl sin a + qn dl cos a
5.87
Fig. 5.23 – Componenti del carico di linea.
Osservando che:
dl cos a = dx
dl sin a = dy
5.88
si ottiene:
dq x = qt dx – qn dy
5.89
dq y = qt dy + qn dx
I termini dx, dy in funzione della coordinata naturale x sono:
dx =
per cui si ha:
162
∂x
dx
∂x
dy =
∂y
dx
∂x
5.90
ELEMENTI ISOPARAMETRICI
∂x
∂y
dq x = Ê q t – q n ˆ d x
Ë ∂x
∂x¯
5.91
∂y
∂x
dq y = Ê q t + q n ˆ d x
Ë ∂x
∂x¯
Le componenti f x , f y del carico nodale equivalente al nodo i (i = 1,2 o
i
i
i = 1,2,3) sono quindi calcolate come:
1
f xi =
f yi =
∂x
- – qn
Ú–1 ni ÊË qt -------∂x
∂y
------ ˆ d x
∂x ¯
5.92
∂y
∂x
n i Ê q t ------ + q n ------ ˆ d x
Ë ∂x
∂x ¯
–1
Ú
1
dove l'integrale, esteso alla lunghezza del contorno, è generalmente valutato
numericamente.
Nel caso di un contorno definito da due nodi è possibile determinare in modo
abbastanza semplice le componenti del vettore del carico nodale equivalente in
forma esplicita; a titolo di esempio si determinerà la componente f x :
1
f x1 =
∂x
∂y
n 1 Ê q t ------ – q n ------ ˆ d x
Ë ∂x
∂x ¯
–1
Ú
1
5.93
dove:
1
n 1 = --- ( 1 – x )
2
5.94
1
1
q t = --- ( 1 – x )q t1 + --- ( 1 + x )q t2
2
2
1
1
q n = --- ( 1 – x )q n1 + --- ( 1 + x )q n2
2
2
5.95
∂n 2
∂n 1
x2 – x1
∂x
=
x1 +
x 2 = --------------∂x
∂x
2
∂x
∂n 1
∂n 2
y2 – y1
∂y
=
y1 +
y 2 = -------------∂x
∂x
∂x
2
5.96
Sostituendo ed integrando per x da –1 a +1, si ha:
1
1
1
1
f x1 = ( x 2 – x 1 ) Ê --- q t1 + --- q t2ˆ + ( y 2 – y 1 ) Ê --- q n1 + --- q n2ˆ
Ë3
Ë3
6 ¯
6 ¯
5.97
163
ELEMENTI ISOPARAMETRICI
5.3.2 Carichi di volume
Il carico nodale equivalente a carichi di volume è dato da:
{ fe }f =
ÚV [ n ]T { f } dV
5.98
dove {} nel caso piano è:
{ f } T = { fx
fy }
5.99
e dV = h det [ J ] d x d . Inoltre nel caso di peso proprio con l'asse Y orientato
nel senso della gravità si ha:
fx = 0
f y = rg
5.100
Il calcolo del vettore dei carichi nodali equivalenti segue le stesse procedure viste
per il calcolo delle matrici di rigidezza; ad esempio nel caso piano è calcolabile
come:
m
{ fe }f =
m
  wi wj [ n ]ijT { f }det [ J ]ij h
5.101
i = 1j = 1
5.3.3 Effetti termici
Il carico nodale equivalente ad effetti termici viene calcolato come:
{ fe }e =
0
ÚV [ b ]T [ E ] { e0 } d V
5.102
dove { e 0 } nel caso di stato di tensione piano è:
{ e 0 } T = a* ( T – T 0 ) { 1
1
0}
5.103
e, nel caso di stato di deformazione piano:
{ e 0 } T = a* ( T – T 0 ) { 1 + n
1+n
0}
5.104
in cui a* è il coefficiente di dilatazione termica del materiale, T 0 è la temperatura di riferimento e T è la temperatura in un punto generico, che può essere
espressa come:
n
T =
 ni Ti
5.105
i=1
dove n è il numero di nodi dell'elemento, n i sono le funzioni di forma associate
all'elemento e T i è il valore della temperatura al nodo i dell'elemento.
164
ELEMENTI ISOPARAMETRICI
Si ha quindi:
m
{ fe }e =
0
m
  wi wj
[ b ] ijT [ E ] { e 0 } ij det [ J ] ij h
5.106
i = 1j = 1
5.4 CALCOLO DELLE TENSIONI
Noto il vettore degli spostamenti { S } della struttura gli spostamenti { s } h
relativi ai singoli elementi possono essere ricavati con una operazione inversa a
quella utilizzata per assemblare gli elementi.
Le deformazioni nell'elemento sono quindi date da:
{e} = [b]{s}
5.107
e le tensioni:
{ s } = [ E ] ( [ b ] { s } – { e0 } ) + { s0 }
5.108
Lo stato di tensione dipende dal punto che si considera all'interno dell'elemento
perché in generale la matrice [ b ] è funzione delle coordinate; questa deve
quindi essere calcolata nel punto dove si vuole conoscere lo stato di tensione.
I punti nei quali si desidererebbe conoscere lo stato di tensione sono i nodi
dell’elemento, perché è nota la loro posizione geometrica nella struttura e quindi
è abbastanza semplice associare uno stato di tensione ad un punto della struttura.
Sfortunatamente i nodi sono i punti nei quali lo stato di tensione risulta meno
preciso, anzi sono i punti nei quali lo scarto tra valore calcolato e valore esatto è
il più alto all'interno dell'elemento.
Generalmente si ha un errore minore se lo stato di tensione viene calcolato in
determinati punti interni dell'elemento; questi punti possono essere utilizzati per
caratterizzare il campo di tensione che poi potrà essere estrapolato per valutare le
tensioni nei nodi dell'elemento.
La determinazione dei punti ottimali per il calcolo delle tensioni può essere fatta
secondo quanto suggerito da Barlow1; l’obiettivo è quello di determinare i punti
nei quali le tensioni sono calcolate con lo stesso grado di precisione con il quale
sono calcolati gli spostamenti nodali.
Si assume per il dato elemento un campo di spostamenti { u' } di un ordine superiore a quello comunemente utilizzato { u } e si cercano i punti nei quali lo stato
di tensione calcolato a partire dai due campi di spostamento è uguale
{ s ' } = { s } . In questi punti lo stato di tensione { s } avrà lo stesso grado di
1.Barlow, J., Optimal stress location in finite element models, International Journal for
Numerical Methods In Engineering, 1976, vol. 10, p. 243-251.
165
ELEMENTI ISOPARAMETRICI
precisione dello spostamento { s } ; detti punti coincidono con i punti di Gauss
di sottointegrazione e sono talvolta chiamati punti di Barlow.
Nel caso monodimensionale della trave inflessa il procedimento seguito da Barlow per la determinazione dei punti ottimali può essere cos riassunto. Il campo
di spostamenti può essere espresso come:
va = 1 x x 2 x 3 { a }
5.109
dove x è una coordinata adimensionale misurata dal centro della trave e variabile
da –1 a 1 e {a} è il vettore delle coordinate generalizzate.
Gli spostamenti nodali { s a } sono dati da:
{ sa } = [ A ] { a }
5.110
Assumendo una funzione di ordine superiore, lo spostamento della trave sarà
dato da:
vb = 1 x x 2 x 3 x 4 { b }
5.111
e gli spostamenti nodali { s b } sono dati da:
{ sb } = [ B ] { b }
5.112
Se l’elemento è capace di rappresentare correttamente il campo di spostamenti
nei due casi considerati, allora gli spostamenti nodali ottenuti nei due casi
saranno approssimativamente uguali:
[A]{a} ª [B]{b}
5.113
{ a } = [ A ] –1 [ B ] { b }
5.114
con:
[ A ] –1 [ B ] =
1
0
0
0
0
1
0
0
0
0
1
0
0
0
0
1
–1
0
2
0
5.115
Le deformazioni flessionali e quindi le tensioni flessionali saranno uguali nei due
casi quando:
2
2
∂ va
∂x2
166
=
∂ vb
∂x2
5.116
ELEMENTI ISOPARAMETRICI
ovvero quando:
[ 0 0 2 6 x ] { a } = 0 0 2 6 x 12 x2 { b }
5.117
e sostituendo in termini di {b}:
0 0 2 6 x 4 { b } = 0 0 2 6x 12x 2 { b }
5.118
Questa uguaglianza è soddisfatta per tutti i valori di { b} solo se:
1
x = ± ------3
5.119
La stessa tecnica può essere applicata all'elemento piano isoparametrico a quattro
nodi; si trova che i punti ottimali per il calcolo delle tensioni sono:
x = 0
h = 0
5.120
mentre per l'elemento piano isoparametrico a otto nodi si ha:
1
x = ± ------3
1
h = ± ------3
5.121
Le tensioni ai nodi dell'elemento possono quindi essere ottenute estrapolando i
valori calcolati nei punti di Gauss (altrimenti detti di Barlow). Nel caso monodimensionale (trave inflessa o elemento asta a tre nodi) si ottiene:
1+ 3
---------------Ï s1 ¸
2
Ì ý =
Ó s2 þ
1– 3
---------------2
1– 3
---------------- Ï s ¸
I
2
Ì
ý
1 + 3 Ó s II þ
---------------2
5.122
mentre nel caso bidimensionale (elemento piano isoparametrico a 8 nodi) si
ottiene:
3
1 + ------2
Ï s1 ¸
1
Ô Ô
– --Ô s2 Ô
2
Ì ý =
s
Ô 3Ô
3
1 – ------Ô Ô
2
Ó s4 þ
1
– --2
1
– --2
3
1 – ------2
1
– --2
3
1 + ------2
1
– --2
3
1 – ------2
1
– --2
3
1 + ------2
1
– --2
3
1 – ------2
1
– --2
3
1 + ------2
Ï sI ¸
Ô
Ô
Ô s II Ô
Ì
ý
Ô s III Ô
Ô
Ô
Ó s IV þ
5.123
e nel caso tridimensionale (elemento solido isoparametrico a 20 nodi) si ottiene:
167
ELEMENTI ISOPARAMETRICI
Ï s1 ¸
a3
Ô
Ô
Ô s2 Ô
a2b
Ô
Ô
Ô s3 Ô
ab 2
Ô
Ô
Ô s4 Ô
a2b
Ì
ý =
Ô s5 Ô
a2b
Ô
Ô
ab 2
Ô s6 Ô
Ô
Ô
b3
Ô s7 Ô
Ô
Ô
ab 2
Ó s8 þ
a = 1+ 3
a2b
ab 2
a2b
a2b
ab 2
b3
a3
a2b
ab 2
ab 2
a2b
ab 2
a2b
a3
a2b
b3
ab 2
a2b
ab 2
a2b
a3
ab 2
b3
ab 2
ab 2
b3
ab 2
a3
a2b
ab 2
a2b
ab 2
b3
a2b
a3
a2b
ab 2
a2b
ab 2
ab 2
a2b
a3
b3
ab 2
a2b
a2b
ab 2
a2b
Ï s ¸
ab 2 Ô I Ô
Ô s Ô
b 3 Ô II Ô
ab 2 ÔÔ s III ÔÔ
a 2 b Ô s IV Ô
Ì
ý
a 2 b Ô sV Ô
Ô
Ô
ab 2 Ô s VI Ô
Ô
Ô
a 2 b Ô s VII Ô
Ô
a 3 ÔÓ s
VIII þ
5.124
b = 1– 3
Nella tabella 5.11 sono indicati i punti ottimali per alcuni elementi.
Tab. 5.11 – No di punti e coordinate per il calcolo ottimale delle tensioni
ELEMENTO
TRAVE INFLESSA
NO PUNTI
COORDINATE NATURALI
2
x = ±a
ISOPARAMETRICO PIANO A
4 NODI
1
x = 0;
ISOPARAMETRICO PIANO A
8 NODI
2¥2
x = ± a;
x = 0;
ISOPARAMETRICO SOLIDO A
8 NODI
1
ISOPARAMETRICO SOLIDO A
20 NODI
2¥2¥2
h = 0
h = ±a
h = 0
x = ± a; h = ± a; z = ± a
a = ±1 § 3
Le tensioni così calcolate ai nodi di un elemento saranno in generale diverse dalle
tensioni calcolate negli stessi nodi appartenenti però ad elementi adiacenti; ciò
risulta dal fatto che le tensioni ai nodi dipendono dagli spostamenti di tutti i
nodi dell'elemento.
Ciò deriva anche dal fatto che nell’equazione dei lavori virtuali è garantito l’equilibrio globale (equilibrio delle forze nodali), ma non quello locale punto per
punto (equilibrio delle tensioni).
Ad un nodo allora sono attribuiti diversi valori di tensione, ciascuno derivante
dai vari elementi che concorrono in quel nodo; è usuale calcolare, a partire da dati
valori, una tensione media nodale ed attribuire al nodo quest'unico valore di tensione media, valore che sarà probabilmente più esatto dei singoli valori calcolabili.
168
ELEMENTI ISOPARAMETRICI
5.5 PROBLEMI RELATIVI AGLI ELEMENTI ISOPARAMETRICI
5.5.1 Scelta dell'ordine di integrazione
La scelta dell'ordine di integrazione m dipende dal grado della funzione integranda; nel caso della matrice di rigidezza di un elemento piano l'ordine di integrazione si può valutare determinando il grado massimo della funzione:
f = [ b ] T [ E ] [ b ]det [ J ]h
5.125
tenendo conto che un polinomio di grado 2m – 1 è integrato esattamente da m
punti di Gauss.
Per elementi di forma regolare è abbastanza semplice determinare il grado della
funzione integranda, perché det [ J ] è costante; non altrettanto si può dire per
elementi di forma complessa, dove il determinante della matrice jacobiana varia
in modo irregolare. Di conseguenza l’ordine di integrazione viene anche stabilito
in base a verifiche numeriche sull'elemento.
Una integrazione esatta della matrice di rigidezza, in special modo per elementi
di grado elevato, può portare ad avere lunghi tempi di calcolo, per cui è a volte
preferibile una integrazione ridotta o sottointegrazione; si utilizzano cioè meno
punti di campionamento di quelli necessari per integrare esattamente la funzione
integranda.
A tal proposito si ricorda ora che con la formulazione degli elementi finiti a spostamenti assegnati si ha sempre una stima in eccesso della rigidezza del sistema
(avendo limitato la cinematica del sistema a quella espressa dalle funzioni di
forma). Quindi sottointegrare, cioè sottovalutare la rigidezza dell'elemento, può,
di fatto, portare a risultati migliori a patto che l'errore nell'integrazione numerica compensi la sovrastima della rigidezza strutturale dovuta al tipo di formulazione.
Esiste però un limite alla sottointegrazione perché è necessario poter calcolare
esattamente il volume dell’elemento. Questa condizione deriva dal fatto che al
limite, per suddivisioni molto fitte, nell'elemento si verifica uno stato di deformazione costante (cosa che l'elemento deve in ogni caso essere in grado di rappresentare) e quindi l'energia di deformazione può, in tali condizioni assumere la
forma, limitandosi al caso piano:
U µ [ k* ]
ÚVh det [ J ] d x d h
5.126
Quindi si deve integrare esattamente il prodotto det [ J ] h; nel caso di spessore
costante e di elemento lineare a quattro nodi la funzione integranda è una funzione lineare in x e e quindi è necessario un solo punto di campionamento; nel
caso di elemento a otto nodi a spessore costante si ha invece una dipendenza
cubica da x e e quindi è necessario uno schema di integrazione 2 x 2.
169
ELEMENTI ISOPARAMETRICI
In pratica però una sottointegrazione troppo spinta può portare a matrici singolari quando in tutti i punti di campionamento l'energia di deformazione è nulla.
Si hanno in tal caso nell'elemento modi spuri o meccanismi nascosti, ovvero modi
di deformazione con energia di deformazione nulla.
Per un singolo elemento il numero di meccanismi può essere calcolato mediante
la seguente relazione:
M = GN – L – rn
5.127
dove G è il numero di gradi di libertà per nodo, N il numero di nodi dell'elemento, L il numero di gradi di libertà di moto rigido, r il rango della matrice di
elasticità [E ] e n il numero di punti di campionamento. Se M non si hanno
meccanismi.
Si può pensare infatti di intendere l'integrazione numerica nel seguente modo:
Integrare numericamente significa campionare la funzione integranda in alcuni
punti nel dominio dell'elemento. Nel nostro caso significa valutare la rigidezza
dell'elemento in alcuni punti ed attribuire questa alla rigidezza reale.
Ciò equivale, fisicamente, a sostituire una rigidezza distribuita sull'elemento con
delle molle spia concentrate nei punti di integrazione.
Il numero delle molle da considerare per ogni punto dipende dal numero delle
componenti del vettore deformazione (o tensione) e cioè dal rango della matrice
[E ] (ad esempio nel caso piano si può considerare una rigidezza associata alla
componente e xx , una alla e yy ed una terza associata alla g xy .
Affinché la matrice di rigidezza [k], depurata dai moti rigidi, non sia singolare, è
necessario che il numero di molle spia sia maggiore del numero di gradi di
libertà di deformazione dell'elemento, cioè:
rn ≥ GN – L
5.128
Ad esempio nel caso di problemi bidimensionali in stato di tensione piano o
deformazione piano il numero di gradi di libertà per nodo è G = 2, il rango della
matrice [E] è r = 3, ed il numero di gradi di libertà di moto rigido è L = 3 per cui
si ha:
M = 2N – 3 – 3n
5.129
Per un elemento a 4 nodi, bilineare, sottointegrato con un punto di campionamento (n = 1) si hanno due meccanismi:
M = 2◊4–3–3◊1 = 2
5.130
mentre nel caso di un elemento parabolico a otto nodi (Serendipity) con una sottointegrazione con quattro punti di campionamento (2 x 2) si ha un meccanismo nascosto.
In quest’ultimo caso però se si hanno più elementi adiacenti il meccanismo non
170
ELEMENTI ISOPARAMETRICI
si propaga perché gli spostamenti corrispondenti ad energia nulla in un elemento
sono incompatibili con gli spostamenti di elementi adiacenti per cui un insieme
di elementi può non avere modi di deformazione ad energia nulla, anche se questi sono presenti in un elemento, così che una integrazione 2 x 2 può essere utilizzata efficacemente.
La tabella 5.12 riassume gli ordini di integrazione per elementi piani, distinguendo tra ordini di integrazione massimi e ordini di integrazione ridotti.
Tab. 5.12 – Ordine di integrazione per elementi quadrangolari
ELEMENTO
ORDINE DI INTEGRAZIONE
MASSIMO
ORDINE DI INTEGRAZIONE
RIDOTTO
4 NODI RETTANGOLO
2¥2
2¥2
4 NODI DISTORTO
2¥2
2¥2
8 NODI RETTANGOLO
3¥3
2¥2
8 NODI DISTORTO
3¥3
3¥3
Tuttavia nel caso di sottointegrazione è consigliabile effettuare verifiche numeriche dell'elemento, quali il calcolo degli autovalori della matrice di rigidezza ed il
patch test.
5.5.2 Distorsione degli elementi
Gli elementi forniscono i risultati migliori quando la loro forma geometrica
coincide con quella dell'elemento di riferimento, come triangoli equilateri e quadrati nel piano, o cubi nello spazio.
In pratica comunque è quasi impossibile che tutti gli elementi abbiano una
forma geometrica ideale; si parla allora di elementi distorti (rispetto alla forma di
riferimento). Un limitato valore di distorsione può essere comunque accettato.
Oltre un determinato limite di distorsione, però, la rigidezza dell'elemento viene
valutata con scarsa precisione; oltre un limite critico diventa addirittura impossibile calcolare la matrice di rigidezza (la matrice jacobiana diventa singolare).
È difficile determinare quello che costituisce il limite accettabile di distorsione;
oltretutto il comportamento è diverso da elemento a elemento.
I tipi di distorsione di un elemento possono essere classificati in tre categorie
fondamentali: il rapporto lunghezza altezza, il rapporto di forma e la distorsione
171
ELEMENTI ISOPARAMETRICI
angolare (fig. 5.24).
Non esiste un modo generale per analizzare questi tipi di distorsioni; è necessario
verificare numericamente l'elemento in diverse condizioni di distorsione e di
carichi.
Fig. 5.24 – Distorsione di un elemento.
Con riguardo al primo tipo di distorsione (fig. 5.24a) si può accettare un rapporto 1:2 per elementi situati in zone dove si vuole conoscere lo stato di tensione
e di 1:10 per elementi lontani da zone critiche.
Nel caso della figura 5.24b un rapporto accettabile è di 1:4 in zone critiche e di
1:10 in zone non critiche.
Per ciò che riguarda la distorsione angolare (fig. 5.24c) è in genere consigliabile
non avere angoli interni minori di 45˚ per elementi quadrangolari e minori di
15˚ per elementi triangolari.
Una verifica della distorsione può essere fatta anche mediante il calcolo della
matrice jacobiana e del suo determinante. Il determinante della matrice jacobiana è infatti costante se l’elemento è un parallelepipedo con i nodi sui lati
equamente spaziati, ma varia se esiste altro tipo di distorsione. Esso può al limite
diventare nullo o negativo nelle vicinanze di una zona ad eccessiva distorsione,
con la conseguenza di impedire il calcolo della matrice di rigidezza dell'elemento. Causa di ciò è la mancanza di corrispondenza biunivoca con conseguente
non invertibilità della matrice jacobiana.
Ad esempio nel caso dell’elemento asta parabolico a tre nodi la matrice jacobiana
è data da:
[ J ] = Ê – 1--- + xˆ
Ë 2
¯
172
Ï x1 ¸
Ô Ô
1
Ê --- + xˆ ( – 2x ) Ì x ý
2
Ë2
¯
Ô Ô
Ó x3 þ
5.131
ELEMENTI ISOPARAMETRICI
Fig. 5.25 – Distorsione di un elemento asta a 3 nodi.
Se il nodo centrale non è in mezzeria dell’elemento, supponendo x1 = 0, x2 = l
con x3 compresa tra 0 e l (fig. 5.25) il determinante dello jacobiano è:
1
J = Ê --- + xˆ l – 2xx 3
Ë2
¯
5.132
La posizione x3 del nodo centrale che rende singolare la matrice jacobiana può
essere calcolata annullando il determinante dello jacobiano e tenendo presente
che – 1 £ x £ 1 . Per x = – 1 si ha x3 = l/4 mentre per x = 1 si ha x3 = 3l/4.
I valori ammissibili di x3 per i quali [ J ] non è singolare sono quindi:
3l
l
--- £ x 3 £ ---4
4
5.133
Fig. 5.26 – Distorsione di un elemento quadrangolare.
Nel caso di elementi quadrangolari a quattro e otto nodi Zienkiewicz ha proposto, per non avere una matrice jacobiana singolare, che gli angoli interni siano
minori di 180 gradi e per gli elementi parabolici i nodi sui lati debbano essere
compresi nella terza parte centrale (fig. 5.26).
Cook ha proposto il calcolo di un fattore di distorsione DP, dato rispettivamente
per elementi quadrangolari e elementi solidi:
173
ELEMENTI ISOPARAMETRICI
minJ
DP = 4 Ê -----------ˆ
Ë Area ¯
minJ
DP = 8 Ê ------------------ˆ
Ë Volume¯
5.134
parallelepipedo; se DP £ 0.2 si ha una distorsione eccessiva.
È da notare che nei casi della figura 5.24 il fattore di distorsione è sempre 1, cioè
il determinante della matrice jacobiana è costante e quindi non può essere preso
in considerazione quale misura di quel tipo di distorsione.
5.5.3 Integrazione selettiva e modi incompatibili
L’elemento isoparametrico lineare (piano a 4 nodi e solido a 8 nodi) non è in
grado di rappresentare correttamente uno stato di tensione con presenza di gradienti.
Ad esempio nel caso di sollecitazione a flessione pura l’elemento presenta un
comportamento molto più rigido della realtà e questo può essere imputato
all'effetto del taglio, che teoricamente deve essere nullo.
Nel caso di integrazione esatta (n = 2 x 2) infatti si ha, a causa dei vincoli cinematici imposti dalle funzioni di forma adottate, una deformazione a taglio nei
punti di campionamento utilizzati (fig. 5.27a).
Se invece si procede con una sottointegrazione (n = 1) si nota che nel punto di
campionamento (punto centrale) non si ha deformazione ha taglio e quindi
l'energia di deformazione a taglio è calcolata in modo esatto (fig. 5.27b).
Fig. 5.27 – Punti di integrazione per elemento a 4 nodi.
Il campo di spostamenti all'interno dell’elemento è definito da:
u = u* xh
ed essendo:
174
v = 0
5.135
ELEMENTI ISOPARAMETRICI
2x
x = ----l
2y
h = ----h
5.136
si ha:
4
u = ----u* xy
lh
v = 0
5.137
Il campo corretto di spostamenti invece prevede:
4
u = ----u* xy
lh
lu*
hu* Ê
4x 2ˆ
4y 2
--------- 1 – --------ˆ
v = ------- Ê 1 – ------+
n
2h Ë
2l Ë
l2 ¯
h2 ¯
5.138
Nel caso di sollecitazione derivante da momento flettente puro si deve avere che
la deformazione a taglio xy sia nulla dovunque; ciò è verificato dal campo di
spostamenti corretto:
g xy =
4
4
∂u ∂v
+
= ----u* x – ----u* x = 0
lh
lh
∂y ∂x
5.139
ma non da quello associato all'elemento quadrangolare:
g xy =
4
∂u ∂v
+
= ----u* x
lh
∂y ∂x
5.140
La stessa cosa si ottiene considerando un momento flettente applicato lungo i
lati di lunghezza l; il campo di spostamenti approssimato risulterà:
4
v = ----v* xy
lh
5.141
Si può notare che i due modi flessionali forniscono uno stato di deformazione
esatto solo nel punto centrale dell'elemento x = 0; = 0; questo suggerisce
quindi la possibilità di ridurre la sovrastima di rigidezza effettuando una integrazione selettiva in cui le deformazioni normali sono integrate esattamente
mediante una griglia 2 x 2, mentre i termini corrispondenti alla deformazione
dovuta al taglio sono sottointegrati, calcolati cioè nel punto centrale dell'elemento.
Un modo per ottenere risultati più soddisfacenti è quello di ricorrere ad elementi
quadratici (piano a 8 nodi e solido a 20 nodi); questo modo di procedere se può
essere accettabile nel caso bidimensionale (si passa da 8 a 16 gradi di libertà per
elemento) è però molto gravoso nel caso tridimensionale, dove si passa da 24 a
60 gradi di libertà, quasi triplicando l'ordine delle matrici in esame.
È quindi di grande importanza riuscire a migliorare il comportamento flessionale degli elementi lineari; tra i diversi procedimenti proposti si hanno il metodo
dell'integrazione selettiva o quello dell'aggiunta di modi incompatibili.
175
ELEMENTI ISOPARAMETRICI
Integrazione selettiva
Nell’integrazione selettiva si utilizzano ordini di integrazione diversi per le
diverse forme di energia di deformazione (flessionale e di taglio). Si ha infatti, da:
[k] =
ÚA [ b ]T [ E ] [ b ]h dA
5.142
separando la deformazione dovuta al taglio dalla deformazione normale:
[b] =
[E] =
[ b1 ]
5.143
[ b2 ]
[ E 11 ] [ E 12 ]
5.144
[ E 21 ] [ E 22 ]
essendo, ad esempio nel caso di stato di tensione piano:
E
[ E 11 ] = --------------2 1 n
1–n n 1
[ E 22 ] = [ G ]
[ E 12 ] = [ E 21 ] T = 0
0
5.145
La matrice di rigidezza dell’elemento può allora essere scritta come somma del
contributo dovuto alle deformazioni normali e di quello dovuto alle deformazioni di taglio:
[k] =
=
ÚV { e }T { s } dV + ÚV { g }T { t } dV=
ÚA [ b1 ] [ E11 ] [ b1 ]h dA + ÚA [ b2 ] [ E22 ] [ b2 ]h dA
5.146
Modi incompatibili
Un altro modo per migliorare la risposta dell’elemento lineare è quello di correggere le funzioni di forma considerando i termini flessionali mancanti (1 – x2),
(1 – 2).
Questo tipo di formulazione è stato presentato inizialmente da Wilson (1973)
(elemento Q6 ); successivamente (Taylor, 1976) propose uno schema di integrazione diverso per l'elemento Q6, formulando quindi l’elemento QM6. La differenza tra i due elementi è che il secondo supera il patch test per qualsiasi
configurazione geometrica, mentre il primo solo se di forma rettangolare. Nel
caso bidimensionale si ha:
176
ELEMENTI ISOPARAMETRICI
4
u =
 ni (x ,h )ui + ( 1 – x 2 )a1 + ( 1 – h 2 )a2
i=1
4
v =
5.147
 ni (x ,h )vi + ( 1 – x 2 )a3 + ( 1 – h 2 )a4
i=1
mentre nel caso tridimensionale:
8
u =
 ni (x ,h ,z)ui + ( 1 – x 2 )a1 + ( 1 – h 2 )a2 + ( 1 – z 2 )a3
i=1
8
v =
 ni (x ,h ,z)vi + ( 1 – x 2 )a4 + ( 1 – h 2 )a5 + ( 1 – z 2 )a6
5.148
i=1
8
w =
 ni (x ,h ,z)wi + ( 1 – x 2 )a7 + ( 1 – h 2 )a8 + ( 1 – z 2 )a9
i=1
dove i coefficienti ai sono definiti come spostamenti generalizzati e possono
essere pensati come gradi di libertà interni dell'elemento; è da notare che un elemento così definito è incompatibile essendo gli ai definiti univocamente per ciascun elemento. In generale si ha:
Ï {s} ¸
{ u } = [ ns ] [ na ] Ì
ý
Ó {a} þ
5.149
Ï {s} ¸
[ b ] = [ bs ] [ ba ] Ì
ý
Ó{ a } þ
5.150
e:
Mentre per l’elemento lineare compatibile la matrice [b] ha dimensioni 3 x 8,
per l’elemento incompatibile si hanno 4 colonne in più, dovute ai termini ai:
177
ELEMENTI ISOPARAMETRICI
ºº
∂n i
∂n 1 ∂n 2
0 ºº
0
∂x
∂x ∂x
[ b ] = º º 0 ∂n i º º 0
∂y
ºº
0
∂n 1 ∂n 2
0
∂y ∂y
∂n i ∂n i
∂n 1 ∂n 2 ∂n 1 ∂n 2
ºº
∂y ∂x
∂y ∂y ∂x ∂x
Ï
Ô
Ô
Ô
Ô
Ô
Ô
Ô
Ì
Ô
Ô
Ô
Ô
Ô
Ô
Ô
Ó
º
º
ui
vi
º
º
a1
a2
a3
a4
¸
Ô
Ô
Ô
Ô
Ô
Ô
Ô
ý
Ô
Ô
Ô
Ô
Ô
Ô
Ô
þ
5.151
La matrice di rigidezza avrà ordine 12 x 12 e può essere scritta separando le varie
sottomatrici, come:
[k] =
[ k ss ] [ k sa ]
[ k as ] [ k aa ]
5.152
Poche i gradi di libertà ai sono interni a ciascun elemento essi possono essere eliminati o condensati a livello elemento; questo procedimento, noto come condensazione statica, porta ad ottenere una matrice di rigidezza condensata di
dimensioni 8 x 8. Dalle equazioni elemento si ha:
[ k ss ] [ k sa ] Ï { s } ¸
Ï{ f }¸
=
Ì
ý
Ì
ý
[ k as ] [ k aa ] Ó{ a }þ
Ó{ 0 }þ
5.153
non essendoci forze corrispondenti ai gradi di libertà aggiuntivi. Dal secondo
sottosistema si ricava:
[ k as ] { s } + [ k aa ] { a } = { 0 }
5.154
da cui, essendo [kaa ] non singolare:
{ a } = – [ k aa ] –1 [ k as ] { s }
5.155
dal primo sottosistema si può quindi ottenere la matrice di rigidezza condensata
dell'elemento incompatibile:
[ k ] = [ k ss ] – [ k sa ] [ k aa ] –1 [ k as ]
178
5.156
ELEMENTI ISOPARAMETRICI
L’elemento così definito è l'elemento incompatibile a 4 nodi Q6 introdotto da
Wilson2. I risultati ottenuti con questo elemento sono nettamente migliori di
quelli ottenuti con il semplice elemento a 4 nodi; fu tuttavia notato che nel caso
di elementi distorti (non rettangolari) i risultati erano insoddisfacenti. Fu anche
notato che questo elemento superava il patch test nel caso di geometria rettangolare, ma non nel caso di geometria qualsiasi.
Successivamente Taylor3 modificò la formulazione dell'elemento in modo che il
patch test fosse superato per qualsiasi configurazione geometrica; questo risultato
fu ottenuto semplicemente calcolando il contributo dei modi incompatibili
all'integrale dell'operatore [B] nel centro dell'elemento.
L’elemento così formulato, detto QM6, supera il patch test ed è quello generalmente utilizzato, in quanto fornisce risultati migliori anche dell'elemento con
integrazione selettiva.
L’idea di Taylor fu la seguente: affinché l’elemento passi il patch test (cioè si comporti come un elemento compatibile) il termini vettore {a } deve essere nullo
quando { s } corrisponde ad un moto rigido o ad uno stato di deformazione
costante { s } = { s 0 } . Dovrà cioè essere:
{ a } = – [ k aa ] –1 [ k as ] { s 0 } = 0
5.157
ed essendo [kaa ] definita positiva:
[ k as ] { s } = 0
5.158
ÚV [ ba ]T [ E ] [ bs ] { s0 } dV = ÚV [ ba ]T { s0 } dV = 0
5.159
cioè:
dove { s 0 } è il vettore delle tensioni calcolato in base ai soli spostamenti { s 0 }
e, rappresentando questi un moto rigido, è una costante. Condizione necessaria
per il superamento del patch test è quindi:
ÚV [ ba ] dV = 0
5.160
2. Wilson E. L., Taylor R. L., Doherty W. P., Ghaboussi J., Incompatible displacement
models, in Numerical and Computer Method in Structural Mechanics, 1973,
Academic Press, New York, p. 43.
3. Taylor R. L., Beresford P. J., Wilson E. L., A non conforming element for stress analysis,
International Journal for Numerical Methods in Engineering, 1976, vol 10, p. 12111219.
179
ELEMENTI ISOPARAMETRICI
In generale questa condizione non è soddisfatta, a meno che l’elemento non sia
un parallelogramma. Si può forzare questa condizione valutando numericamente
l’integrale nel centro dell'elemento (x = 0; = 0).
5.5.4 Autovalori della matrice di rigidezza
Gli autovalori della matrice di rigidezza sono proporzionali all’energia potenziale
elastica corrispondente ai vari modi di deformazione e quindi possono rappresentare una spia del numero di moti rigidi associati all'elemento.
L’energia di deformazione di un elemento può essere espressa come:
1
U = --- { s } T [ k ] { s }
2
5.161
Si può pensare di scrivere le equazioni di rigidezza riferendole ad un sistema definito dagli autovettori della matrice [k] antiche allo spazio fisico degli spostamenti {s}. Nel sistema:
[k]{s} = { f }
5.162
il vettore degli spostamenti { s } può essere espresso come combinazione lineare
degli autovettori della matrice di rigidezza [k]:
m
{s} =
 ai { f }i
= [F]{a}
5.163
i=1
dove m è il numero di gradi di libertà dell'elemento, {} è il vettore delle coordinate modali o dei fattori di partecipazione modali e [] è la matrice le cui
colonne rappresentano gli autovettori della matrice [k]:
[k]{f} = l{f}
5.164
Se gli autovettori {} sono normalizzati in modo che:
{ f } iT { f } i = 1
5.165
dalla equazione precedente, scritta per l'autovalore i-esimo, si ha:
{ f } iT [ k ] { f } i = l i { f } iT { f } i
5.166
{ f } iT [ k ] { f } i = l i
5.167
e inoltre essendo gli autovettori ortogonali:
[ F ]T[ k ][ F ] = [ L ]
dove [ L ] è una matrice diagonale contenente gli autovalori l i .
180
5.168
ELEMENTI ISOPARAMETRICI
Il sistema risolutivo può allora essere scritto come:
[k][F]{a} = { f }
5.169
e premoltiplicando per []T ambo i membri:
[ F ]T[ k ][ F ]{ a } = [ F ]T{ f }
5.170
[ L ]{ a } = [ F ]T{ f }
5.171
cioè:
ed essendo la matrice [
] diagonale il sistema può essere scomposto in m equazioni disaccoppiate:
l i a i = { f } iT { f }
5.172
dalle quali si calcolano direttamente gli m fattori di partecipazione modali:
{ f } iT { f }
a i = ----------------------li
5.173
e quindi gli spostamenti nodali {s}.
Gli autovettori { f } i possono quindi essere identificati come i modi di deformazione principali dell'elemento e gli autovalori l i come rigidezze modali relative
a tali modi di deformazione. Essi sono inoltre proporzionali all'energia potenziale elastica U i relativa all'i-esimo modo di deformazione:
{ f } iT [ k ] { f } i = 2U i = l i
5.174
Si ricorda che il sistema di equazioni di rigidezza comprende le equazioni di
equilibrio del sistema, equazioni a cui corrispondono i gradi di libertà di moto
rigido. In queste condizioni l'energia potenziale elastica assorbita dal sistema è
nulla e quindi il corrispondente autovalore sarà nullo e l'autovettore ad esso associato rappresenterà uno dei moti rigidi del sistema.
Se il numero di autovalori nulli è minore del numero di gradi di libertà di moto
rigido del sistema allora il campo di spostamenti assunto non è adeguato alla
rappresentazione della cinematica del sistema.
Se invece il numero di autovalori nulli è maggiore del numero di gradi di libertà
di moto rigido del sistema allora sono presenti dei moti rigidi spuri o meccanismi nascosti, cioè modi di deformazione associati ad una energia di deformazione nulla.
Se inoltre gli autovalori non variano al variare dell'orientamento dell'elemento
rispetto al sistema di riferimento globale, il campo degli spostamenti assunto è
invariante.
Infine modi simili di deformazione devono corrispondere ad autovalori multipli.
181
ELEMENTI ISOPARAMETRICI
ESEMPIO 5.3
Calcolare gli autovalori della matrice di rigidezza dell'elemento asta a due
nodi. Assumendo EA § l = 1 la matrice di rigidezza è data da:
[k] =
1 –1
–1 1
5.175
e gli autovalori della matrice di rigidezza si ricavano da:
det
1 –1 – l 0
–1 1
0 l
5.176
= 0
Il polinomio caratteristico è:
l 2 – 2l = 0
5.177
e gli autovalori:
l1 = 0
l2 = 2
5.178
ed i corrispondenti autovettori normalizzati sono:
1
{ f } 1T = ------- { 1
2
1}
1
{ f } 2T = ------- { 1
2
– 1}
5.179
L'elemento asta possiede un grado di libertà di moto rigido, lo spostamento
assiale; di conseguenza c'è un autovalore nullo a cui corrisponde un autovettore che rappresenta appunto un moto rigido. Il secondo autovettore rappresenta invece il modo di deformazione proprio di questo elemento
(allungamento/accorciamento) (fig. 5.28).
Fig. 5.28 – Modi di deformazione dell’elemento asta a 2 nodi.
182
ELEMENTI ISOPARAMETRICI
ESEMPIO 5.4
Calcolare gli autovalori della matrice di rigidezza dell'elemento asta a tre
nodi. Assumendo EA § l = 1 , la matrice di rigidezza, è data da:
[k] =
1
Ú–1
2
Ê x – 1---ˆ
Ë
2¯
Ê x 2 – 1---ˆ – 2x Ê x – 1---ˆ
Ë
Ë
4¯
2¯
Ê x 2 – 1---ˆ
Ë
4¯
2
Ê x + 1---ˆ – 2x Ê x + 1---ˆ d x
Ë
¯
Ë
2
2¯
1
1
– 2x Ê x – ---ˆ – 2x Ê x + ---ˆ
Ë
Ë
2¯
2¯
5.180
4x 2
che integrata analiticamente o numericamente in modo completo con il
metodo di Gauss con due punti di campionamento (si ha infatti 2m – 1 = 2 ;
m = 3 § 2 Þ 2 ) da:
[k] =
7
--3
1
--3
8
– --3
1
--3
7
--3
8
– --3
8
– --3
8
– --3
16
-----3
5.181
mentre sottointegrata, con un punto di campionamento, da:
1 –1 0
[ k ] = –1 1 0
0 0 0
5.182
Nel caso di integrazione analitica o numerica completa il polinomio caratteristico della matrice di rigidezza è:
l ( l 2 – 10l + 16 ) = 0
5.183
gli autovalori risultano quindi:
l1 = 0
l2 = 2
5.184
l3 = 8
ed i corrispondenti autovettori normalizzati sono:
1
{ f } 1T = ------- { 1 1 1}
3
1
{ f } 2T = ------- { 1 – 1 0}
2
1
{ f } 3T = ------- { 1 1 – 2}
6
5.185
183
ELEMENTI ISOPARAMETRICI
Nel caso di sottointegrazione numerica il polinomio caratteristico della
matrice di rigidezza è:
l 3 – 2l 2 = 0
5.186
gli autovalori risultano quindi:
l1 = 0
l2 = 0
l3 = 2
5.187
i primi due autovettori, corrispondenti agli autovalori nulli, hanno la terza
componente, corrispondente al nodo centrale, indeterminata, mentre le due
componenti, corrispondenti ai nodi di estremità, sono uguali; il terzo autovettore normalizzato è:
1
{ f } 3T = ------- { 1 – 1 0}
2
5.188
Si nota che nel caso di integrazione esatta si ha un solo autovettore nullo,
corrispondente all'effettivo grado di libertà di moto rigido dell'elemento.
Nel caso della sottointegrazione, avendo valutato in modo errato l'energia di
deformazione dell'elemento, si calcola invece un secondo autovettore nullo,
corrispondente ad un moto rigido spurio.
La figura 5.29 illustra i modi di deformazione relativi all'elemento asta a tre
nodi, in cui la matrice di rigidezza è integrata esattamente:
Fig. 5.29 – Modi di deformazione dell’elemento asta a 3 nodi.
La figura 5.30 illustra i modi di deformazione relativi ad un elemento trave
inflessa, calcolato con integrazione esatta; si notano i due modi relativi ai gradi di
libertà di moto rigido ed i due modi di deformazione, uno di flessione ed uno di
taglio.
La figura 5.31 illustra i modi di deformazione relativi ad un elemento isoparametrico a quattro nodi in condizioni di stato di deformazione piano, in cui la
matrice di rigidezza è stata calcolata con integrazione esatta utilizzando 2 x 2
punti campione; si notano i tre modi relativi ai gradi di libertà di moto rigido ed
i cinque modi di deformazione, due di flessione, uno di taglio, uno di estensione/compressione ed uno di estensione uniforme.
184
ELEMENTI ISOPARAMETRICI
Fig. 5.30 – Modi di deformazione dell’elemento trave inflessa.
Fig. 5.31 – Modi di deformazione dell’elemento a 4 nodi (integrazione esatta).
La figura 5.32 illustra invece i modi di deformazione relativi ad un elemento isoparametrico a quattro nodi in condizioni di stato di deformazione piano, in cui
la matrice di rigidezza è stata calcolata con sottointegrazione utilizzando 1 punto
campione; si notano i tre modi relativi ai gradi di libertà di moto rigido, due
modi spuri e tre modi di deformazione, uno di taglio, uno di estensione/compressione ed uno di estensione uniforme.
185
ELEMENTI ISOPARAMETRICI
Fig. 5.32 – Modi di deformazione dell’elemento a 4 nodi (sottointegrazione).
La figura 5.33 illustra i modi di deformazione relativi ad un elemento isoparametrico a quattro nodi in condizioni di stato di deformazione piano, in cui la
matrice di rigidezza è stata calcolata con integrazione selettiva utilizzando 2
punti campione per le componenti normali della deformazione ed 1 punto campione per la componente di taglio.
Fig. 5.33 – Modi di deformazione dell’elemento a 4 nodi (integrazione selettiva).
186
ELEMENTI ISOPARAMETRICI
5.5.5 Patch test
Il patch test, ideato da Irons, è un metodo per verificare se nell'elemento, o
meglio in un insieme di elementi di forma qualsiasi e comunque assemblati,
sono rispettate le condizioni di moto rigido e di stato di deformazione costante.
Questa verifica consiste nel assemblare diversi elementi, di forma qualsiasi, in
modo tale che almeno un nodo sia compreso all'interno della struttura così formata come ad esempio nella figura 5.34.
Se un elemento supera il patch test è garantita la convergenza al risultato esatto
all'aumentare della suddivisione in elementi.
Nella regione così definita si assuma un campo di spostamenti tale da produrre
uno stato di deformazione costante.
Il problema così definito viene ora risolto numericamente applicando su tutti i
nodi esterni spostamenti congruenti con il campo imposto. Se il risultato del calcolo fornisce, nella regione definita, uno stato di tensione costante ed uguale a
quello teorico allora l'elemento in questione ha superato il patch test.
Fig. 5.34 – Patch test.
È da notare che il superare o meno il patch test può dipendere fortemente dalla
geometria dell'elemento utilizzata; ad esempio l'elemento a quattro nodi con
modi incompatibili supera il patch test se ha forma di un parallelogramma, mentre non lo supera con forme diverse. Di conseguenza un solo patch test non è
sufficiente per garantire la convergenza dell'analisi.
187
ELEMENTI ISOPARAMETRICI
ESEMPIO 5.5
Nel caso della figura 5.34 sia definito il seguente campo di spostamenti:
1
y
u = ------ Ê x + ---ˆ
10 Ë
2¯
5.189
1 x
v = ------ Ê --- + yˆ
10 Ë 2 ¯
e le coordinate dei nodi siano:
NODO
1
2
3
4
5
6
7
8
X
0
40
40
0
10
30
25
12
Y
0
0
40
40
10
15
30
28
siano inoltre E = 100.000 , n = 0.3 e h = 1 .
Gli spostamenti teorici dei nodi interni sono:
NODO
5
6
7
8
U
1.50
3.75
4.00
2.60
V
1.50
3.00
4.25
3.40
e le deformazioni teoriche:
e xx = e yy = g xy = 0.1
5.190
Le condizioni al contorno per la soluzione numerica sono ricavate dalla legge
degli spostamenti definita ed applicate ai nodi 1, 2, 3, 4:
188
NODO
1
2
3
4
U
0.00
4.00
6.00
2.00
V
0.00
2.00
6.00
4.00
6. SOLIDI ASSIALSIMMETRICI
6.1 INTRODUZIONE
Un solido di rivoluzione è un corpo tridimensionale generato per rotazione di
una sezione piana attorno ad un asse (fig. 6.1). Un solido di rivoluzione è assialsimmetrico se le proprietà geometriche e del materiale sono indipendenti dalla
anomalia J . Se il carico è anch’esso assialsimmetrico il problema si riduce al
caso bidimensionale; ogni punto della sezione è caratterizzato dagli spostamenti
u (radiale) e w (assiale); l’analisi è essenzialmente coincidente con quella del problema piano con una variazione, dovuta alla presenza delle componenti circonferenziali della deformazione e della tensione.
Fig. 6.1 – Solido assialsimmetrico.
189
SOLIDI ASSIALSIMMETRICI
Nel caso di carico non assialsimmetrico, il problema può ancora essere ricondotto ad un caso bidimensionale, scomponendo il carico in serie di Fourier. Data
l’ortogonalità della serie di Fourier il problema è risolto ricorrendo al principio
di sovrapposizione degli effetti, sommando le n soluzioni ottenute applicando le
componenti simmetriche e antisimmetriche del carico, corrispondenti agli n termini dell’espressione in serie di Fourier del carico.
Un tipico elemento finito che descrive un solido di rivoluzione è caratterizzato
da circonferenze nodali, piuttosto che da punti nodali.
6.2 CARICO ASSIALSIMMETRICO
Fig. 6.2 – Elemento triangolare assialsimmetrico.
La figura 6.2 illustra un elemento triangolare assialsimmetrico. Nel caso di materiale e carico assialsimmetrico le uniche componenti del vettore spostamento { u }
in un punto generico P sono quella radiale, u e quella assiale, w, essendo la componente circonferenziale, v, nulla per simmetria:
{ u }T = { u w }
6.1
Il vettore spostamento {u } è espresso in funzione degli spostamenti nodali come:
n1 n2 n3 0 0 0
Ïu¸
{s}
Ì ý =
0 0 0 n1 n2 n3
Óvþ
190
6.2
SOLIDI ASSIALSIMMETRICI
dove le funzioni di forma dell’elemento triangolare, funzione delle coordinate
nodali dell’elemento sono date dalla 4.111:
1
n 1 = ------ [ ( x 2 y 3 – x 3 y 2 ) + ( y 2 – y 3 )x + ( x 3 – x 2 )y ]
2A
1
n 2 = ------ [ ( x 3 y 1 – x 1 y 3 ) + ( y 3 – y 1 )x + ( x 1 – x 3 )y ]
2A
1
n 3 = ------ [ ( x 1 y 2 – x 2 y 1 ) + ( y 1 – y 2 )x + ( x 2 – x 1 )y ]
2A
6.3
Le deformazioni dell'elemento sono definite dalla 3.94:
{ e } T = { e rr
e zz
e JJ
g rz }
6.4
con:
e rr =
e zz =
u
2p ( r + u ) – 2pr
e JJ = -------------------------------------- = -r
2pr
∂u
∂r
∂w
∂z
g rz
∂u ∂w
=
+
∂z ∂r
6.5
e in forma matriciale:
∂
∂r
{e} =
0
1
--r
∂
∂z
0
∂
----∂z { u } = [ ∂ ] [ n ] { s }
6.6
0
∂
----∂r
e derivando:
{e} = [b]{s}
6.7
Il vettore delle tensioni è:
{ s } T = { s rr
s zz
s JJ
t rz }
6.8
e la relazione tensioni deformazioni:
191
SOLIDI ASSIALSIMMETRICI
1 –n
Ï s rr ¸
Ô
Ô
n
Ô s zz Ô
E
- n
Ì
ý = (------------------------------------1 + n ) ( 1 – 2n )
Ô s JJ Ô
Ô
Ô
0
Ó t rz þ
n
n
1–n
n
n
1–n
0
0
0
0
0
1 – 2n
--------------2
Ï e rr ¸
Ô
Ô
Ô e zz Ô
Ì
ý
Ô e JJ Ô
Ô
Ô
Ó g rz þ
6.9
La matrice di rigidezza è data da:
[k] =
ÚV [ b ]T [ E ] [ b ] dV
6.10
l'elemento di volume dV è espresso da: dV = r dr dJ dz = r dr dA. Si ha:
[k] =
2p
ÚA Ú0
[ b ] T [ E ] [ b ]r dJ dA = 2p
ÚA [ b ]T [ E ] [ b ]r dA
6.11
dove A è l’area della sezione retta dell’elemento assialsimmetrico. L’integrale è di
solito valutato numericamente in quanto i termini di [b] sono funzioni logaritmiche.
Nel caso di elementi isoparametrici la matrice di rigidezza si ricava in modo analogo a quanto visto in 5.2. Ad esempio per elementi quadrangolari assialsimmetrici (fig. 6.3).
Fig. 6.3 – Elemento assialsimmetrico a 4 nodi.
gli spostamenti u e w e le funzioni di forma ni sono date dalle eq. 5.38 e 5.22:
192
SOLIDI ASSIALSIMMETRICI
4
 ni (x ,h )ui
u =
i=1
4
6.12
 ni (x ,h )wi
w =
i=1
1
n 1 = --- ( 1 – x ) ( 1 – h )
4
1
n 3 = --- ( 1 + x ) ( 1 + h )
4
1
n 2 = --- ( 1 + x ) ( 1 – h )
4
1
n 4 = --- ( 1 – x ) ( 1 + h )
4
6.13
Le deformazioni 6.5, sostituendo agli spostamenti u e v le espressioni delle funzioni di forma, sono espresse come:
n
e rr =
Â
∂n i
u
∂r i
Â
∂n i
w
∂z i
i=1
n
e zz =
i=1
n
e JJ =
Â
ni
---- u i
r
Â
Ê ∂n i u + ∂n i w ˆ
Ë ∂ z i ∂ r i¯
i=1
n
g rz =
i=1
6.14
e in forma matriciale, esplicitando il solo termine i-esimo:
{e} =
º
º
∂n i
∂r
º
º
0
º
º
ni
---r
º
º
∂n i
∂z
0 ºº Ï
Ô
Ô
∂n i
Ô
ºº Ô
∂z
Ì
Ô
0 ºº Ô
Ô
Ô
∂n i
Ó
ºº
∂r
¸
Ô
Ô
Ô
Ô
ý
wi Ô
Ô
º Ô
º Ôþ
º
º
ui
6.15
e la matrice di rigidezza [k]:
[k] =
1
1
p
1
1
Ú–1 Ú–1 Ú–p [ b ]T [ E ] [ b ]r dJ det [ J ] dx dh
= 2p
6.16
Ú–1 Ú–1 [ b ]T [ E ] [ b ]r det [ J ] d x d h
193
SOLIDI ASSIALSIMMETRICI
Adottando lo schema di integrazione numerica di Gauss:
m
[ k ] = 2p Â
m
 wi wj [ b ]ij [ E ] [ b ]ij rij det [ J ]ij
T
6.17
i = 1j = 1
dove m è l'ordine di integrazione, wi, wj sono i pesi di integrazione, [b]ij, rij e
det [ J ]ij sono calcolati nei punti di integrazione xi , hj. Alcuni termini sotto integrale sono funzione diretta di 1/ r. Nel caso di integrazione numerica con il
metodo di Gauss questi termini rimangono finiti poiché i punti di integrazione
non comprendono il caso r = 0. Tuttavia si può avere una forma indeterminata
nel calcolo delle deformazioni circonferenziali e JJ = u § r e delle tensioni per
punti giacenti sull’asse z dove r = 0. Un modo di ovviare a questo inconveniente
è quello di estrapolare i valori calcolati nei punti di integrazione di Gauss.
Un’altra opzione, che deriva dal fatto che per r = 0 e rr = e JJ , è quella di
sostituire, per il calcolo dei valori di tensione per punti giacenti sull’asse z, la riga
corrispondente a e JJ con quella corrispondente a err nella matrice [b].
6.3 CARICO NON ASSIALSIMMETRICO
In presenza di solidi assialsimmetrici soggetti a carichi non assialsimmetrici è
ancora possibile risolvere il problema in modo quasi-bidimensionale scomponendo i carichi e gli spostamenti nodali in serie di Fourier lungo la circonferenza.
Nel caso di carico non assialsimmetrico la componente tangenziale dello spostamento, v non può più essere considerata nulla, per cui:
{ u }T = { u w v }
6.18
Le componenti radiale circonferenziale e assiale del carico, F, sono:
{ f } T = { fr fz fc }
6.19
Queste possono essere espresse in serie di Fourier come:
m
fr =
s
 ( fz
s
 ( fc
s
h=0
m
fz =
h=0
m
fc =
h=0
194
m
 ( fr
)h cos ( hJ ) +
 ( fr
h=0
m
)h cos ( hJ ) +
 ( fz
h=0
m
)h sin ( hJ ) +
A
 ( fc
h=0
) h sin ( hJ )
A
A
)h sin ( hJ )
)h cos ( hJ )
6.20
SOLIDI ASSIALSIMMETRICI
Fig. 6.4 – Componenti simmetriche e antisimmetriche.
dove ciascun termine della serie è anche chiamata armonica ed h è l’ordine di ciascuna armonica. L’esponente s ed A denotano le componenti simmetriche ed
antisimmetriche del carico (fig. 6.4), dove simmetria e antisimmetria sono riferite rispetto ad un piano di riferimento J = 0 . La sommatoria è estesa agli m
termini necessari per descrivere correttamente la distribuzione di carico. In
modo analogo gli spostamenti sono espressi in serie di Fourier:
m
u ( r, z, J ) =
Â
h=0
m
w ( r, z, J ) =
v ( r, z, J ) =
m
u hs cos ( hJ )
+
 u hA sin ( hJ )
h=0
m
 whs cos ( hJ ) +  whA sin ( hJ )
h=0
m
h=0
m
h=0
h=0
6.21
 vhs sin ( hJ ) +  v hA cos ( hJ )
195
SOLIDI ASSIALSIMMETRICI
dove u hs , u hA, w hs , w hA, v hs e v hA sono i coefficienti di Fourier per l’armonica di
ordine h.
In forma matriciale la 6.21 è:
m
{u} =
Â
m
[ q hs ] { u hs }
+
h=0
 [ q hA ] { u hA}
6.22
s
Ï { u h } ¸ˆ
Ì A ý˜
Ó { u h } þ¯
6.23
h=0
ovvero:
m
Ê
{ u } = Â Á [ q hs ]
h = 0Ë
[ q hA ]
con:
s
s s
{ u hs } T = { u h w h v h }
{ u hA} T = { u hA w hA v hA}
6.24
e:
[ q hs ] =
[ q hA ] =
cos ( hJ )
0
0
0
cos ( hJ )
0
0
0
sin ( hJ )
6.25
sin ( hJ )
0
0
0
sin ( hJ )
0
0
0
cos ( hJ )
Per l’armonica h-esima gli spostamenti di un punto dell’elemento sono espressi
in funzione degli spostamenti nodali attraverso le funzioni di forma [n]:
{ u hs } = [ n ] { s hs }
{ u hA } = [ n ] {s hA}
6.26
e gli spostamenti complessivi:
m
{u} =
Â
h=0
m
[ q hs ] [ n ] { s hs }
+
 [ q A ] [ n ] {s hA }
6.27
h=0
ovvero:
m
Ê
{ u } = Â Á [ q hs ] [ n ] [ q hA ] [ n ]
h = 0Ë
196
Ï {s h } ¸ˆ
Ì A ý˜
Ó{s h }þ¯
s
6.28
SOLIDI ASSIALSIMMETRICI
Le deformazioni dell'elemento sono definite da:
e rr =
∂u
∂r
u
e JJ = -r
1 ∂u
∂v
g rJ = --- Ê
– vˆ +
r Ë ∂J ¯ ∂r
g zJ =
∂v 1 ∂w
+ --∂z r ∂J
e zz =
∂w
∂z
g rz =
∂u ∂w
+
∂z ∂r
6.29
e in forma matriciale:
∂
----∂r
1
--r
{e} =
0
0
0
1∂
--- ------r ∂J
∂
----∂z
0
1∂
--- ------r ∂J
0
0
1∂
--- ------r ∂J
∂
----∂r
0
∂
----∂z
{u} = [∂]{u}
∂ 1
----- – --∂r r
∂
----∂z
6.30
0
con:
{ e } T = { e rr
e JJ
e zz
g rJ
g zJ
g rz }
6.31
 [ ∂ ] [ qhA ] [ n ] {s hA }
6.32
e, sostituendo la 6.27:
m
{e} =
Â
m
[ ∂ ] [ q hs
] [ n ] {s hs }
h=0
+
h=0
ovvero:
m
{e} =
Â
m
[ b hs
] { s hs }
h=0
+
 [ bhA ] { s hA }
h=0
s
Ê
Ï { s h } ¸ˆ
s
A
{ e } = Â Á [ b h ] [ b h ] Ì A ý˜ =
Ó{ s h }þ¯
h = 0Ë
m
6.33
m
 [ b h ]{ s h }
6.34
h=0
Il vettore delle tensioni è:
197
SOLIDI ASSIALSIMMETRICI
{ s } T = { s rr
s JJ
s zz
t rJ
t zJ
t rz }
6.35
e la relazione tensioni deformazioni:
Ï s rr ¸
Ô
Ô
Ô s JJ Ô
Ô
Ô
Ô s zz Ô
E
Ì
ý = ------------------------------------( 1 + n ) ( 1 – 2n )
Ô t rJ Ô
Ô
Ô
Ô t zJ Ô
Ô
Ô
Ó t rz þ
1 –n
n
n
n
n
1–n
n
n
1–n
0
0
0
1 – 2n
--------------2
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Ï e rr ¸
Ô
Ô
Ô e JJ Ô
Ô
Ô
e
0
0 ÔÌ zz Ôý
Ô e rJ Ô
1 – 2n
--------------- 0 ÔÔ e ÔÔ
zJ
2
Ô
Ô
1 – 2n g
0 --------------- Ó rz þ
2
6.36
e, sostituendo per {e} la 6.33:
m
m
Ê
ˆ
{ s } = [ E ] Á Â [ b hs ] { s hs } + Â [ b hA ] { s hA } ˜
Ëh=0
¯
h=0
6.37
ovvero:
m
s
Ê
Ï { s n } ¸ˆ
s
A
{ s } = [ E ] Â Á [ b h ] [ b h ] Ì A ý˜ = [ E ] Â [ b h ]{ s h }
Ó{ s n }þ¯
h = 0Ë
h=0
m
6.38
La matrice di rigidezza è data da:
m
[k] =
m
  [ khl ]
6.39
h = 0l = 0
avendo indicato con [ k hl ] il contributo alla rigidezza delle armoniche di
ordine h e l:
[ k hl ] =
[ k hl ] =
ÚV [ bh ]T [ E ] [ bl ] d V
ÚV
[ b hs ]T
[ b hA ]T
[ E ] [ b ls ] [ b lA ] d V
l'elemento di volume dV è espresso da: dV = r dr dJ dz = r dr dA. Si ha:
198
6.40
6.41
SOLIDI ASSIALSIMMETRICI
[ k hl ] =
2p
ÚA Ú0
[ b hs ] T
[ b hA ] T
s
A
[ E ] [ bl ] [ bl ] r d J d A
6.42
dove A è l’area della sezione retta dell’elemento assialsimmetrico. La 6.42 si può
elaborare come:
[ k hl ] =
2p
ÚA Ú0
[ b hs ]T [ E ] [ b ls ] [ b hs ] T [ E ] [ b lA ]
[ b hA ]T [ E ] [ b ls ]
[ b hA ] T [ E ]
[ b lA ]
rdJdA
6.43
rdJdA
6.44
La matrice [k] è quindi:
[k] =
2p m
m
[ b hs ]T [ E ] [ b ls ] [ b hs ] T [ E ] [ b lA]
ÚA Ú0 hÂ= 0 lÂ= 0 [ bhA ]T [ E ] [ bls ]
[ b hA ] T [ E ] [ b lA]
Sono quindi da calcolare i seguenti integrali:
2p m
m
s
ÚA Ú0 hÂ= 0 lÂ= 0 [ bh ]T [ E ] [ bls ] r d J d A
2p m
m
2p m
m
2p m
m
6.45
A
A
ÚA Ú0 hÂ= 0 lÂ= 0 [ b h ]T [ E ] [ b l ] r d J d A
A
ÚA Ú0 hÂ= 0 lÂ= 0 [ bhs ]T [ E ] [ bl ] r d J d A
6.46
s
A
ÚA Ú0 hÂ= 0 lÂ= 0 [ b h ]T [ E ] [ bl ] r d J d A
Tutti i termini delle matrici in 6.46 contengono i seguenti integrali, il cui risultato è sempre nullo quando integrato tra 0 e 2p:
2p
Ú0
sin ( hJ ) cos ( hJ ) dJ = 0
h, l = 0, 1, 2, º, m
6.47
I termini simmetrici e antisimmetrici dell’espressione in serie sono quindi disaccoppiati, come d’altronde era chiaro visto che le serie di Fourier sono ortogonali.
In modo analogo si può constatare che gli integrali in 6.45 sono riconducibili a
due tipi:
199
SOLIDI ASSIALSIMMETRICI
2p
Ú0
2p
Ú0
Ï0
Ô
cos ( hJ ) cos ( hJ ) dJ = Ì p
Ô
Ó 2p
hπl
l = nπ0
l = n = 0
6.48
Ï0
Ô
sin ( hJ ) sin ( ( hJ ) dJ ) = Ì p
Ô
Ó 2p
hπl
l = nπ0
l = n = 0
La matrice [k] sarà quindi composta solo da sottomatrici sulla diagonale principale. Le singole armoniche possono quindi essere analizzate separatamente ed il
risultato finale può essere ottenuto sovrapponendo i risultati relativi alle m armoniche, come indicato in 6.27, 6.33, 6.37.
[k] =
ss ]
[ k 00
0
0
0
0
0
0
0
0
0
0
0
0
ss ]
[ k 11
0
0
0
0
0
0
0
0
0
0
0
0
0
º
0
0
0
0
0
0
0
0
0
0
0 [ k iiss ] 0
0
0
0
0
0
0
0
0
0
0
0
º
0
0
0
0
0
0
0
0
0
0
0
ss ]
0 [ k mm
0
0
0
0
0
0
0
0
0
0
0
0
AA ]
[ k 00
0
0
0
0
0
0
0
0
0
0
0
0
AA ] 0
[ k 11
0
0
0
0
0
0
0
0
0
0
º
0
0
0
0
0
º
0
0
AA ]
[ k mm
0
0
0
0
0
0
0
0
0
[ k iiAA ]
0
0
0
0
0
0
0
0
0
0
0
200
0
0
0
0
0
0
0
0
0
0
6.49
7. PIASTRE INFLESSE
7.1 RICHIAMI DELLA TEORIA DELLE PIASTRE INFLESSE
Questo capitolo si riferisce solamente alle piastre piane in campo lineare caricate
perpendicolarmente alla loro superficie media; se si suppone (fig. 7.1) che il
piano medio coincide con il piano di riferimento xy si avrà:
s xx = s yy = t xy = 0
per
z=0
7.1
Verranno inoltre considerate solo piastre sottili (una delle due dimensioni è trascurabile rispetto alle altre due) e quindi è possibile fare le seguenti ipotesi:
1. La componente di tensione szz normale al piano medio può essere trascurata rispetto alle altre componenti sxx, syy, txy (stato di tensione
piano)
2. Lo spostamento trasversale w non varia lungo lo spessore. (questa ipotesi è in contrasto con la prima)
3. Si ipotizza altresì, così come fatto per la teoria delle travi inflesse, che
sezioni piane perpendicolari al piano medio restino piane dopo deformazione. Nell’ipotesi di Kirchoff si ammette inoltre che il contributo
del taglio trasversale tyz , tzx sia nullo ed allora segmenti inizialmente
perpendicolari alla superficie media restano tali dopo la deformazione.
Nell’ipotesi di Mindlin si considera invece il contributo del taglio trasversale e quindi segmenti inizialmente perpendicolari alla superficie
media in generale non restano tali dopo deformazione
7.1.1 Sistemi di riferimento
L’elemento piastra è caratterizzato, nel suo piano medio, da uno spostamento
verticale w, e da due rotazioni ax , ay. È possibile definire spostamenti e rotazioni
secondo due sistemi di riferimento, il primo di tipo convenzionale ed il secondo
più utile per la caratterizzazione dell’elemento piastra (fig. 7.1a e fig. 7.1b). Si ha:
a x = – â y
a y = â x
7.2
201
PIASTRE INFLESSE
Fig. 7.1 – Sistemi di riferimento.
7.1.2 Tensioni e carichi per unità di lunghezza
È conveniente, così come per le travi, usare i carichi ed i momenti per unità di
lunghezza al posto delle tensioni, utilizzare cioè le caratteristiche delle sollecitazioni (in tal modo si effettua una preintegrazione lungo lo spessore della piastra).
Fig. 7.2 – Tensioni e carichi per unità di lunghezza.
Detto h lo spessore della piastra (fig. 7.2) si ha:
Mx =
Fx =
202
h§2
Ú–h § 2
h§2
sx z dz
Ú–h § 2 txz dz
My =
Fx =
h§2
Ú–h § 2
h§2
s y z dz M xy =
Ú–h § 2 tyz dz
h§2
Ú–h § 2 txy z dz
7.3
PIASTRE INFLESSE
Le tensioni espresse in funzione dei carichi per unità di lunghezza sono:
12M y
-z
s y = -----------h3
12M x
-z
s x = -----------h3
12M xy
-z
t xy = -------------h3
7.4
ed i valori massimi, in corrispondenza delle superficie di intradosso ed estradosso, valgono in valore assoluto:
6M x
s x = --------h2
6M y
s y = --------h2
6M xy
t xy = ----------h2
7.5
Normalmente gli sforzi di taglio txz e tyz sono trascurabili; essi raggiungono il
valor massimo per z = 0:
Fx
t xz = 1.5 ----h
Fy
t yz = 1.5 ----h
7.6
7.2 IPOTESI DI KIRCHOFF
7.2.1 Spostamenti e deformazioni
Nell’ipotesi di Kirchoff le sezioni inizialmente perpendicolari alla superficie
media restano tali dopo deformazione (fig. 7.3). Utilizzando il sistema di riferimento illustrato in figura 7.1b, si ha:
ax =
∂w
∂x
u = –z
ay =
∂w
= – za x
∂x
∂w
∂y
7.7
v = –z
∂w
= – za y
∂y
7.8
Fig. 7.3 – Spostamenti e rotazioni: ipotesi di Kirchoff.
203
PIASTRE INFLESSE
Le deformazioni sono espresse da:
e xx
2
∂a x
∂w
∂u
=
= –z
= –z 2 = –z kx
∂x
∂x
∂x
e yy
2
∂a y
∂w
∂v
=
= –z
= –z 2 = –z ky
∂y
∂y
∂y
g xy =
7.9
2
∂a ∂a
∂w
∂u ∂v
= – z Ê x + yˆ = – z 2
= – z k xy
+
Ë∂y
∂ xy
∂y ∂x
∂x ¯
avendo definito le curvature kx , ky , kxy come:
kx =
2
∂a x
∂w
=
∂x
∂x 2
ky =
2
∂a y
∂w
=
∂y
∂y 2
2
∂a ∂a
∂w
k xy = Ê ---------x + yˆ = 2
Ë ∂y ∂ x ¯
∂ xy
7.10
7.2.2 Relazioni tensioni-deformazioni
Per una piastra isotropa si ha, nell’ipotesi di Kirchoff:
Ï s xx ¸
Ô
Ô
E
Ì s yy ý = -------------21–n
Ô
Ô
t
Ó xy þ
1
n
n
1
0
0
0
0
1–n
-----------2
Ê Ï e xx ¸ Ï a x* T ¸ˆ
ÁÔ Ô Ô
Ô˜
Á Ì e yy ý – Ì a * T ý˜
Á Ô Ô Ô y Ô˜
Ë Ó g xy þ Ó 0 þ¯
7.11
in funzione delle curvature e dei carichi per unità di lunghezza si ha:
Ï e xx ¸
Ï c xx ¸
Ô Ô
Ô
Ô
Ì e yy ý = – z Ì c yy ý
Ô Ô
Ô
Ô
Ó g xy þ
Ó c xy þ
Ï s xx ¸
Ï Mx ¸
Ô
Ô
Ô
12z Ô
M
Ì s yy ý = ------Ì
ý
h 3 Ô xy Ô
Ô
Ô
Ó t xy þ
Ó M xy þ
e per il termine della temperatura:
2z
T = -----T 0
h
Fig. 7.4 – Gradiente di temperatura.
204
7.12
PIASTRE INFLESSE
Ï a x* T ¸
Ï a* ¸
Ï k xo ¸
2zT 0Ô x Ô
Ô
Ô
Ô
Ô
-Ì a * ý = – z Ì k yo ý
Ì a y* T ý = – ----------h Ô yÔ
Ô
Ô
Ô
Ô
Ó 0 þ
Ó 0þ
Ó 0 þ
2T 0
k xo = a x* --------h
2T 0
k yo = a y* --------h
7.13
Si ha quindi:
Ï Mx ¸
Ô
Ô
Eh 3
Ì M y ý = ------------------------12 ( 1 – n 2 )
Ô
Ô
Ó M xy þ
1
n
n
1
0
0
0
0
1–n
-----------2
Ê Ï k x ¸ Ï k xo ¸ˆ
ÁÔ
Ô˜
Ô Ô
Á Ì k y ý – Ì k yo ý˜ 7.14
ÁÔ
Ô˜
Ô Ô
Ë Ó k xy þ Ó 0 þ¯
7.3 IPOTESI DI MINDLIN
7.3.1 Spostamenti e deformazioni
Nell’ipotesi di Mindlin si ammette che la deformazione a taglio non sia più
trascurabile per cui una sezione piana perpendicolare alla superficie media
prima della deformazione rimane piana ma non più perpendicolare alla
superficie media (fig. 7.5). Utilizzando il sistema di riferimento illustrato in
figura 7.1b, si ha:
ax =
∂w
– g zx
∂x
ay =
∂w
– g zy
∂y
7.15
Fig. 7.5 – Spostamenti e rotazioni: ipotesi di Mindlin.
205
PIASTRE INFLESSE
u = – za x
v = – za y
7.16
Le deformazioni sono espresse da:
e xx =
e yy =
∂a ∂a
g xy = – z Ê x + yˆ = – zk xy
Ë∂y
∂x ¯
∂a x
∂u
= –z
= – zk x
∂x
∂x
∂a y
∂v
= –z
= – zk y
∂y
∂y
g yz =
∂w
– a y = – k yz
∂y
g zx =
∂w
– a x = – k zx
∂x
7.17
avendo definito le curvature kx , ky , kxy , kyz, kzx come:
kx =
k xy
∂a x
∂x
ky =
∂a ∂a
= Ê x + yˆ
Ë ∂y
∂x ¯
k yz =
∂a y
∂y
7.18
∂w
– ay
∂y
k zx =
∂w
– ax
∂x
7.3.2 Relazioni tensioni-deformazioni
Per una piastra isotropa si ha, nell’ipotesi di Mindlin:
Ï s xx ¸
Ô
Ô
E
Ì s yy ý = -------------21–n
Ô
Ô
t
Ó xy þ
1
n
n
1
0
0
0
0
1–n
-----------2
Ê Ï e xx ¸ Ï a x* T ¸ˆ
ÁÔ Ô Ô
Ô˜
Á Ì e yy ý – Ì a * T ý˜
Á Ô Ô Ô y Ô˜
Ë Ó g xy þ Ó 0 þ¯
7.19
in funzione delle curvature e dei carichi per unità di lunghezza si ha:
Ï e xx ¸
Ï kx ¸
Ô
Ô
Ô
Ô
Ì e yy ý = – z Ì k y ý
Ô
Ô
Ô
Ô
Ó g xy þ
Ó k xy þ
Ï s xx ¸
Ï Mx ¸
Ô
Ô
Ô
12z Ô
- Ì My ý
Ì s yy ý = ------3
h Ô
Ô
Ô
Ô
t
Ó xy þ
Ó M xy þ
Ï g yz ¸
Ï k yz ¸
Ì
ý = Ì
ý
Ó g zx þ
Ó k zx þ
Ï t yz ¸
1.5 Ï F y ¸
Ì
ý = -----h ÌÓ F ýþ
Ó t zx þ
x
Si ha quindi:
206
7.20
PIASTRE INFLESSE
1 n 0
0
0
Ï Mx ¸
n 1
0
0
0
Ô
Ô
1–n
Ô My Ô
0 0 -----------0
0
Ô
Ô
Eh 3
2
Ì M xy ý = ------------------------12 ( 1 – n 2 )
6(1 – n)
Ô
Ô
0
0 0
0 ------------------Ô Fy Ô
h2 c
Ô
Ô
Ó Fx þ
6(1 – n)
-------------------0 0
0
0
h2 c
Ê Ï k x ¸ Ï k ¸ˆ
Á Ô Ô Ô xoÔ˜
Á Ô k y Ô Ô k Ô˜
Á Ô Ô Ô yo Ô˜
Á Ì k xyý – Ì 0 ý˜
Á Ô Ô Ô Ô˜
Á Ô k yz Ô Ô 0 Ô˜
Á Ô Ô Ô Ô˜
Ë Ó k þ Ó 0 þ¯
7.21
zx
7.4 FORMULAZIONE DI RIGIDEZZA NELL’IPOTESI DI KIRCHOFF
Nell’ipotesi in cui si trascura il contributo del taglio alla deformazione (ipotesi di
Kirchoff ), appare chiaro che se viene assegnato il campo degli spostamenti w, è
necessario calcolarne le derivate seconde per ricavare le deformazioni, e che inoltre si dovrà formulare un elemento in classe C1, ovvero con continuità sulle derivate prime dello spostamento, dovendo assicurare la congruenza (anche lungo lo
spessore dell’elemento).
La cosa non è impossibile, ma diventa immediatamente piuttosto complicata.
Sia ammetta a tale scopo di avere un elemento a quattro nodi rappresentato in
figura 7.6 nello spazio naturale x, h.
Fig. 7.6 – Elemento piastra quadrangolare.
Dovendo interpolare in classe C1 una funzione f (x, h) si definirà innanzitutto
una funzione interpolatrice sul lato 1-2:
f 12 = G 1 ( x ) f 1 + G 2 ( x )
∂f 1
∂f 2
+ G3 ( x ) f2 + G4 ( x )
∂x
∂x
7.22
207
PIASTRE INFLESSE
essendo Gi ( x ) le funzioni di forma cubiche già viste per le travi inflesse, fi i valori
della funzione f nel nodo i ed ∂f i § ∂h il valore che la derivata di f secondo la
variabile h assume nel nodo i. Sul lato 1-2 occorre anche definire indipendentemente la derivata rispetto a h:
2
2
∂f 1
∂ f1
∂f 2
∂ f2
∂f 12
= G1 ( x )
+ G2 ( x )
+ G3 ( x )
+ G4 ( x )
∂h
∂ x ∂h
∂h
∂ x ∂h
∂h
7.23
Analogamente sul lato opposto, 4-3 si ha:
f 43 = G 1 ( x )f 4 + G 2 ( x )
∂f 4
∂f 3
+ G 3 ( x )f 3 + G 4 ( x )
∂x
∂x
2
7.24
2
∂f 4
∂f 3
∂ f4
∂ f3
∂f 43
= G1 ( x )
+ G2 ( x )
+ G3 ( x )
+ G4 ( x )
7.25
∂h
∂h
∂ x ∂h
∂ x ∂h
∂h
Interpolando secondo h:
f ( x, h ) = G 1 ( h ) f 12 + G 2 ( h )
∂f 12
∂f 43
+ G 3 ( h ) f 43 + G 4 ( h )
∂h
∂h
7.26
Sostituendo le 7.22, 7.23, 7.24, 7.25 nella 7.26 si ottiene la funzione di interpolazione bicubica nella quale appaiono come coefficienti i valori, calcolati nei nodi,
della funzione, di tutte le sue derivate prime e delle derivate seconde miste; in
totale 16 coefficienti.
Per il calcolo delle rotazioni si ha:
ax =
∂w ∂x ∂w ∂h
∂w
◊ + ------- ◊ -----=
∂ x ∂ x ∂h ∂x
∂x
∂w ∂x ∂w ∂h
∂w
=
◊ + ------- ◊ -----ay =
∂ x ∂ y ∂h ∂y
∂y
7.27
I valori delle derivate di x, h rispetto a x e y si possono ricavare numericamente
invertendo lo jacobiano, così come si è già visto nel caso degli elementi isoparametrici membranali.
Le deformazioni si ricavano dalla 7.10:
208
PIASTRE INFLESSE
2
2
2
2
2
2
2
2
2
2
2
2
cx =
∂ w ∂x ∂h ∂w ∂ h ∂ w Ê ∂hˆ 2
∂ w Ê ∂xˆ 2 ∂w ∂ x
------ ------ + ------2
+
+
+
∂x ∂x 2
∂ x∂h ∂x ∂x ∂h ∂ x 2 ∂ h 2 Ë ∂ x ¯
∂ x 2 Ë ∂ x¯
cy =
∂ w ∂x ∂h ∂w ∂ h ∂ w Ê ∂hˆ 2
∂ w Ê ∂xˆ 2 ∂w ∂ x
------ ------ + ------2
+
+
+
∂x ∂y 2
∂ x∂h ∂y ∂y ∂h ∂ y 2 ∂ h 2 Ë ∂ y ¯
∂ x 2 Ë ∂ y¯
c xy =
2
2
7.28
2
∂ w Ê ∂x ∂h ∂x ∂hˆ ∂w ∂ h
∂ w ∂x ∂x ∂w ∂ x
∂ w ∂h∂h
------ ------ +
-------------+
+2
+ ------+
∂ x ∂h Ë ∂ x ∂ y ∂ y ∂ x ¯ ∂h ∂ x∂y ∂ h 2 ∂x∂y
∂ x 2 ∂x ∂y ∂ x ∂ x∂y
Si intuisce, senza procedere oltre, che le cose si complicano assai rapidamente:
infatti questa non è una strada effettivamente seguita.
Continuità degli spostamenti implica congruenza di w e∂w § ∂n (derivata direzionale); questo rende difficile costruire funzioni di forma congruenti. Se si
usano come gradi di libertàw, ∂w § ∂x, ∂w § ∂y, vi sarà discontinuità di
∂ 2 w § ∂x∂y .
Si escogitano allora elementi con sottocampi aventi differenti funzioni di forma.
Una delle soluzioni possibili si ottiene rinunciando a lavorare in classe C 1; si possono definire funzioni interpolatrici indipendenti per lo spostamento e per le
rotazioni, lavorando così in classe C0. Occorre però che rotazioni e spostamento
siano effettivamente indipendenti; bisogna quindi abbandonare l’ipotesi di Kirchoff e questo può essere fatto tenendo conto dell’effetto del taglio sulla deformata della piastra.
Nei paragrafi seguenti si tratterà prima il problema, più semplice, relativo alla
trave per estendere successivamente il tutto al caso bidimensionale.
7.5 LA TRAVE DI TIMOSHENKO
Nella letteratura anglosassone si chiama trave di Timoshenko la trave inflessa
nella quale la rotazione della sezione è svincolata da quella della normale all’asse
neutro.
Si ha:
e xx =
g xy =
∂a z
∂u
= –y
∂x
∂x
∂v
– az
∂x
7.29
Lo spostamento v e la rotazione az sono interpolati indipendentemente così
l’elemento può rappresentare la deformazione dovuta al taglio. È tuttavia da
notare che questo elemento, se lineare, può simulare un momento costante ma
non uno variabile linearmente.
209
PIASTRE INFLESSE
Fig. 7.7 – Effetto del taglio sulla rotazione della sezione.
7.5.1 Elemento lineare
La figura 7.8 rappresenta l’elemento trave di Timoshenko lineare.
Fig. 7.8 – Elemento trave a 2 nodi.
Lo spostamento v e la rotazione az sono interpolati indipendentemente come:
Ïv ¸
v = n1 n2 Ì 1 ý
Ó v2 þ
Ïa ¸
a z = n 1 n 2 Ì z1 ý
Ó a z2 þ
con
210
7.30
PIASTRE INFLESSE
x
n1 = 1 – l
x
n2 = l
7.31
Definendo il vettore degli spostamenti come:
{ s } T = { v 1 a z1 v 2 a z2 }
7.32
le deformazioni exx e gxy sono espresse come:
e xx = – y
∂a z
∂n 2
∂n 1
= –y 0
{ s } = [ bF ] { s }
0
∂x
∂x
∂x
7.33
con
∂n 2
1
= --∂x
l
∂n 1
1
= – --l
∂x
7.34
La matrice di rigidezza dell’elemento può essere calcolata come somma di una
parte flessionale [kF ] ed una di taglio [kT ], essendo i due fenomeni disaccoppiati.
[ k ] = [ kF ] + [ kT ]
7.35
con:
[ kF ] =
[ kT ] =
ÚV [ bF ]T [ E ] [ bF ] dV
ÚV [ bT
]T[ E ][ b
7.36
T ]dV
Si ha:
[ kF ] = E Jz
Úl
0
0
0
0
Ê ∂n 1ˆ
Ë∂x ¯
0
0
0
0
∂n 1 ∂n 2
◊
∂x ∂x
0
2
0
0
∂n 1 ∂n 2
◊
∂ x ∂ x dx
0
Ê ∂n 1ˆ
Ë∂x ¯
7.37
2
211
PIASTRE INFLESSE
e:
Ê ∂n 1ˆ
Ë∂x ¯
GA
[ k T ] = -------c
Úl
2
– n1
Ê ∂n 1ˆ
Ë∂x ¯
∂n 1
–n1
∂x
∂n 1 ∂n 2
◊
∂x ∂x
–n2
∂n 1
∂x
∂n 1
∂x
–n1
∂n 1 ∂n 2
◊
∂x ∂x
2
∂n 2
∂x
∂n 1 ∂n 2
◊
∂x ∂x
–n1
∂n 2
∂x
Ê ∂n 2ˆ
Ë∂x ¯
–n2
–n2
∂n 1
∂x
n1 n2
2
–n2
∂n 2
∂x
∂n 2
∂x
Ê ∂n 2ˆ
Ë∂x ¯
dx
7.38
2
Integrando analiticamente (o numericamente con due punti di Gauss) si ha:
E Jz
[ kF ] = -------l
GA
[ k T ] = -------cl
0 0
0 1
0 0
0 –1
0
0
0
0
1
l
--2
l
--2
–1
l
--2
0
–1
0
1
7.39
–1
l2
---3
l
– --2
l
– --2
l2
---6
l
– --2
1
l
--2
l2
---6
l
– --2
7.40
l2
---3
Se si sottostima la componente di rigidezza dovuta al taglio, sottointegrando la
matrice di rigidezza [kT ] con un solo punto di Gauss si ottiene:
212
PIASTRE INFLESSE
1
l
--2
l
--GA 2
[ kT ] = -------cl
–1
l2
---4
l
– --2
l
– --2
l
--2
l2
---4
l
– --2
–1
1
l
--2
l2
---4
l
– --2
7.41
l2
---4
I vettori dei carichi nodali equivalenti ad un carico distribuito qv e quello dovuto
ad una differenza di temperatura 2DT fra intradosso ed estradosso sono rispettivamente:
Ïl
l ¸
{ f e } qT = q v Ì --- 0 --- 0 ý
2
2
Ó
þ
T
{ f e } DT
7.42
2E Jz a* DT
= --------------------------- { 0 1 0 – 1}
h
7.5.2 Elemento parabolico
La figura 7.9 rappresenta l’elemento trave di Timoshenko parabolico. Le leggi
dello spostamento v e della rotazione az sono in tal caso date da funzioni di
secondo grado, così come nel caso dell’elemento asta parabolico visto nel cap. 5,
eq. 5.65 e 5.66:
Ï v1 ¸
Ô Ô
v = n1 n2 n3 Ì v2 ý
Ô Ô
Ó v3 þ
7.43
con:
1
n 1 = – --- x ( 1 – x )
2
∂n 1
1
= – --- + x
2
∂x
1
n 2 = --- x ( 1 – x )
2
∂n 2
1
= --- + x
∂x
2
n3 = 1 – x 2
∂n 3
= – 2x
∂x
7.44
213
PIASTRE INFLESSE
Fig. 7.9 – Elemento trave a 3 nodi.
Supponendo il nodo centrale in mezzeria ( x3 = (x1 + x2)/2), lo Jacobiano ed il
suo determinante valgono, eq. 5.68:
l
[ J ] = --2
2
[ J ] –1 = --l
l
det [ J ] = --2
7.45
Ordinando il vettore degli spostamenti come:
{ s } T = { v 1 a z1 v 2 a z2 v 3 a z3 }
7.46
le deformazioni exx e gxy sono date da:
e xx = – y
∂n 1
∂n 2
∂n 3
∂a
0
0
= –y 0
{s}
∂x
∂x
∂x
∂x
∂n 1
∂n 2
∂n 3
2
= – y --- 0 ∂ x 0 ∂ x 0 ∂ x { s }
l
7.47
2
= – y --- 0 Ê – 1--- + xˆ 0 Ê 1--- + xˆ 0 ( – 2x ) { s }
l
Ë 2
¯ Ë2
¯
= [ bF ] { s }
g xy =
∂v
–a =
∂x
= 2--- Ê – 1--- + xˆ
¯
lË 2
= [ bT ] { s }
214
∂n 1
∂n 2
∂n 3
–n1
–n2
–n3 { s }
∂x
∂x
∂x
1
2 1
1
4
--- x ( 1 – x ) --- Ê --- + xˆ – --- x ( 1 + x ) – --- x – ( 1 – x 2 ) { s }
¯
2
l Ë2
2
l
7.48
PIASTRE INFLESSE
La matrice di rigidezza si calcola come somma dei due contributi flessionale e di
taglio così come per l’elemento lineare:
[ k ] = [ kF ] + [ kT ]
[ kF ] =
[ kT ] =
ÚV [ bF ]T [ E ] [ bF ] ( dA ) dx
7.49
ÚV [ bT ]T [ E ] [ bT ] ( dA ) dx
Integrando analiticamente i due termini si ottiene per la parte flessionale (stessi
termini dell’elemento asta parabolica 5.71):
0 0
7
0 --6
2E Jz 0 0
[ k F ] = -----------1
l
0 --6
0 0
4
0 – --3
0 0
1
0 --6
0 0
7
0 --6
0 0
4
0 – --3
0 0
4
0 – --3
0 0
4
0 – --3
0 0
8
0 --3
7.50
e per il contributo dovuto al taglio:
14
-------3l 2
1
--l
2
-------GAl 3 l 2
[ kT ] = --------2c
1
– ---3l
16
– -------3l 2
4
---3l
1
--l
4
-----15
1
---3l
1
– -----15
4
– ---3l
2
-----15
2
-------23l
1
---3l
14
-------3l 2
1
– --l
16
– -------23l
4
– ---3l
1
– ---3l
1
– -----15
1
– --l
4
-----15
4
---3l
2
-----15
16
– -------23l
4
– ---3l
16
– -------3l 2
4
---3l
32
-------23l
0
4
---3l
2
-----15
4
– ---3l
2
-----15
7.51
0
16
-----15
mentre una sottointegrazione (con due punti di Gauss) fornisce:
215
PIASTRE INFLESSE
1
2
14
16 4
1
-------2- --- -------2- – ---- – -------- ---3 l 3l 2 3 l
l
3l
3l
2
1
2
1
4
1
---- – --- – ---- ------3l
9
l
3l 9
9
2
1
14
4
1
16
-------- ---- -------- – --- – -------- – ---2
GAl 3l 2 3 l 3l 2
3
l
l
3l
[ kT ] = --------2c
2
4
2
1
1
1
-------– ---- – --- – --9
3l
9
3l
9
l
32
16 4
16
4
– -------2- – ---- – -------2- ---- -------2- 0
3
3 l 3l
l 3l
3l
2
8
2
4
4
----- – ---- -----0
9
9
3l
3l 9
7.52
I vettori dei carichi nodali equivalenti ad un carico uniformemente distribuito qv
e ad una distribuzione di temperatura 2 DT tra intradosso ed estradosso sono
rispettivamente:
qv l Ï 1 1 4 ¸
{ f e } qT = ------ Ì --- 0 --- 0 --- 0 ý
2Ó 3 3 3 þ
2E Jz a* DT
T = -------------------------- { 0 1 0 –1 0 0 }
{ f e } DT
h
ESEMPIO 7.1
Calcolare lo spostamento dell’estremo libero della mensola illustrata in figura
7.10, modellando la trave sia con un solo elemento lineare che con un solo
elemento parabolico.
Fig. 7.10 – Mensola.
216
7.53
PIASTRE INFLESSE
La soluzione esatta fornisce per la freccia v* e la rotazione a* dell’estremo
libero:
Pl
Pl 3
v* = – ------------ – c -------GA
3E Jz
Pl 2
a* = – -----------2E Jz
7.54
La formulazione di rigidezza nel caso dell’elemento lineare è:
l
l
1 --- – 1 --2
2
EJ
-------zl
0
0
0
0
0
1
0
–1
0
0
0
0
0
– 1 + GA
-------cl
0
1
Ï F v1 ¸
l l 2 l l 2 ÏV ¸
--- ----- – --- ----- Ô 1 Ô
Ô
Ô
2 3 2 6 ÔA Ô
Ô A z1 Ô
1
= Ì
ý
l
l Ì ý
Ô F v2 Ô
– 1 – --- 1 – --- Ô V 2 Ô
2
2 Ô Ô
Ô
Ô
Ó A2 þ
Ó A z2 þ
l l2 l l2
--- ----- – --- ----2 6 2 3
7.55
e, imponendo le condizioni al contorno V 1 = 0 ; A z1 = 0 ; F v2 = – P ; M z2 = 0
si ha:
E J 0 0 GA 1
-------z+ -------l 01
cl
l
– --2
l
– --2
l2
----3
Ï V 2¸
Ï –P ¸
Ì ý = Ì ý
Ó A 2þ
Ó0þ
7.56
Risolvendo il sistema si ricava:
4Pl c 3E Jz c + l 2 GA
V 2 = – ------------ ◊ -------------------------------------GA 12E Jz c + l 2 GA
A z2
3E Jz c + l 2 GA ˆ
2P c Ê
-˜
= ---------- Á 1 – 4 --------------------------------------GA Ë
12E Jz c + l 2 GA¯
7.57
2
Ponendo J z = r A , con r raggio giratore d’inerzia della sezione, si ha:
4Pl c 3E c ( r § l ) 2 + G
V 2 = – ------------- ----------------------------------------GA 12E c ( r § l ) 2 + G
3E c ( r § l ) 2 + G ˆ
2P c
A z2 = ---------- Ê 1 – 4 ----------------------------------------GA Ë
12E c ( r § l ) 2 + G¯
7.58
e nel caso di una trave snella r § l « 1 si ha:
4Pl c
V 2 = – -----------GA
6Pl c
A z2 = – -----------GA
7.59
Si nota che lo spostamento è dovuto esclusivamente al contributo del taglio,
mancando completamente quello flessionale; il modello della trave è più
rigido della realtà e si ha il fenomeno noto come locking.
217
PIASTRE INFLESSE
D’altronde in un caso del genere, essendo v e a z espressi da una legge lineare, ed essendo il taglio costante su tutta la lunghezza della trave, si ha che:
g xy =
∂v
– a z = cost
∂x
7.60
ma se lo spostamento v è lineare allora
∂v
= cost e allora anche a z = cost ;
∂x
ma a z = 0 all’estremo sinistro, quindi a z = 0 dovunque.
La componente dovuta al taglio prevale quindi su quella flessionale; se si sottostima il contributo alla rigidezza dovuto al taglio, utilizzando la matrice di
rigidezza [ k T ] sottointegrata, il sistema risolutivo 7.56 diventa:
E J 0 0 GA 1
--------z
+ -------cl
l 01
l
– --2
l
– --2
2
l
----4
Ï –P ¸
Ï V2 ¸
Ì ý = Ì ý
Ó0þ
Ó A2 þ
7.61
Si ricava:
Pl 3 Pl c
V 2 = – ------------ – ---------4E Jz GA
Pl 2
A z2 = – -----------2E Jz
7.62
risposta che approssima meglio il risultato esatto della trave incastrata, fornendo anzi per ciò che riguarda la rotazione il risultato esatto; con la sottointegrazione della parte di matrice di rigidezza dovuta al taglio (integrazione
selettiva) si evita, almeno in questo caso, il fenomeno del locking.
Utilizzando un solo elemento parabolico e imponendo le condizioni al contorno ( V 1 = 0 ; A z1 = 0 ; F v2 = – P ; M z2 = 0 ) la formulazione di rigidezza
con la matrice dovuta al taglio integrata esattamente da:
0 0
7
0 --2E Jz
6
-----------l
0 0
4
0 – --3
14
1
16
-------- – --- – -------0 0
l 3l 2
3l 2
4
4
4
1
0 – --- ---3 GAl – --l- ----15
3l
+ --------0 0
2c
32
16 4
– -------- ---- --------2
2 3l
8
3l
3l
-0
3
2
4
– ---- -----0
3l 15
4
– ---3l
2
-----15
0
16
-----15
Ï ¸
Ï V2 ¸
Ô –P Ô
Ô Ô
Ô0Ô
Ô A z2Ô
Ì ý = Ì ý 7.63
V
Ô0Ô
Ô 3Ô
Ô0Ô
Ô Ô
Ó þ
Ó A z3þ
Risolvendo il sistema si ricava:
2 2 2
2
4 2 2
Pl 240E Jz c + 84l GAE Jz c + l A G
V 2 = – ----------- --------------------------------------------------------------------------------------------4GA
E Jz ( 60E Jz c + l 2 GA )
A z2
218
Pl 2
= – -----------2E Jz
7.64
PIASTRE INFLESSE
2
Ponendo J z = r A dove r è il raggio giratore d’inerzia della sezione, si ha:
Pl 240E 2( r § l ) 4 c 2 + 84l 2GAE ( r § l ) 2 c + l 4 A 2 G 2
V 2 = – ----------- ----------------------------------------------------------------------------------------------------------------4GA
E ( r § l ) 2 ( 60E ( r § l ) 2 c + l 2 GA )
A z2
Pl 2
= – ----------2E Jz
7.65
Se si sottostima la componente di rigidezza dovuta al taglio il sistema risolutivo 7.2 diventa:
0 0
7
0 --2E Jz
6
-----------l
0 0
4
0 – --3
14
-------0 0
3l 2
4
1
0 – --3 GAl – --l+ --------0 0
2c
16
– -------2
8
3l
0 --3
4
– ---3l
1
16
4
– --- – -------- – ---l 3l 2 3l
2
4
2
-------9
3l
9
4
32
---- -------- 0
3l 3l 2
2
8
----0
9
9
Ï ¸
Ï V2 ¸
Ô –P Ô
Ô Ô
Ô0Ô
Ô A z2Ô
Ì ý = Ì ý 7.66
V
Ô0Ô
Ô 3Ô
Ô0Ô
Ô Ô
Ó þ
Ó A z3þ
Si ricava:
Pl 3 Pl c
V 2 = – ------------ – --------3E Jz GA
Pl 2
A z2 = – -----------2E Jz
7.67
risposta che coincide con il risultato esatto della trave incastrata.
Utilizzando i seguenti dati:
4
2
2
4
E = 210000 N/mm , n = 0.3 , G = 80770 N/mm , J = 8.3 ¥ 10 mm ,
P = 100 , l = 1000 mm , A = 100 mm2, e facendo variare il rapporto l § h
tra la lunghezza l e l’altezza h della trave tra 1 e 10, si ottengono (mantenendo costante il valore dell’area A ) i risultati illustrati nella figura 7.11 sia
per l’elemento lineare che per quello parabolico.
219
PIASTRE INFLESSE
Fig. 7.11 – Trave incastrata con carico concentrato.
ESEMPIO 7.2
Calcolare lo spostamento dell’estremo libero della mensola illustrata in figura
7.12, modellando la trave sia con un solo elemento lineare che con un solo
elemento parabolico.
Fig. 7.12 – Mensola.
La soluzione esatta fornisce per la freccia v* e la rotazione a* dell’estremo
libero:
Ml 2
v* = -----------2E Jz
220
Ml
a* = -------E Jz
7.68
PIASTRE INFLESSE
Utilizzando un solo elemento lineare integrato esattamente e risolvendo il
sistema si ricava:
6Ml 2 c
V 2 = --------------------------------------12E Jz c + GAl 2
12Ml c
A z2 = --------------------------------------12E Jz c + GAl 2
7.69
mentre un solo elemento lineare con integrazione selettiva fornisce, al pari di
un solo elemento parabolico integrato sia esattamente sia in modo selettivo:
Ml 2
V 2 = ------------2E Jz
Ml
A z2 = --------E Jz
7.70
Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nella
figura 7.13.
Fig. 7.13 – Trave incastrata con momento.
ESEMPIO 7.3
Calcolare lo spostamento dell’estremo libero della mensola illustrata in figura
7.14, modellando la trave sia con un solo elemento lineare che con un solo
elemento parabolico.
Fig. 7.14 – Mensola con carico distribuito.
221
PIASTRE INFLESSE
La soluzione esatta fornisce per la freccia v* e la rotazione a* dell’estremo
libero:
ql 2
ql 4
v* = ------------ + ----------- c
8E Jz 2GA
ql 3
a* = -----------6E Jz
7.71
Utilizzando un solo elemento lineare integrato esattamente e risolvendo il
sistema si ricava:
2ql 2 c 3E Jz c + GAl 2
V 2 = --------------- ◊ --------------------------------------2GA 12E Jz c + GAl
7.72
ql 3
q GAl 5
A z2 = ------------ – -------------------------------------------------------4E Jz 4E Jz ( 12E Jz c + GAl 2 )
mentre un solo elemento lineare con integrazione selettiva fornisce:
ql 2
ql 4
V 2 = ------------ + ----------- c
8E Jz 2GA
ql 3
A z2 = -----------4E Jz
7.73
Utilizzando un solo elemento parabolico integrato esattamente e risolvendo il
sistema si ricava:
360 ( E Jz c ) 2 + 96E Jz GAl 2 c + ( GAl 2 ) 2
V 2 = ql 2 ----------------------------------------------------------------------------------------------12E Jz GA ( 60E Jz c + GAl 2 )
ql 3
A z2 = -----------6E Jz
7.74
mentre un solo elemento parabolico con integrazione selettiva fornisce:
ql 2
ql 4
V 2 = ------------ + ----------- c
8E Jz 2GA
ql 3
A z 2 = -----------6E Jz
7.75
Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nella
figura 7.15.
222
PIASTRE INFLESSE
Fig. 7.15 – Trave incastrata con carico distribuito.
ESEMPIO 7.4
Si consideri ora il seguente esempio, illustrato in figura 7.16, e relativo al
problema iperstatico di una trave caricata in mezzeria e incastrata agli
estremi. Calcolare lo spostamento del punto centrale modellando la trave sia
con elementi lineari che con elementi parabolici.
Fig. 7.16 – Trave incastrata con carico in mezzeria.
223
PIASTRE INFLESSE
La soluzione esatta fornisce per la freccia v* e la rotazione a* dell’estremo
libero:
Pl 3
Pl
v* = – --------------- – ----------- c
24E Jz 2GA
a* = 0
7.76
Utilizzando due elementi lineari con integrazione esatta o selettiva e risolvendo il sistema si ricava:
Pl
V 2 = – ----------- c
2GA
A z2 = 0
7.77
mentre un solo elemento parabolico (a parità quindi di gradi di libertà
rispetto alla soluzione precedente) con integrazione esatta o selettiva fornisce:
3 Pl
V 2 = – ------ ◊ -------- c
16 GA
Az2 = 0
7.78
Si nota che manca totalmente la parte dovuta alla flessione pura. Nel caso di
una trave a sezione rettangolare di altezza h e larghezza b si hanno i seguenti
rapporti rispetto alla soluzione esatta, rispettivamente per gli elementi lineari
e parabolici:
G l 2
v*
------ = 1 + ------- Ê -- ˆ
V2
E c Ë h¯
7.79
8
G l 2
v*
------ = --- Ê 1 + ------c- Ê -- ˆ ˆ
E Ë h¯ ¯
3Ë
V2
7.80
e si ha quindi il risultato esatto solo nel caso trave a snellezza nulla e di due
elementi lineari.
Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nella
figura 7.17.
Fig. 7.17 – Trave incastrata con carico concentrato in mezzeria.
224
PIASTRE INFLESSE
La figura 7.18 illustra i risultati nel caso di 4 elementi TLE o TLS e 2 elementi
TPE o TPS.
Fig. 7.18 – Trave incastrata con carico concentrato in mezzeria.
7.6 MODELLI DISCRETI (LEGAME DISCRETO TRA FLESSIONE E TAGLIO)
Si è visto negli esempi precedenti che nel caso di travi sottili, dove l’effetto del
taglio può essere trascurato rispetto al fenomeno flessionale, l’elemento a tre
nodi, con ordine di interpolazione parabolico, fornisce risultati migliori dell’elemento lineare, a meno che non si adottino tecniche particolari nel caso di questi
ultimi (ad esempio l’integrazione selettiva).
Una spiegazione di ciò può essere ottenuta considerando l’energia potenziale elastica della trave nella sua formulazione di trave di Timoshenko; separando la
parte flessionale da quella a taglio si ha, dalle 7.35 e 7.36 o 7.49:
E Jz
P = -------2
2
Ê ∂a zˆ dx + GmA
-----------Ë
¯
2c
0 ∂x
Ú
l
2
Ê ∂v – a ˆ dx
z¯
Ë
0 ∂x
Ú
l
7.81
e dividendo per 1/2EJz:
225
PIASTRE INFLESSE
P =
GA
Ê ∂a zˆ dx + ----------Ë
¯
E Jz c
0 ∂x
Ú
l
2
2
Ê ∂v – a ˆ dx
z¯
Ë
0 ∂x
Ú
l
7.82
Il termine GA/EJz che moltiplica la quota parte di energia potenziale dovuta al
taglio assume valori molto grandi nel caso di una trave snella ( l § h Þ • ) può
essere interpretato come un fattore di penalty. È quindi importante che il
modello utilizzato (il tipo di elemento e le funzioni di forma scelte) sia in grado
di simulare piccole deformazioni a taglio per grandi valori di GA/EJz in modo
da annullare il contributo di energia elastica dovuto al taglio. Se invece, il
secondo integrale a secondo membro non è piccolo, allora nasce il problema del
locking, perché il termine GA/EJz funziona da moltiplicatore dell’energia di
taglio. Questo è quello che succede nel caso dell’elemento lineare a due nodi
integrato esattamente che quindi non può e non deve essere usato nel caso di
travi snelle. Queste conclusioni sono applicabili anche alle piastre inflesse.
Un modo di superare questo inconveniente è quello di sottointegrare la energia
dovuta al taglio, rendendo così singolare la matrice di rigidezza al taglio. Esiste
un ulteriore metodo, detto Kirchoff discreto, per superare questo problema. Questo metodo consiste nell’imporre l’uguaglianza della rotazione totale della
sezione e derivata dello spostamento in un numero discreto di punti.
Si consideri la trave illustrata in figura 7.19.
Fig. 7.19 – Elemento trave a 3 nodi.
Le funzioni di interpolazione siano le 7.44. Si possono eliminare le due variabili
v3 e a3 relative al nodo centrale e imporre quindi due condizioni di Kirchoff.
Scelti due punti x = ± q si impone in essi che:
Ï v1 ¸
2 ∂n 1 ∂n 2 ∂n 3 Ô Ô
∂v
= --- -------- ◊ -------- ◊ -------- Ì v 2 ý
l ∂x ∂x ∂x Ô Ô
∂x
Ó v3 þ
sia uguale a:
226
7.83
PIASTRE INFLESSE
Ï a1 ¸
Ô Ô
a = n1 n2 n3 Ì a2 ý
Ô Ô
Ó a3 þ
7.84
Scrivendo le due equazioni e separando le variabili v3 e a3, si ottiene:
∂n 1 l
∂n 3 l
∂n
--- n 1 2
– --- n 3
∂x 2 ∂x
2 Ï v3 ¸
∂x
Ì ý =
∂n 3 l
∂n 1 l
∂n 2
Ó a3 þ
– --- n 3
-n
1
2
∂x
∂x 2 ∂x
Ï v1 ¸
l
--- n 2 Ô Ô
2 Ô a1 Ô
Ì ý
Ô v2 Ô
l
--- n 2 Ô Ô
2 Ó a2 þ
7.85
ovvero:
Ï v1 ¸
Ô Ô
Ôa Ô
Ï v3 ¸
[a]Ì ý = [b]Ì 1 ý
Ô v2 Ô
Ó a3 þ
Ô Ô
Ó a2 þ
7.86
Se la matrice [a] è invertibile:
Ï v1 ¸
Ô Ô
Ô a1 Ô
Ï v3 ¸
Ì ý = [ a ] –1 [ b ] Ì ý
Ô v2 Ô
Ó a3 þ
Ô Ô
Ó a2 þ
7.87
e, mediante la 7.87 si potrà definire una matrice [t ] che permette di eseguire la
trasformazione:
{ v1 a1 v2 a2 v3 a3 } T = [ t ] { v1 a1 v2 a2 } T
7.88
ovvero:
{ s } = [ t ] { sr }
7.89
dove { s r } è il vettore delle variabili ridotto.
L’equazione dei lavori virtuali diventa allora:
{ ds r } [ t ]T [ k F ] [ t ] { s r } = { ds r } [ t ] T { f }
7.90
227
PIASTRE INFLESSE
Calcolando la 7.85 in x = ± q :
1
--- – q
2
1
--- + q
2
l
– 2q – --- ( 1 – q 2 ) Ï v ¸
2
3
Ì ý =
a
l
2q – --- ( 1 – q 2 ) Ó 3 þ
2
l
1
--- ( q 2 – q ) – q – --4
2
l 2
1
--- ( q + q ) q – --4
2
l
--- ( q 2 + q )
4
l
--- ( q 2 – q )
4
Ï v1 ¸
Ô Ô
Ô a1 Ô
Ì ý
Ô v2 Ô
Ô Ô
Ó a2 þ
7.91
che produce:
Ï v3 ¸
Ì ý =
Ó a3 þ
Ï v1 ¸
Ô Ô
Ô a1 Ô
Ì ý
1
1
q2
q2
- – ---------------------- Ô v 2 Ô
– --------------------- – ---------------------- -------------------l( 1 – q2 ) 2( 1 – q2 ) l( 1 – q2 ) 2( 1 – q2 ) Ô Ô
Ó a2 þ
1
--2
l
--8
1
--2
l
– --8
7.92
la matrice [t ] diventa:
[t] =
1
0
1
--2
0
1
l
--8
0
0
1
--2
0
0
l
– --8
1
1
q2
q2
- – ---------------------– --------------------- – ---------------------- -------------------2
2
2
l( 1 – q ) 2( 1 – q ) l( 1 – q ) 2( 1 – q2 )
0
0
1
0
0
0
0
1
7.93
e la matrice [kF ] è la 7.50. La matrice di rigidezza [k] coniugata al vettore ridotto
delle variabili nodali è:
[ k ] = [ t ] T [ kF ] [ t ]
7.94
Il coefficiente k11 della matrice [k ] è:
16EJ
k 11 = ----------------------------3
2l ( 1 – q 2 ) 2
7.95
mentre quello della matrice di rigidezza dell’elemento hermitiano (a due nodi):
12EJ
k 11 = ----------l3
7.96
Questi due coefficienti sono uguali se:
36 ( 1 – q 2 ) 2 = 16
che comporta:
228
7.97
PIASTRE INFLESSE
1
q = ± ------3
7.98
È un fatto notevole che questi siano proprio i valori delle ascisse dei punti di
integrazione di Gauss. Si può verificare che in tali ascisse anche tutti gli altri
coefficienti, e non solo l’elemento k11, coincidono con quelli della matrice hermitiana.
Se ne deduce che, accettando di trascurare gli effetti del taglio, i risultati migliori
si ottengono imponendo le condizioni di Kirchoff proprio nei due punti di integrazione di Gauss.
7.7 L’ELEMENTO PIASTRA DI MINDLIN
7.7.1 Matrice di rigidezza
Analogamente a quanto visto per la trave di Timoshenko, la teoria di Mindlin
considera la rotazione del piano medio separata dalla rotazione della sezione. Si
hanno le relazioni viste in 7.3.1:
e xx
∂a x
∂u
=
= –z
= – zk x
∂x
∂x
e yy
∂a y
∂v
=
= –z
= – zk y
∂y
∂y
∂a ∂a
g xy = – z Ê x + yˆ = – zk xy
Ë ∂y
∂x ¯
g yz =
∂w
– a y = – k yz
∂y
g zx =
∂w
– a x = – k zx
∂x
7.99
Lo spostamento w e le rotazioni ax e ay sono descritti indipendentemente come:
n
n
w =
Â
n i ( x, h )w i
i=1
ax =
Â
i=1
n
n i ( x, h )a x
i
ay =
 ni ( x, h )ay
i
7.100
i=1
ovvero:
{u} = [n]{s}
7.101
La matrice di rigidezza [k ] si può ottenere come somma della parte flessionale e
di quella dovuta al taglio, poiché i due fenomeni sono disaccoppiati:
[ k ] = [ kF ] + [ kT ] =
ÚV [ bF ]T [ EF ] [ bF ] dV + ÚV [ bT ] T [ ET ] [ bT ] dV
7.102
È conveniente riferirsi a momenti per unità di lunghezza e curvature, cioè alle
caratteristiche di sollecitazione (momenti flettenti e forze di taglio). Le deformazioni sono espresse da:
229
PIASTRE INFLESSE
∂a x
kx =
=
∂x
ky =
k xy
∂a y
=
∂y
n
Â
i=1
n
∂n i
a
∂ x xi
∂n i
 ∂ y ay
∂a x ∂a y
=
+
=
∂y
∂x
k yz =
k zx =
i
i=1
∂w
+ ay =
∂y
∂w
+ ax =
∂x
n
∂n i
∂n i
 ÊË ∂ y ax + ∂ x ay ˆ¯
i
7.103
i
i=1
n
∂n i
 ÊË ∂ y wi + ni ay ˆ¯
i
i=1
n
∂n i
 ÊË ∂ x wi + ni ax ˆ¯
i
i=1
e in forma matriciale:
Ï
Ô
Ô
Ô
Ì
Ô
Ô
Ô
Ó
kx ¸
Ô
ky Ô
Ô
k xy ý =
Ô
k yz Ô
Ô
k zx þ
∂n i
∂x
0
0
0
º
0
º
º
º
º
º
º
∂n i
∂y
º
º
∂n i
∂y
∂n i
∂x
º
º
∂n i
∂y
0
–ni
º
º
∂n i
∂x
–ni
0
º
º
º
º
0
º
º
º
Ï
Ô
Ô
Ô
Ô
Ô
Ì
Ô
Ô
Ô
Ô
Ô
Ó
¸
Ô
Ô
Ô
Ô
Ô
a xi ý
Ô
a yi Ô
Ô
ºÔ
º Ôþ
º
º
wi
7.104
Ï º ¸
Ô º Ô
º º [ k F, i ] º º Ô
[ bF ]
Ô
Ï {k} ¸
{s}
Ì { si } ý =
Ì
ý =
[ bT ]
º º [ k T, i ] º º Ô
Ô
Ó {g } þ
Ô º Ô
Ó º þ
230
7.105
PIASTRE INFLESSE
e le derivate
∂n i ∂n i
,
si calcolano come in 5.51.
∂x ∂y
La relazione carichi curvature è:
Eh 3
[ E F ] = -----------------------12 ( 1 – n 2 )
1
n
n
1
0
1 – n2
0 -------------2
0
0
7.106
Eh
[ E T ] = ------------------------ 1 0
2(1 + n) c 0 1
e le matrici di rigidezza sono normalmente calcolate mediante integrazione
numerica:
m
[ kF ] =
m
  wi wj [ bF ]Tij [ EF ] [ bF ]ij hij det [ J ]ij
i = 1j = 1
m m
[ kT ] =
7.107
  wi wj [ bT ]ijT [ ET ] [ bT ]ij hij det [ J ]ij
i = 1j = 1
7.7.2 Carico nodale equivalente a carico superficiale distribuito
Sia qv il carico per unità di superficie agente perpendicolarmente alla superficie
media dell’elemento piastra. Il vettore del carico nodale equivalente è:
{ fe }q =
ÚA [ n ]T { q } dA
7.108
con:
{ f e } qT = { F z M x M y }
i
i
i
{ q }T = { q
0
7.109
0}
e [n] è la matrice delle funzioni di forma dell’elemento 7.101. Il calcolo del
carico nodale equivalente si esegue numericamente come:
m
{ fe }q =
m
  [ nij ] { qij }wi wj
det [ J ] ij
7.110
i = 1j = 1
231
PIASTRE INFLESSE
ESEMPIO 7.5
È data una piastra circolare non forata, incastrata al bordo esterno, di raggio
esterno R = 50 mm e spessore costante h = 2mm . La piastra realizzata in
2
acciaio, con modulo elastico E = 200000 N/mm e coefficiente di Poisson
n = 0.3 , è soggetta ad un carico uniformemente distribuito pari a
2
p = 1 N/mm .
Calcolare la deformata ed i momenti flettenti radiale e circonferenziale. La
soluzione classica del problema fornisce i seguenti risultati:
4h 2
12 ( 1 – n 2 )p
2 – r 2 ) 2 + ---------------- ( R2 – r2 )
w = ---------------------------(
R
(1 – n)
64Eh 3
p
m r = ------ [ R 2 ( 1 + n ) – r 2 ( 3 + n ) ]
16
p
m c = ------ [ R 2 ( 1 + n ) – r 2 ( 1 + 3n ) ]
16
7.111
rispettivamente per la freccia, il momento radiale e quello circonferenziale.
Data la simmetria del problema è stato analizzato solo un quarto di piastra. È
stato utilizzato l’elemento piastra di Mindlin a 8 nodi e sono state esaminate
due suddivisioni, rispettivamente con 5 e 9 elementi e 72 e 120 gradi di
libertà (fig. 7.20)
Fig. 7.20 – Schematizzazione di piastra circolare: a. 5 elementi, b. 9 elementi.
La tabella 7.1 riassume i risultati ottenuti dal calcolo numerico confrontati
con quelli analitici.
Tab. 7.1 – Valori massimi della freccia e dei momenti flettenti
232
soluzione
w
mr ( mc )
ANALITICO
0.669
203.1
FEM
5 EL.
0.659
205.8
FEM
9 EL.
0.667
203.4
PIASTRE INFLESSE
La figura 7.21 illustra la deformata calcolata con le due suddivisioni, mentre
in figura 7.22 è rappresentato l’andamento dei momenti flettenti.
Fig. 7.21 – Deformata (esempio 7.1).
Fig. 7.22 – Momento radiale e circonferenziale (esempio 7.1).
Per la stessa piastra è stata infine esaminata l’influenza del rapporto h § R
sulla precisione del calcolo; la figura 7.23 riporta i risultati ottenuti per h § R
variabile da 0.01 a 0.4.
Fig. 7.23 – Andamento della freccia in funzione del rapporto h/R.
233