Elementi finiti - Parte III

progettodidattica in rete

prog

etto

dida

ttica

in re

tePolitecnico di Torino, maggio 2002

Dipartimento di Meccanica

Elementi finitiParte III

A. Gugliotta

otto editore

otto

IPERTESTI Questo file contiene elementi di ipertesto (quando si ritiene possano essere d'aiuto alla lettura) segnalati mediante l'uso di oggetti colorati: riquadro blu: figure, tabelle sottolineato arancio: rimando a note, esempi sottolineato verde: rimando a equazioni Una selezione sul riferimento (es. n. equazione o figura) visualizza l'informazione a cui si fa riferimento (es. equazione o figura stesse) Una selezione sull'informazione (es. equazione o figura) riporta alla pagina di lettura. FORMATO Il formato della pagina e' 17x24 cm, come indicato dai crocini di taglio sulla copertina. Per stampare il documento mantenendo le corrette dimensioni dell'area stampata bisogna selezionare l'opzione di stampa che permette di non adattare l'area di stampa alle dimensioni della pagina. La procedura dipende dal tipo di stampante utilizzata.

ELEMENTI FINITIParte III

A. GUGLIOTTA

POLITECNICO DI TORINO

WWW.POLITO.IT

INDICE – III

5. ELEMENTI ISOPARAMETRICI ........................................135

5.1 SISTEMA DI RIFERIMENTO LOCALE E FUNZIONI DI FORMA .... 135

Elementi monodimensionali ...........................................................136Elementi triangolari ........................................................................139Elementi quadrangolari ...................................................................142Elementi solidi tetraedri ..................................................................147Elementi solidi parallelepipedi ........................................................ 149

5.2 CALCOLO DELLA MATRICE DI RIGIDEZZA ............................... 152

5.3 CARICHI NODALI EQUIVALENTI ............................................. 160

Carichi di linea ...............................................................................160Carichi di volume ...........................................................................164Effetti termici.................................................................................. 164

5.4 CALCOLO DELLE TENSIONI .................................................... 165

5.5 PROBLEMI RELATIVI AGLI ELEMENTI ISOPARAMETRICI ........... 169

Scelta dell'ordine di integrazione .....................................................169Distorsione degli elementi............................................................... 171Integrazione selettiva e modi incompatibili .....................................174Autovalori della matrice di rigidezza ................................................180Patch 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 ......................................................................201Tensioni e carichi per unità di lunghezza .........................................202

7.2 IPOTESI DI KIRCHOFF ............................................................. 203

Spostamenti e deformazioni ............................................................203Relazioni tensioni-deformazioni ......................................................204

7.3 IPOTESI DI MINDLIN .............................................................. 205

Spostamenti e deformazioni ............................................................205Relazioni tensioni-deformazioni ...................................................... 206

7.4 FORMULAZIONE DI RIGIDEZZA NELL’IPOTESI DI KIRCHOFF .. 207

7.5 LA TRAVE DI TIMOSHENKO.................................................... 209

Elemento lineare .............................................................................210Elemento parabolico .......................................................................213

7.6 MODELLI DISCRETI (LEGAME DISCRETO TRA FLESSIONE E TAGLIO) ...................... 225

7.7 L’ELEMENTO PIASTRA DI MINDLIN ........................................ 229

Matrice di rigidezza......................................................................... 229Carico nodale equivalente a carico superficiale distribuito ............... 231

ii

135

5. ELEMENTI ISOPARAMETRICI

5.1 S

ISTEMA

DI

RIFERIMENTO

LOCALE

E

FUNZIONI

DI

FORMA

Per semplificare la caratterizzazione di elementi di forma geometrica complessa siintroduce 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 unelemento di forma qualsiasi nel sistema cartesiano (

x , y , z

) e un elemento diforma semplice nel sistema

naturale

, o

locale

, (

x

,

�

,

z

).

Un generico elemento nel sistema fisico (

x, y, z

) viene cioè pensato come la tra-sformazione di un elemento di forma semplice, e sempre lo stesso, nel sistemalocale (

x

,

�

,

z

); nei casi mono-, bi- e tri-dimensionale si hanno rispettivamentele seguenti relazioni:

5.1

dove

x

i

,

y

i

,

z

i

(

i

= 1,

n

) sono le coordinate cartesiane degli

n

nodi dell'elemento,

n

i

sono le funzioni di trasformazione o funzioni di forma definite nel sistema dicoordinate naturale (

x

,

�

,

z

) dell'elemento.

Le funzioni di forma

n

i

sono in numero pari ai nodi dell'elemento e variano tra 0e 1; esse assumono valore 1 in corrispondenza del nodo

i

di coordinate (

x

i

,

y

i

,

z

i

)e valore nullo in corrispondenza degli altri nodi dell'elemento.

x ni x( )xii 1=

n

Â=

x ni x h,( )xii 1=

n

Â=

y ni x h,( )yii 1=

n

Â=

x ni x h z,,( )xii 1=

n

Â=

y ni x h z,,( )yii 1=

n

Â=

z ni x h z,,( )zii 1=

n

Â=

ELEMENTI

ISOPARAMETRICI

5.1.1 Elementi monodimensionali

Si consideri un elemento a due nodi di lunghezza

l

(fig.

5.1

):

la posizione di un punto

P

sul segmento può essere individuata dalla sua coordi-nata fisica

x

P

oppure dalle misure dei segmenti

l

1

ed

l

2

o dalle loro misure adi-mensionali:

5.2L1

l1l---= L2

l2l---=

Fig. 5.1 – Elemento a 2 nodi: sistema di riferimento naturale.

essendo naturalmente:

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:

5.4

L1 L2+ 1=

x L1x1 L2x2+=

Fig. 5.2 – Elemento a 2 nodi: sistema di riferimento naturale.

136

Volendo esprimere la coordinata fisica x in funzione di una sola coordinata natu-rale x, variabile tra 0 e 1 (fig. 5.2) si effettua il cambio di variabile:

5.5L1 1 x–= L2 x=

ELEMENTI

ISOPARAMETRICI

e quindi:

5.6

È usuale riferire l'origine del sistema naturale nel baricentro geometrico dell'ele-mento (fig.

5.3

), con la coordinata naturale

x

variabile tra –1 e 1. Si avrà allora:

5.7

Fig. 5.3 –

Elemento a 2 nodi: sistema di riferimento naturale.

x 1 x–( )x1 xx2+=

L11 x–

2-----------= L2

1 x+2

------------=

5.8x1 x–

2-----------x1

1 x+2

------------x2+=

5.9

avendo definito le funzioni di forma n1 e n2:

5.10

La figura 5.4 mostra l'andamento delle funzioni di forma n1 e n2.

x n1x1 n2x2+=

n11 x–

2-----------= n2

1 x+2

------------=

Fig. 5.4 – Funzioni di forma per l'elemento monodimensionale a 2 nodi.

137

Per un elemento a tre nodi (fig. 5.5) la coordinata x è data da:

5.11x n1x1 n2x2 n3x3+ +=

ELEMENTI ISOPARAMETRICI

Fig. 5.5 – Elemento monodimensionale a 3 nodi.

dove:

5.12

e ponendo l'origine del sistema di riferimento naturale nel baricentro geometricodell'elemento si ha:

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.2riassumono le funzioni di forma e le loro derivate rispettivamente per gli ele-menti monodimensionali a 2 nodi e a 3 nodi.

n1 L1 2L1 1–( )= n2 L2 2L2 1–( )= n3 4L1L2=

n112---– x 1 x–( )= n2

12---x 1 x+( )= n3 1 x2–=

Tab. 5.1 – Funzioni di forma e loro derivate per l'elemento monodimensionale a 2 nodi

NODO 1 2

n 12--- 1 x–( ) 1

2--- 1 x+( )

∂n ∂x§ 12---– 1

2---

138


Tab. 5.2 – Funzioni di forma e loro derivate per l'elemento monodimensionale a 3 nodi

NODO 1 2 3

n 12---– x 1 x–( ) 1

2---x 1 x+( ) 1 x–( )2

∂n ∂x§ 12---– x+ 1

2--- x+ 2x–

5.1.2 Elementi triangolari

Nel caso di elementi triangolari le funzioni interpolatrici sono definite in unsistema 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 verticiopposti.

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 totaledel triangolo, A.

5.14L1

A1

A------= L2

A2

A------= L3

A3

A------=

Fig. 5.7 – Elemento triangolare a 3 nodi.

139

I tre rapporti L1, L2, L3 rappresentano le coordinate naturali o d'area per untriangolo (fig. 5.8). Chiaramente deve essere:

5.15L1 L2 L3+ + 1=


Fig. 5.8 – Elemento triangolare a 3 nodi: sistema di riferimento naturale.

140

Come già visto per il caso monodimensionale, le coordinate (x, y) di un punto Ppotranno essere descritte da:

5.16

Effettuando un cambio di variabile:

5.17

le coordinate x, y divengono:

5.18

5.19

avendo definito le funzioni di forma n1, n2 e n3:

5.20

La figura 5.9 mostra l'andamento della funzione di forma n1.

x L1x1 L2x2 L3x3+ +=

y L1y1 L2 y2 L3 y3+ +=

L1 1 x– h–= L2 x= L3 h=

x 1 x– h–( )x1 x x2 h x3+ +=

y 1 x– h–( )y1 x y2 h y3+ +=

x n1x1 n2x2 n3x3+ +=

y n1 y1 n2 y2 n3 y3+ +=

n1 1 x– h–= n2 x= n3 h=


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:

5.21

n1 L1 2L1 1–( ) 1 3 x h+( ) 2 x h+( )2+–= =

n2 L2 2L2 1–( ) x 2x 1–( )= =

n3 L3 2L3 1–( ) h 2h 1–( )= =

n4 4L1L2 4x 1 x– h–( )= =

n5 4L2L3 4xh= =

n6 4L1L3 4h 1 x– h–( )= =

Fig. 5.10 – Elemento triangolare a 6 nodi.

La figura 5.11 mostra l'andamento delle funzioni di forma relative rispettivamente adun 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


Le tabelle 5.3 e 5.4 illustrano le funzioni di forma e le loro derivate rispettiva-mente 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

Tab. 5.4 – funzioni di forma e loro derivate per l'elemento triangolare a 6 nodi

NODO 1 2 3

NODO 1 2 3 4 5 6

n 1 x– h– x h

∂n ∂x§ 1– 1 0

∂n ∂h§ 1– 0 1

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 quadran-golare. L'origine del sistema di riferimento naturale è nel baricentro geometricodell'elemento, mentre gli assi locali x ed � intersecano i lati del quadrilatero nellaloro mezzeria.

Fig. 5.12 – Elemento quadrangolare a 4 nodi.

142


Le funzioni di forma per un elemento quadrangolare a quattro nodi sono date da:

5.22

n114--- 1 x–( ) 1 h–( )= n2

14--- 1 x+( ) 1 h–( )=

n314--- 1 x+( ) 1 h+( )= n4

14--- 1 x–( ) 1 h+( )=

143

Esse sono il risultato di una interpolazione bilineare, che si può descrivere comesegue. 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:

5.23

Sul lato 3-4, allo stesso modo:

5.24

Considerando allora una data ascissa , e interpolando linearmente la funzione secondo l'ordinata :

5.25

Sostituendo le 5.23 e 5.24 nella 5.25 e definendo un unico vettore dei valorinodali fi, si ha:

5.26

5.27

La figura 5.13 mostra l'andamento della funzione di forma n1.

f121 x–

2----------- 1 x+

2------------

f1

f2Ó þÌ ýÏ ¸

=

f341 x+

2------------ 1 x–

2-----------

f3

f4Ó þÌ ýÏ ¸

=

cf h

f 1 h–2

------------ 1 h+2

-------------f12

f34Ó þÌ ýÏ ¸

=

f 1 h–2

------------ 1 h+

2-------------

1 x–2

----------- 1 x+2

------------ 0 0

0 01 x+

2------------ 1 x–

2-----------

f1

f2

f3

f4Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=

f 1 x–( ) 1 h–( )4

---------------------------------- 1 x+( ) 1 h–( )4

----------------------------------- 1 x+( ) 1 h+( )4

----------------------------------- 1 x–( ) 1 h+( )4

-----------------------------------

f1f2f3f4Ó þ

Ô ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=


Fig. 5.13 – Elemento quadrangolare: funzione di forma .n1

Fig. 5.14 – Elementi quadrangolari parabolici.

144

La figura 5.14 illustra due elementi quadrangolari a lati curvi, entrambi con trenodi per lato; il primo elemento è detto elemento Serendipity; il secondo, aventeun nodo anche nel baricentro geometrico, è detto elemento Lagrange.

Nel secondo caso si utilizza un polinomio completo di secondo grado, mentrenel primo viene a mancare il termine quadratico misto x 2y2.

Le funzioni di forma dei due elementi possono essere ricavate mediante sovrap-posizione delle singole componenti, come illustrato in figura 5.15.

La funzione di forma relativa al nodo 5 (fig. 5.15b) può essere ottenuta comeprodotto di una parte quadratica (1 – x2) nella direzione del lato e di una partelineare (1 – �)/2 nella direzione perpendicolare al lato:

5.28

In modo analogo si ricava la funzione di forma n8 relativa al nodo 8 (fig. 5.15c):

5.29

n512--- 1 x2–( ) 1 h–( )=

n612--- 1 x–( ) 1 h2–( )=


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 essereottenuta sottraendo alla funzione bilineare (fig. 5.15a), che fornisce valori pari a0.5 in corrispondenza dei nodi 5 e 8, metà delle funzioni relative ai nodi 5 e 8,avendosi:

5.30

La figura 5.15d mostra la funzione di forma relativa al nodo 1 per un elemento a8 nodi di tipo Serendipity. In modo analogo si ottengono le altre funzioni diforma relative ai nodi d'angolo. Questo modo di procedere suggerisce inoltre lapossibilità di ottenere elementi con numero di nodi variabile per lato, come adesempio indicato in figura 5.16.

n114--- 1 x–( ) 1 h–( ) 1

2---n5–

12---n8–

14--- 1 x–( ) 1 h–( ) x– h– 1–( )[ ]= =

Fig. 5.16 – Elemento quadrangolare con numero di nodi variabile.

145


La tabella 5.5 riassume le funzioni di forma per elementi con numero di nodivariabile tra 4 e 8. Le tabelle 5.6 e 5.7 riassumono le derivate delle funzioni diforma 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) + (SE C'È 6) + (SE C'È 7) + (SE C'È 8)

1

2

3

4

5

6

7

8

14--- 1 x–( ) 1 h–( ) 1

2---n5–

12---n8–

14--- 1 x+( ) 1 h–( ) 1

2---n5– 1

2---n6–

14--- 1 x+( ) 1 h+( ) 1

2---n6–

12---n7–

14--- 1 x–( ) 1 h+( ) 1

2---n7– 1

2---n8–

14--- 1 x2–( ) 1 h–( )

14--- 1 x+( ) 1 h2–( )

14--- 1 x2–( ) 1 h+( )

14--- 1 x–( ) 1 h2–( )

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


Tab. 5.7 – Derivate delle funzioni di forma per l'elemento quadrangolare a 8 nodi

NODO 1 2 3 4

NODO 5 6 7 8

∂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

-------------------------------------

∂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:

5.31L1

V1

V------= L2

V2

V------= L3

V3

V------= L4

V4

V------=

Fig. 5.17 – Elemento tetraedro lineare.

147

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.


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

148

La tabella 5.9 riassume le funzioni di forma e le loro derivate per l'elementotetraedro quadratico a 10 nodi (fig. 5.18).

Tab. 5.9 – Funzioni di forma e loro derivate per l'elemento triangolare a 10 nodi

NODO

1

2

3

4

5

6

7

8

9

10

n ∂n ∂x§ ∂n ∂h§ ∂h ∂z§

l 2l 1–( ) 1 4l– 1 4l– 1 4l–

x 2x 1–( ) 4x 1– 0 0

h 2h 1–( ) 0 4h 1– 0

z 2z 1–( ) 0 0 4z 1–

4xl 4 l x–( ) 4x– 4x–

4xh 4h 4x 0

4hl 4h– 4 l h–( ) 4h–

4zl 4z– 4z– 4 l z–( )

4xz 4z 0 4x

4hl 0 4z 4h

l 1 x– h– z–=


Fig. 5.18 – Elemento tetraedro parabolico.

5.1.5 Elementi solidi parallelepipedi

Le funzioni di forma per gli elementi solidi parallelepipedi possono essere rica-vate a partire da quelle per gli elementi quadrangolari piani aggiungendo la terzacoordinata �.

Fig. 5.19 – Elemento solido (brick).

149

In modo compatto le funzioni di forma e le loro derivate per l'elemento lineare a8 nodi (fig. 5.19) possono essere scritte come:


5.32

dove i valori di sono riportati in tabella (tab. 5.10).

ni18--- 1 xxi+( ) 1 hhi+( ) 1 zzi+( )=

x∂∂ni 1

8---xi 1 hhi+( ) 1 zzi+( )=

h∂∂ni 1

8---h

i1 xxi+( ) 1 zzi+( )=

z∂∂ni 1

8---zi 1 xxi+( ) 1 zzi+( )=

xi hi zi,,

Fig. 5.20 – Elemento solido a 20 nodi.

Tab. 5.10 – Valori di per i nodi d'angolo

NODO 1 2 3 4 5 6 7 8

x i h i z i,,

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

150

Per l'elemento quadratico a 20 nodi (fig. 5.20) le funzioni di forma e le loro deri-vate sono:

– per i nodi sui vertici dell'elemento (1, 2, 3, 4, 5, 6, 7, 8), (tab. 5.10):


151

5.33

– per i nodi con x = 0 (9, 11, 13, 15):

5.34

– per i nodi con � = 0 (10, 12, 14, 16):

NODO 9 11 13 15

NODO 10 12 14 16

ni18--- 1 xxi+( ) 1 hhi+( ) 1 zzi+( ) 2xxi– hhi zzi+ +( )=

x∂∂ni 1

8---xi 1 hhi+( ) 1 zzi+( ) 1– 2xxi hhi zzi+ + +( )=

h∂∂ni 1

8---h

i1 xxi+( ) 1 zzi+( ) 1– xxi 2hhi zzi+ + +( )=

z∂∂ni 1

8---zi 1 xxi+( ) 1 zzi+( ) 1– xxi hhi 2zzi+ + +( )=

h 1– 1 1– 1

z 1– 1– 1 1

ni14--- 1 x2–( ) 1 hhi+( ) 1 zzi+( )=

x∂∂ni 1

2---– x 1 hhi+( ) 1 zzi+( )=

h∂∂ni 1

4---hi 1 x2–( ) 1 zzi+( )=

z∂∂ni 1

4---zi 1 x2–( ) 1 hhi+( )=

x 1 1– 1 1–

z 1– 1– 1 1


152

5.35

– per i nodi con � = 0 (17, 18, 19, 20):

5.36

5.2 CALCOLO DELLA MATRICE DI RIGIDEZZA

Le leggi dello spostamento sono definite, così come per le coordinate, nel sistemadi riferimento naturale; nei casi mono-, bi- e tri-dimensionale si ha, rispettiva-mente:

NODO 1 18 19 20

ni14--- 1 xxi+( ) 1 h2–( ) 1 zzi+( )=

x∂∂ni 1

4---xi 1 h2–( ) 1 zzi+( )=

h∂∂ni 1

2---– h 1 xxi+( ) 1 zzi+( )=

z∂∂ni 1

4---zi 1 xxi+( ) 1 h2–( )=

x 1 1 1– 1–

h 1– 1 1 1–

ni14--- 1 xxi+( ) 1 hhi+( ) 1 z2–( )=

x∂∂ni 1

4---xi 1 hhi+( ) 1 z2–( )=

h∂∂ni 1

4---hi 1 xxi+( ) 1 z2–( )=

z∂∂ni 1

2---– z 1 xxi+( ) 1 hhi+( )=


5.37

dove ui , vi , wi (i = 1,m) sono gli spostamenti nodali dell'elemento, ni sono lefunzioni di forma definite nel sistema di coordinate naturali (�, �, � ) dell'ele-mento.

È da notare che le leggi per definire gli spostamenti possono essere diverse daquelle utilizzate per definire la geometria dell'elemento; nelle espressioni deglispostamenti la sommatoria è estesa all'indice m, mentre nelle espressioni dellageometria la sommatoria è estesa all'indice n.

È tuttavia utilizzata la formulazione isoparametrica, nella quale è utilizzato lostesso numero di parametri (n = m) per definire sia la geometria sia gli sposta-menti.

Altri tipi di formulazione danno luogo a elementi subparametrici e superparame-trici; nei primi la legge per definire la geometria ha un grado minore di quellautilizzata per definire gli spostamenti (n < m), nei secondi ha un grado maggiore(n > m).

Riferendosi per semplicità al caso bidimensionale si avrà quindi:

u ni x( )uii 1=

m

Â=

u ni x h,( )uii 1=

m

Â=

v ni x h,( )vii 1=

m

Â=

u ni x h z,,( )uii 1=

m

Â=

v ni x h z,,( )vii 1=

m

Â=

w ni x h z,,( )wii 1=

m

Â=

5.38

u ni x h,( )uii 1=

m

Â=

v ni x h,( )vii 1=

m

Â=

153

Per il calcolo della matrice di rigidezza [k] dell'elemento:

5.39

è necessario calcolare la matrice [b] che lega le deformazioni agli spostamenti:

5.40

avendo definito il vettore {e} come:

5.41

k[ ] b[ ]T E[ ] b[ ] VdVÚ b[ ]T E[ ] b[ ]h Ad

AÚ= =

e{ } b[ ] s{ }=

e{ }T exx eyy g xy{ }=


154

con:

5.42

Sostituendo agli spostamenti u e v le espressioni delle funzioni di forma si ha:

5.43

e in forma matriciale, esplicitando il solo termine i-esimo:

5.44

Sono quindi da calcolare le derivate delle funzioni di forma ni rispetto alle varia-bili x, y , mentre dette funzioni dipendono esplicitamente da x, �. Per la genericafunzione si può scrivere:

5.45

e in forma matriciale:

exx x∂∂u= eyy y∂

∂v= g xy y∂∂u

x∂∂v+=

exx x∂∂niui

i 1=

n

Â=

eyy y∂∂nivi

i 1=

n

Â=

g xy y∂∂niui x∂

∂nivi+Ë ¯Ê �

i 1=

n

Â=

e{ }

º º x∂

∂ni 0 º º

º º 0 y∂

∂ni º º

º º y∂

∂ni x∂

∂ni º º

ººui

vi

ººÓ þ

Ô ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

x∂∂ni

x∂∂ni

x∂∂x

y∂∂ni

x∂∂y

+=

h∂∂ni

x∂∂ni

h∂∂x

y∂∂ni

h∂∂y

+=


155

5.46

5.47

dove [ J ] è detta matrice jacobiana e mette in relazione le derivate nelle coordi-nate naturali a quelle nel piano fisico:

5.48

e i coefficienti della matrice sono dati da:

5.49

Ammesso che l'inversa della matrice jacobiana esista, cioè che si abbia una rela-zione biunivoca tra le coordinate naturali e quelle fisiche dell'elemento, si ha:

5.50

x∂∂ni

h∂∂ni

Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

x∂∂x

x∂∂y

h∂∂x

h∂∂y

x∂∂ni

y∂∂ni


=

x∂∂ni

h∂∂ni


J[ ] x∂∂ni

y∂∂ni


=

J ][ x∂∂x

x∂∂y

h∂∂x

h∂∂y

=

x∂∂x

x∂∂ni xi

i 1=

n

Â=

h∂∂x

h∂∂ni xi

i 1=

n

Â=

x∂∂y

x∂∂ni yi

i 1=

n

Â=

h∂∂y

h∂∂ni yi

i 1=

n

Â=

x∂∂ni

y∂∂ni


J ][ 1– x∂∂ni

h∂∂ni


=


5.51x∂

∂ni

y∂∂ni


1det J ][---------------- h∂

∂y ∂y∂x------ –

∂x∂h------– x∂

∂x

x∂∂ni

h∂∂ni


=

È così possibile valutare i vari coefficienti della matrice [b] e procedere quindi alcalcolo della matrice di rigidezza [k], che non può, in generale, essere fatto utiliz-zando una integrazione analitica. Si deve quindi ricorrere ad uno schema di inte-grazione numerica.

Prima però di eseguire l'integrazione bisogna esprimere anche il termine dA infunzione delle variabili naturali x, �; la figura 5.21 mostra un elemento di areadA, compreso tra le linee di coordinate x e x + dx e � e � + d�. I lati del paralle-logramma così individuato valgono:

5.52

Fig. 5.21 – Elemento di area dA.

a{ } x∂∂x

x∂∂y


d x= b{ }h∂

∂x

h∂∂y


dh=

156

L'area dA è data dal prodotto vettoriale di {a} per {b}:

5.53

e la matrice di rigidezza [k]:

Ad a{ } b{ }∑ x∂∂x

h∂∂y

x∂∂y

h∂∂x

–Ë ¯Ê � x hdd det J ][ x hdd= = =


157

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:

5.55

dove m è l'ordine di integrazione, wi, wj sono i pesi di integrazione, [b]ij edet [ J ]ij sono calcolati nei punti di integrazione xi, �j.

ESEMPIO 5.1

Calcolare le matrici di rigidezza nel caso di elementi monodimensionali di tipolineare e parabolico.

Elemento asta lineareNel caso di un elemento asta isoparametrico a due nodi le espressioni dellageometria e degli spostamenti sono date da:

5.56

La deformazione nell'elemento è data da:

5.57

e:

5.58

dove è la matrice jacobiana della trasformazione. Le derivate delle fun-zioni di forma rispetto a sono:

5.59

essendo inoltre la lunghezza dell'elemento si ha:

5.60

k[ ] b[ ]T E[ ] b[ ]h AdAÚ b[ ]T E[ ] b[ ]det J[ ] x hdd

1–

1

Ú1–

1

Ú= =

k[ ] wiwj b[ ]ijT E[ ] b[ ]ij det J[ ]ij h

j 1=

m

Âi 1=

m

Â=

x1 x–

2-----------x1

1 x+2

------------x2+=

u1 x–

2-----------u1

1 x+2

------------u2+=

exx xddu

xd

dn1

xd

dn2u1

u2Ó þÌ ýÏ ¸

b[ ] s{ }= = =

xd

dni

xddx 1–

xd

dni J[ ] 1– xd

dni= =

J[ ]x

xd

dn1 12---–=

xd

dn2 12---=

l x2 x1–=

xddx

xd

dn1x1 xd

dn2x2+12---x1–

12---x2+ l

2---= = =


per cui:

5.61

e la matrice :

5.62

La matrice di rigidezza, considerando costanti l'area ed il modulo elastico è:

5.63

5.64

Elemento asta parabolicoLo spostamento è espresso da:

J[ ] l2---= J[ ] 1– 2

l---= det J[ ] l

2---=

b[ ]

b[ ] 1l---–

1l---=

AE

k[ ] EA b[ ]T b[ ] xd0

l

Ú EA b[ ]T b[ ] l2--- xd

1–

1

Ú= =

k[ ] EAl

-------1 1–

1– 1=

u

5.65

Le derivate delle funzioni di forma rispetto a sono:

5.66

u 12---x 1 x–( )–

12---x 1 x–( ) 1 x2–( )

u1

u2

u3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=

x

xd

dn1 12---– x+=

xd

dn2 12--- x+=

xd

dn3 2x–=

per cui la matrice jacobiana è:

5.67

funzione lineare di . Se il nodo centrale 3 è in mezzeria dell'elemento lamatrice jacobiana è una costante e si ha:

J[ ] 12---– x+Ë ¯

Ê � 12--- x+Ë ¯

Ê � 2x–

x1

x2

x3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=

c

5.68J[ ] l2---= J[ ] 1– 2

l---= det J[ ] l

2---=

158

per cui la matrice è:

5.69

b[ ]

b[ ] 2l--- 1

2---– x+Ë ¯

Ê � 12--- x+Ë ¯

Ê � 2x–=


Supponendo costanti l'area ed il modulo elastico , si ha:

5.70

che integrata analiticamente da:

A E

k[ ] 2EAl

----------

x 12---–Ë ¯

Ê � 2 x2 1

4--- – 2x x 1

2---–Ë ¯

Ê � –

x2 14--- – x 1

2---+Ë ¯

Ê � 2 2x x 1

2---+Ë ¯

Ê � –

2x x 12---–Ë ¯

Ê � – 2x x 12---+Ë ¯

Ê � – 4x2

xd

1–

1

Ú=

5.71k[ ] 2EAl

----------

76---

16---

43--- –

16---

76---

43--- –

43--- –

43--- –

83---

=

159

ESEMPIO 5.2

Calcolare la matrice di rigidezza per l'elemento piano a tre nodi in stato ditensione piano o deformazione piano.Gli spostamenti e sono dati da:

5.72

le derivate delle funzioni di forma rispetto alle coordinate locali sono:

5.73

e gli elementi della matrice jacobiana:

5.74

5.75

u v

u 1 x– h–( )u1 xu2 hu3+ +=

v 1 x– h–( )v1 xu2 hv3+ +=

x∂∂n1 1–=

h∂∂n1 1–=

x∂∂n2 1=

h∂∂n2 0=

x∂∂n3 0=

h∂∂n3 1=

x∂∂x

x∂∂nixi

i 1=

3

Â x2 x1–= =

h∂∂x

h∂∂nixi

i 1=

3

Â x3 x1–= =

x∂∂y

x∂∂ni yi

i 1=

3

Â y2 y1–= =

h∂∂y

h∂∂ni yi

i 1=

3

Â y3 y1–= =

J ][x2 x1–( ) y2 y1–( )

x3 x1–( ) y3 y1–( )=


160

e l'inversa:

5.76

dove:

5.77

essendo l'area dell'elemento triangolare. Le derivate delle funzioni di formarispetto alle coordinate fisiche sono quindi:

5.78

e infine la matrice :

5.79

simile a quella ottenuta nell'esempio 4.6 mediante le coordinate generaliz-zate, tenendo conto del fatto che qui si è ordinato il vettore degli sposta-menti nodali nel seguente modo:

5.80

5.3 CARICHI NODALI EQUIVALENTI

5.3.1 Carichi di linea

È qui ricavata l'espressione per il carico nodale equivalente a carichi distribuitilungo un contorno nel caso di elementi piani isoparametrici. Si considerano latiaventi 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, 2nel caso parabolico individuano il contorno su cui agiscono i carichi per unità dilunghezza e stabiliscono il verso della coordinata naturale x sulla linea.

In entrambi i casi illustrati in figura 5.22 sono definiti un carico distribuito tan-

J ][ 1– 1det J ][---------------

y3 y1–( ) y2 y1–( )–

x3 x1–( )– x2 x1–( )=

det J ][ x2 x1–( ) y3 y1–( ) x3 x1–( ) y2 y1–( )– 2A= =

A

x∂∂ni

y∂∂ni


12A------

yi 1+ yi 2+–

xi 2+ xi 1+–Ó þÌ ýÏ ¸

= i 1 2 3, ,=

b[ ]

b[ ]

12A------

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–( )

=

s{ }

s{ }T u1 v1 u2 v2 u3 v3{ }=


genziale qt , positivo se diretto nel verso individuato dai nodi 1 e 2, e un caricodistribuito normale qn , positivo se agisce nel senso che si ottiene ruotando insenso antiorario il vettore tangente alla linea.

Fig. 5.22 – Carichi di linea per elementi lineari e parabolici.

161

Le coordinate di un punto qualsiasi sulla linea sono definite da:

5.81

dove n è il numero di nodi lungo la linea e le ni sono le funzioni di forma defi-nite che per un contorno definito da due nodi sono:

5.82

mentre per un lato definito da tre nodi:

5.83

I carichi distribuiti tangenziale e normale al contorno possono essere espressi infunzione dei valori che tali carichi assumono in corrispondenza dei nodi, in unmodo del tutto analogo alla definizione della geometria o della cinematica:

5.84

Si nota che per lati definiti da due nodi è possibile considerare al massimo carichidistribuiti variabili linearmente, mentre per contorni definiti da tre nodi si potràavere al massimo una legge parabolica.

Il vettore del carico nodale equivalente è:

x ni xii 1=

n

Â= y ni yii 1=

n

Â=

n11 x–

2-----------= n2

1 x+2

------------=

n112---– x 1 x–( )= n2

12---x 1 x+( )= n3 1 x2–=

qt ni qti

i 1=

n

Â= qn ni qni

i 1=

n

Â=


5.85

dove {dq}, vettore delle forze elementari agenti secondo gli assi coordinati x e y ,è dato da:

5.86

Le componenti dqx, dqy possono essere ricavate considerando un tratto dl dilinea (fig. 5.23):

5.87

fe{ }q n[ ]T q{ } ldlÚ n[ ]T qd{ }

lÚ= =

qd{ }T qxd qyd{ }=

qxd qt l acosd qn l asind–=

qyd qt l asind qn l acosd+=

Fig. 5.23 – Componenti del carico di linea.

162

Osservando che:

5.88

si ottiene:

5.89

I termini dx, dy in funzione della coordinata naturale x sono:

5.90

per cui si ha:

l acosd xd= l asind yd=

qxd qt xd qn yd–=

qyd qt yd qn xd+=

xd x∂∂x xd= yd x∂

∂y xd=


163

5.91

Le componenti del carico nodale equivalente al nodo i (i = 1,2 oi = 1,2,3) sono quindi calcolate come:

5.92

dove l'integrale, esteso alla lunghezza del contorno, è generalmente valutatonumericamente.

Nel caso di un contorno definito da due nodi è possibile determinare in modoabbastanza semplice le componenti del vettore del carico nodale equivalente informa esplicita; a titolo di esempio si determinerà la componente :

5.93

dove:

5.94

5.95

5.96

Sostituendo ed integrando per x da –1 a +1, si ha:

5.97

qxd qt x∂∂x

qn x∂∂y

–Ë ¯Ê � xd=

qyd qt x∂∂y

qn x∂∂x

+Ë ¯Ê � xd=

fxifyi

,

fxini qt

∂x ∂x--------- qn

∂y∂x------–Ë ¯

Ê � xd1–

1

Ú=

fyini qt

∂y∂x------ qn

∂x∂x------+Ë ¯

Ê � xd1–

1

Ú=

fx1

fx1n1 qt

∂x∂x------ qn

∂y∂x------–Ë ¯

Ê � xd1–

1

Ú=

n112--- 1 x–( )=

qt12--- 1 x–( )qt1

12--- 1 x+( )qt2

+=

qn12--- 1 x–( )qn1

12--- 1 x+( )qn2

+=

x∂∂x

x∂∂n1x1 x∂

∂n2x2+x2 x1–

2----------------= =

x∂∂y

x∂∂n1y1 x∂

∂n2y2+y2 y1–

2---------------= =

fx1x2 x1–( ) 1

3---qt1

16---qt2

+Ë ¯Ê � y2 y1–( ) 1

3---qn1

16---qn2

+Ë ¯Ê �+=


164

5.3.2 Carichi di volume

Il carico nodale equivalente a carichi di volume è dato da:

5.98

dove {�} nel caso piano è:

5.99

e dV = h det [ J ] d x d �. Inoltre nel caso di peso proprio con l'asse Y orientatonel senso della gravità si ha:

5.100

Il calcolo del vettore dei carichi nodali equivalenti segue le stesse procedure visteper il calcolo delle matrici di rigidezza; ad esempio nel caso piano è calcolabilecome:

5.101

5.3.3 Effetti termici

Il carico nodale equivalente ad effetti termici viene calcolato come:

5.102

dove nel caso di stato di tensione piano è:

5.103

e, nel caso di stato di deformazione piano:

5.104

in cui è il coefficiente di dilatazione termica del materiale, è la tempera-tura di riferimento e è la temperatura in un punto generico, che può essereespressa come:

5.105

dove è il numero di nodi dell'elemento, sono le funzioni di forma associateall'elemento e è il valore della temperatura al nodo i dell'elemento.

fe{ }f n[ ]T f{ } VdVÚ=

f{ }T fx fy{ }=

fx 0= fy rg=

fe{ }f wiwj n[ ]ijT f{ }det J ][ ij h

j 1=

m

Âi 1=

m

Â=

fe{ }e 0b[ ]T E[ ] e0{ } Vd

VÚ=

e0{ }

e0{ }T a* T T0–( ) 1 1 0{ }=

e0{ }T a* T T0–( ) 1 n+ 1 n+ 0{ }=

a* T0T

T niTii 1=

n

Â=

n niTi


Si ha quindi:

5.106

5.4 CALCOLO DELLE TENSIONI

Noto il vettore degli spostamenti della struttura gli spostamenti relativi ai singoli elementi possono essere ricavati con una operazione inversa aquella utilizzata per assemblare gli elementi.

Le deformazioni nell'elemento sono quindi date da:

5.107

e le tensioni:

5.108

Lo stato di tensione dipende dal punto che si considera all'interno dell'elementoperché in generale la matrice è funzione delle coordinate; questa devequindi 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 nodidell’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 menopreciso, 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 indeterminati punti interni dell'elemento; questi punti possono essere utilizzati percaratterizzare il campo di tensione che poi potrà essere estrapolato per valutare letensioni nei nodi dell'elemento.

La determinazione dei punti ottimali per il calcolo delle tensioni può essere fattasecondo quanto suggerito da Barlow1; l’obiettivo è quello di determinare i puntinei quali le tensioni sono calcolate con lo stesso grado di precisione con il qualesono calcolati gli spostamenti nodali.

Si assume per il dato elemento un campo di spostamenti { } di un ordine supe-riore a quello comunemente utilizzato e si cercano i punti nei quali lo statodi tensione calcolato a partire dai due campi di spostamento è uguale

. In questi punti lo stato di tensione avrà lo stesso grado di

fe{ }e0wiwj b[ ]ij

T E[ ] e0{ }ij det J[ ]ij hj 1=

m

Âi 1=

m

Â=

S{ } s{ }h

e{ } b[ ] s{ }=

s{ } E[ ] b[ ] s{ } e0{ }–( ) s0{ }+=

b[ ]

u 'u{ }

s '{ } s{ }= s{ }

1.Barlow, J., Optimal stress location in finite element models, International Journal forNumerical Methods In Engineering, 1976, vol. 10, p. 243-251.

165


166

precisione dello spostamento ; detti punti coincidono con i punti di Gaussdi sottointegrazione e sono talvolta chiamati punti di Barlow.

Nel caso monodimensionale della trave inflessa il procedimento seguito da Bar-low per la determinazione dei punti ottimali può essere cos riassunto. Il campodi spostamenti può essere espresso come:

5.109

dove x è una coordinata adimensionale misurata dal centro della trave e variabileda –1 a 1 e {a} è il vettore delle coordinate generalizzate.

Gli spostamenti nodali sono dati da:

5.110

Assumendo una funzione di ordine superiore, lo spostamento della trave saràdato da:

5.111

e gli spostamenti nodali sono dati da:

5.112

Se l’elemento è capace di rappresentare correttamente il campo di spostamentinei due casi considerati, allora gli spostamenti nodali ottenuti nei due casisaranno approssimativamente uguali:

5.113

5.114

con:

5.115

Le deformazioni flessionali e quindi le tensioni flessionali saranno uguali nei duecasi quando:

5.116

s{ }

va 1 x x2 x3 a{ }=

sa{ }

sa{ } A[ ] a{ }=

vb 1 x x2 x3 x4 b{ }=

sb{ }

sb{ } B[ ] b{ }=

A[ ] a{ } B[ ] b{ }ª

a{ } A[ ] 1– B[ ] b{ }=

A[ ] 1– B[ ]

1 0 0 0 1 –

0 1 0 0 0 0 0 1 0 2 0 0 0 1 0

=

x2

2

∂∂ va

x2

2

∂∂ vb=


167

ovvero quando:

5.117

e sostituendo in termini di {b}:

5.118

Questa uguaglianza è soddisfatta per tutti i valori di {b} solo se:

5.119

La stessa tecnica può essere applicata all'elemento piano isoparametrico a quattronodi; si trova che i punti ottimali per il calcolo delle tensioni sono:

5.120

mentre per l'elemento piano isoparametrico a otto nodi si ha:

5.121

Le tensioni ai nodi dell'elemento possono quindi essere ottenute estrapolando ivalori calcolati nei punti di Gauss (altrimenti detti di Barlow). Nel caso monodi-mensionale (trave inflessa o elemento asta a tre nodi) si ottiene:

5.122

mentre nel caso bidimensionale (elemento piano isoparametrico a 8 nodi) siottiene:

5.123

e nel caso tridimensionale (elemento solido isoparametrico a 20 nodi) si ottiene:

0 0 2 6 x ][ a{ } 0 0 2 6 x 12 x2 b{ }=

0 0 2 6 x 4 b{ } 0 0 2 6x 12x2 b{ }=

x 1

3-------±=

x 0= h 0=

x 1

3-------±= h 1

3-------±=

s1

s2Ó þÌ ýÏ ¸

1 3+2

---------------- 1 3–2

----------------

1 3–2

---------------- 1 3+2

----------------

sI

sIIÓ þÌ ýÏ ¸

=

s1

s2

s3

s4Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

13

2------- +

12--- – 1

32

------- –12--- –

12--- – 1

32

------- +12--- – 1

32

------- –

13

2------- –

12--- – 1

32

------- +12--- –

12--- – 1

32

------- –12--- – 1

32

------- +

sI

sII

sIII

sIVÓ þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=


168

5.124

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

Le tensioni così calcolate ai nodi di un elemento saranno in generale diverse dalletensioni calcolate negli stessi nodi appartenenti però ad elementi adiacenti; ciòrisulta dal fatto che le tensioni ai nodi dipendono dagli spostamenti di tutti inodi dell'elemento.

Ciò deriva anche dal fatto che nell’equazione dei lavori virtuali è garantito l’equi-librio globale (equilibrio delle forze nodali), ma non quello locale punto perpunto (equilibrio delle tensioni).

Ad un nodo allora sono attribuiti diversi valori di tensione, ciascuno derivantedai vari elementi che concorrono in quel nodo; è usuale calcolare, a partire da dativalori, una tensione media nodale ed attribuire al nodo quest'unico valore di ten-sione media, valore che sarà probabilmente più esatto dei singoli valori calcolabili.

ELEMENTO NO PUNTI COORDINATE NATURALI

TRAVE INFLESSA

ISOPARAMETRICO PIANO A 4 NODI

ISOPARAMETRICO PIANO A 8 NODI

ISOPARAMETRICO SOLIDO A 8 NODI

ISOPARAMETRICO SOLIDO A 20 NODI

s1

s2

s3

s4

s5

s6

s7

s8Ó þÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

a3 a2b ab2 a2b a2b ab2 b3 ab2

a2b a3 a2b ab2 ab2 a2b ab2 b3

ab2 a2b a3 a2b b3 ab2 a2b ab2

a2b ab2 a2b a3 ab2 b3 ab2 a2b

a2b ab2 b3 ab2 a3 a2b ab2 a2b

ab2 a2b ab2 b3 a2b a3 a2b ab2

b3 ab2 a2b ab2 ab2 a2b a3 a2b

ab2 b3 ab2 a2b a2b ab2 a2b a3

sI

sII

sIII

sIV

sV

sVI

sVII

sVIIIÓ þÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

a 1 3+= b 1 3–=

2 x a±=

1 x 0;= h 0=

2 2¥ x a;±= h a±=

1 x 0;= h 0=

2 2¥ 2¥ x a;±= h a;±= z a±=

a 1 3§±=


169

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 inte-granda; nel caso della matrice di rigidezza di un elemento piano l'ordine di inte-grazione si può valutare determinando il grado massimo della funzione:

5.125

tenendo conto che un polinomio di grado 2m – 1 è integrato esattamente da mpunti di Gauss.

Per elementi di forma regolare è abbastanza semplice determinare il grado dellafunzione integranda, perché det [ J ] è costante; non altrettanto si può dire perelementi di forma complessa, dove il determinante della matrice jacobiana variain modo irregolare. Di conseguenza l’ordine di integrazione viene anche stabilitoin base a verifiche numeriche sull'elemento.

Una integrazione esatta della matrice di rigidezza, in special modo per elementidi grado elevato, può portare ad avere lunghi tempi di calcolo, per cui è a voltepreferibile una integrazione ridotta o sottointegrazione; si utilizzano cioè menopunti di campionamento di quelli necessari per integrare esattamente la funzioneintegranda.

A tal proposito si ricorda ora che con la formulazione degli elementi finiti a spo-stamenti assegnati si ha sempre una stima in eccesso della rigidezza del sistema(avendo limitato la cinematica del sistema a quella espressa dalle funzioni diforma). Quindi sottointegrare, cioè sottovalutare la rigidezza dell'elemento, può,di fatto, portare a risultati migliori a patto che l'errore nell'integrazione nume-rica compensi la sovrastima della rigidezza strutturale dovuta al tipo di formula-zione.

Esiste però un limite alla sottointegrazione perché è necessario poter calcolareesattamente il volume dell’elemento. Questa condizione deriva dal fatto che allimite, per suddivisioni molto fitte, nell'elemento si verifica uno stato di defor-mazione costante (cosa che l'elemento deve in ogni caso essere in grado di rap-presentare) e quindi l'energia di deformazione può, in tali condizioni assumere laforma, limitandosi al caso piano:

5.126

Quindi si deve integrare esattamente il prodotto det [ J ] h; nel caso di spessorecostante e di elemento lineare a quattro nodi la funzione integranda è una fun-zione lineare in x e � e quindi è necessario un solo punto di campionamento; nelcaso di elemento a otto nodi a spessore costante si ha invece una dipendenzacubica da x e � e quindi è necessario uno schema di integrazione 2 x 2.

f b[ ]T E[ ] b[ ]det J[ ]h=

U k*[ ] h det J[ ] x hddVÚµ


170

In pratica però una sottointegrazione troppo spinta può portare a matrici singo-lari 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 modidi deformazione con energia di deformazione nulla.

Per un singolo elemento il numero di meccanismi può essere calcolato mediantela seguente relazione:

5.127

dove G è il numero di gradi di libertà per nodo, N il numero di nodi dell'ele-mento, L il numero di gradi di libertà di moto rigido, r il rango della matrice dielasticità [E ] e n il numero di punti di campionamento. Se M � � non si hannomeccanismi.

Si può pensare infatti di intendere l'integrazione numerica nel seguente modo:Integrare numericamente significa campionare la funzione integranda in alcunipunti nel dominio dell'elemento. Nel nostro caso significa valutare la rigidezzadell'elemento in alcuni punti ed attribuire questa alla rigidezza reale.

Ciò equivale, fisicamente, a sostituire una rigidezza distribuita sull'elemento condelle molle spia concentrate nei punti di integrazione.

Il numero delle molle da considerare per ogni punto dipende dal numero dellecomponenti del vettore deformazione (o tensione) e cioè dal rango della matrice[E ] (ad esempio nel caso piano si può considerare una rigidezza associata allacomponente , una alla ed una terza associata alla .

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 dilibertà di deformazione dell'elemento, cioè:

5.128

Ad esempio nel caso di problemi bidimensionali in stato di tensione piano odeformazione piano il numero di gradi di libertà per nodo è G = 2, il rango dellamatrice [E] è r = 3, ed il numero di gradi di libertà di moto rigido è L = 3 per cuisi ha:

5.129

Per un elemento a 4 nodi, bilineare, sottointegrato con un punto di campiona-mento (n = 1) si hanno due meccanismi:

5.130

mentre nel caso di un elemento parabolico a otto nodi (Serendipity) con una sot-tointegrazione con quattro punti di campionamento (2 x 2) si ha un meccani-smo nascosto.

In quest’ultimo caso però se si hanno più elementi adiacenti il meccanismo non

M GN L– rn–=

exx eyy g xy

rn GN L–≥

M 2N 3– 3n–=

M 2 4 3– 3 1◊–◊ 2= =


171

si propaga perché gli spostamenti corrispondenti ad energia nulla in un elementosono incompatibili con gli spostamenti di elementi adiacenti per cui un insiemedi elementi può non avere modi di deformazione ad energia nulla, anche se que-sti sono presenti in un elemento, così che una integrazione 2 x 2 può essere uti-lizzata efficacemente.

La tabella 5.12 riassume gli ordini di integrazione per elementi piani, distin-guendo tra ordini di integrazione massimi e ordini di integrazione ridotti.

Tab. 5.12 – Ordine di integrazione per elementi quadrangolari

Tuttavia nel caso di sottointegrazione è consigliabile effettuare verifiche numeri-che dell'elemento, quali il calcolo degli autovalori della matrice di rigidezza ed ilpatch test.

5.5.2 Distorsione degli elementi

Gli elementi forniscono i risultati migliori quando la loro forma geometricacoincide con quella dell'elemento di riferimento, come triangoli equilateri e qua-drati nel piano, o cubi nello spazio.

In pratica comunque è quasi impossibile che tutti gli elementi abbiano unaforma geometrica ideale; si parla allora di elementi distorti (rispetto alla forma diriferimento). Un limitato valore di distorsione può essere comunque accettato.

Oltre un determinato limite di distorsione, però, la rigidezza dell'elemento vienevalutata con scarsa precisione; oltre un limite critico diventa addirittura impossi-bile 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 categoriefondamentali: il rapporto lunghezza altezza, il rapporto di forma e la distorsione

ELEMENTO ORDINE DI INTEGRAZIONE MASSIMO

ORDINE DI INTEGRAZIONE RIDOTTO

4 NODI RETTANGOLO

4 NODI DISTORTO

8 NODI RETTANGOLO

8 NODI DISTORTO

2 2¥ 2 2¥

2 2¥ 2 2¥

3 3¥ 2 2¥

3 3¥ 3 3¥


angolare (fig. 5.24).

Non esiste un modo generale per analizzare questi tipi di distorsioni; è necessarioverificare numericamente l'elemento in diverse condizioni di distorsione e dicarichi.

Fig. 5.24 – Distorsione di un elemento.

172

Con riguardo al primo tipo di distorsione (fig. 5.24a) si può accettare un rap-porto 1:2 per elementi situati in zone dove si vuole conoscere lo stato di tensionee 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 di1:10 in zone non critiche.

Per ciò che riguarda la distorsione angolare (fig. 5.24c) è in genere consigliabilenon avere angoli interni minori di 45˚ per elementi quadrangolari e minori di15˚ per elementi triangolari.

Una verifica della distorsione può essere fatta anche mediante il calcolo dellamatrice jacobiana e del suo determinante. Il determinante della matrice jaco-biana è infatti costante se l’elemento è un parallelepipedo con i nodi sui latiequamente spaziati, ma varia se esiste altro tipo di distorsione. Esso può al limitediventare nullo o negativo nelle vicinanze di una zona ad eccessiva distorsione,con la conseguenza di impedire il calcolo della matrice di rigidezza dell'ele-mento. Causa di ciò è la mancanza di corrispondenza biunivoca con conseguentenon invertibilità della matrice jacobiana.

Ad esempio nel caso dell’elemento asta parabolico a tre nodi la matrice jacobianaè data da:

5.131J[ ] 12---– x+Ë ¯

Ê � 12--- x+Ë ¯

Ê � 2x–( )

x1

x2

x3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=


Fig. 5.25 – Distorsione di un elemento asta a 3 nodi.

Se il nodo centrale non è in mezzeria dell’elemento, supponendo x1 = 0, x2 = lcon x3 compresa tra 0 e l (fig. 5.25) il determinante dello jacobiano è:

5.132

La posizione x3 del nodo centrale che rende singolare la matrice jacobiana puòessere calcolata annullando il determinante dello jacobiano e tenendo presenteche . Per si ha x3 = l/4 mentre per si ha x3 = 3l/4.

I valori ammissibili di x3 per i quali [ J ] non è singolare sono quindi:

5.133

J 12--- x+Ë ¯

Ê � l 2xx3–=

1– x 1£ £ x 1–= x 1=

l4--- x3

3l4----£ £

Fig. 5.26 – Distorsione di un elemento quadrangolare.

173

Nel caso di elementi quadrangolari a quattro e otto nodi Zienkiewicz ha propo-sto, per non avere una matrice jacobiana singolare, che gli angoli interni sianominori di 180 gradi e per gli elementi parabolici i nodi sui lati debbano esserecompresi nella terza parte centrale (fig. 5.26).

Cook ha proposto il calcolo di un fattore di distorsione DP, dato rispettivamenteper elementi quadrangolari e elementi solidi:


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 presoin 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 è ingrado di rappresentare correttamente uno stato di tensione con presenza di gra-dienti.

Ad esempio nel caso di sollecitazione a flessione pura l’elemento presenta uncomportamento molto più rigido della realtà e questo può essere imputatoall'effetto del taglio, che teoricamente deve essere nullo.

Nel caso di integrazione esatta (n = 2 x 2) infatti si ha, a causa dei vincoli cine-matici imposti dalle funzioni di forma adottate, una deformazione a taglio neipunti di campionamento utilizzati (fig. 5.27a).

Se invece si procede con una sottointegrazione (n = 1) si nota che nel punto dicampionamento (punto centrale) non si ha deformazione ha taglio e quindil'energia di deformazione a taglio è calcolata in modo esatto (fig. 5.27b).

DP 4minJArea-----------Ë ¯

Ê �= DP 8minJ

Volume------------------Ë ¯

Ê �=

Fig. 5.27 – Punti di integrazione per elemento a 4 nodi.

174

Il campo di spostamenti all'interno dell’elemento è definito da:

5.135

ed essendo:

u u*xh= v 0=


175

5.136

si ha:

5.137

Il campo corretto di spostamenti invece prevede:

5.138

Nel caso di sollecitazione derivante da momento flettente puro si deve avere chela deformazione a taglio �xy sia nulla dovunque; ciò è verificato dal campo dispostamenti corretto:

5.139

ma non da quello associato all'elemento quadrangolare:

5.140

La stessa cosa si ottiene considerando un momento flettente applicato lungo ilati di lunghezza l; il campo di spostamenti approssimato risulterà:

5.141

Si può notare che i due modi flessionali forniscono uno stato di deformazioneesatto solo nel punto centrale dell'elemento x = 0; � = 0; questo suggeriscequindi la possibilità di ridurre la sovrastima di rigidezza effettuando una integra-zione selettiva in cui le deformazioni normali sono integrate esattamentemediante una griglia 2 x 2, mentre i termini corrispondenti alla deformazionedovuta al taglio sono sottointegrati, calcolati cioè nel punto centrale dell'ele-mento.

Un modo per ottenere risultati più soddisfacenti è quello di ricorrere ad elementiquadratici (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à perelemento) è però molto gravoso nel caso tridimensionale, dove si passa da 24 a60 gradi di libertà, quasi triplicando l'ordine delle matrici in esame.

È quindi di grande importanza riuscire a migliorare il comportamento flessio-nale degli elementi lineari; tra i diversi procedimenti proposti si hanno il metododell'integrazione selettiva o quello dell'aggiunta di modi incompatibili.

x 2xl

-----= h 2yh-----=

u4lh----u*xy= v 0=

u4lh----u*xy= v

lu*2h------- 1 4x2

l2--------–Ë ¯

Ê � nhu*2l

--------- 1 4y2

h2--------–Ë ¯

Ê �+=

g xy y∂∂u

x∂∂v+

4lh----u*x

4lh----u*x– 0= = =

g xy y∂∂u

x∂∂v+

4lh----u*x= =

v4lh----v*xy=


176

Integrazione selettiva

Nell’integrazione selettiva si utilizzano ordini di integrazione diversi per lediverse forme di energia di deformazione (flessionale e di taglio). Si ha infatti, da:

5.142

separando la deformazione dovuta al taglio dalla deformazione normale:

5.143

5.144

essendo, ad esempio nel caso di stato di tensione piano:

5.145

La matrice di rigidezza dell’elemento può allora essere scritta come somma delcontributo dovuto alle deformazioni normali e di quello dovuto alle deforma-zioni di taglio:

5.146

Modi incompatibili

Un altro modo per migliorare la risposta dell’elemento lineare è quello di correg-gere 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 integra-zione diverso per l'elemento Q6, formulando quindi l’elemento QM6. La diffe-renza tra i due elementi è che il secondo supera il patch test per qualsiasiconfigurazione geometrica, mentre il primo solo se di forma rettangolare. Nelcaso bidimensionale si ha:

k[ ] b[ ]T E[ ] b[ ]h AdAÚ=

b[ ]b1[ ]

b2[ ]=

E[ ]E11[ ] E12[ ]

E21[ ] E22[ ]=

E11[ ] E1 n2–-------------- 1 n

n 1= E22[ ] G[ ]= E12[ ] E21[ ]T 0

0= =

k[ ] e{ }T s{ } VdVÚ g{ }T t{ } V=d

VÚ+=

= b1[ ] E11[ ] b1[ ]h AdAÚ b2[ ] E22[ ] b2[ ]h Ad

AÚ+


177

5.147

mentre nel caso tridimensionale:

5.148

dove i coefficienti ai sono definiti come spostamenti generalizzati e possonoessere pensati come gradi di libertà interni dell'elemento; è da notare che un ele-mento così definito è incompatibile essendo gli ai definiti univocamente per cia-scun elemento. In generale si ha:

5.149

e:

5.150

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:

u ni x h,( )uii 1=

4

Â 1 x2–( )a1 1 h2–( )a2+ +=

v ni x h,( )vii 1=

4

Â 1 x2–( )a3 1 h2–( )a4+ +=

u ni x h z,,( )uii 1=

8

Â 1 x2–( )a1 1 h2–( )a2 1 z2–( )a3+ + +=

v ni x h z,,( )vii 1=

8

Â 1 x2–( )a4 1 h2–( )a5 1 z2–( )a6+ + +=

w ni x h z,,( )wii 1=

8

Â 1 x2–( )a7 1 h2–( )a8 1 z2–( )a9+ + +=

u{ } ns[ ] na[ ]s{ }a{ }Ó þ

Ì ýÏ ¸

=

b[ ] bs[ ] ba[ ]s{ }a{ }Ó þ

Ì ýÏ ¸

=


178

5.151

La matrice di rigidezza avrà ordine 12 x 12 e può essere scritta separando le variesottomatrici, come:

5.152

Poche i gradi di libertà ai sono interni a ciascun elemento essi possono essere eli-minati o condensati a livello elemento; questo procedimento, noto come conden-sazione statica, porta ad ottenere una matrice di rigidezza condensata didimensioni 8 x 8. Dalle equazioni elemento si ha:

5.153

non essendoci forze corrispondenti ai gradi di libertà aggiuntivi. Dal secondosottosistema si ricava:

5.154

da cui, essendo [kaa] non singolare:

5.155

dal primo sottosistema si può quindi ottenere la matrice di rigidezza condensatadell'elemento incompatibile:

5.156

b[ ]

º ºx∂

∂ni 0 º ºx∂

∂n1

x∂∂n2 0 0

º º 0y∂

∂ni º º 0 0y∂

∂n1

y∂∂n2

º ºy∂

∂ni

x∂∂ni º º

y∂∂n1

y∂∂n2

x∂∂n1

x∂∂n2

ººui

vi

ººa1

a2

a3

a4Ó þÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

k[ ]kss[ ] ksa[ ]

kas[ ] kaa[ ]=

kss[ ] ksa[ ]

kas[ ] kaa[ ]

s{ }

a{ }Ó þÌ ýÏ ¸ f{ }

0{ }Ó þÌ ýÏ ¸

=

kas[ ] s{ } kaa[ ] a{ }+ 0{ }=

a{ } kaa[ ] 1– kas[ ] s{ }–=

k[ ] kss[ ] ksa[ ] kaa[ ] 1– kas[ ]–=


L’elemento così definito è l'elemento incompatibile a 4 nodi Q6 introdotto daWilson2. I risultati ottenuti con questo elemento sono nettamente migliori diquelli ottenuti con il semplice elemento a 4 nodi; fu tuttavia notato che nel casodi elementi distorti (non rettangolari) i risultati erano insoddisfacenti. Fu anchenotato che questo elemento superava il patch test nel caso di geometria rettango-lare, ma non nel caso di geometria qualsiasi.

Successivamente Taylor3 modificò la formulazione dell'elemento in modo che ilpatch test fosse superato per qualsiasi configurazione geometrica; questo risultatofu ottenuto semplicemente calcolando il contributo dei modi incompatibiliall'integrale dell'operatore [B] nel centro dell'elemento.

L’elemento così formulato, detto QM6, supera il patch test ed è quello general-mente utilizzato, in quanto fornisce risultati migliori anche dell'elemento conintegrazione selettiva.

L’idea di Taylor fu la seguente: affinché l’elemento passi il patch test (cioè si com-porti come un elemento compatibile) il termini vettore {a } deve essere nulloquando { s } corrisponde ad un moto rigido o ad uno stato di deformazionecostante . Dovrà cioè essere:

5.157

ed essendo [kaa] definita positiva:

5.158

cioè:

5.159

dove è il vettore delle tensioni calcolato in base ai soli spostamenti e, rappresentando questi un moto rigido, è una costante. Condizione necessariaper il superamento del patch test è quindi:

5.160

s{ } s0{ }=

a{ } kaa[ ] 1– kas[ ] s0{ }– 0= =

kas[ ] s{ } 0=

ba[ ]T E[ ] bs[ ] s0{ } VdVÚ ba[ ]T s0{ } Vd

VÚ 0= =

s0{ } s0{ }

ba[ ] VdVÚ 0=

2.Wilson E. L., Taylor R. L., Doherty W. P., Ghaboussi J., Incompatible displacementmodels, 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. 1211-1219.

179


180

In generale questa condizione non è soddisfatta, a meno che l’elemento non siaun parallelogramma. Si può forzare questa condizione valutando numericamentel’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 potenzialeelastica corrispondente ai vari modi di deformazione e quindi possono rappre-sentare una spia del numero di moti rigidi associati all'elemento.

L’energia di deformazione di un elemento può essere espressa come:

5.161

Si può pensare di scrivere le equazioni di rigidezza riferendole ad un sistema defi-nito dagli autovettori della matrice [k] antiche allo spazio fisico degli sposta-menti { s}. Nel sistema:

5.162

il vettore degli spostamenti { s } può essere espresso come combinazione linearedegli autovettori della matrice di rigidezza [k]:

5.163

dove m è il numero di gradi di libertà dell'elemento, {�} è il vettore delle coordi-nate modali o dei fattori di partecipazione modali e [] è la matrice le cuicolonne rappresentano gli autovettori della matrice [k]:

5.164

Se gli autovettori {�} sono normalizzati in modo che:

5.165

dalla equazione precedente, scritta per l'autovalore i-esimo, si ha:

5.166

5.167

e inoltre essendo gli autovettori ortogonali:

5.168

dove è una matrice diagonale contenente gli autovalori .

U12--- s{ }T k[ ] s{ }=

k[ ] s{ } f{ }=

s{ } ai f{ }ii 1=

m

Â F[ ] a{ }= =

k[ ] f{ } l f{ }=

f{ }iT f{ }i 1=

f{ }iT k[ ] f{ }i li f{ }i

T f{ }i=

f{ }iT k[ ] f{ }i li=

F[ ]T k[ ] F[ ] L[ ]=

L[ ] li


181

Il sistema risolutivo può allora essere scritto come:

5.169

e premoltiplicando per []T ambo i membri:

5.170

cioè:

5.171

ed essendo la matrice [] diagonale il sistema può essere scomposto in m equa-zioni disaccoppiate:

5.172

dalle quali si calcolano direttamente gli m fattori di partecipazione modali:

5.173

e quindi gli spostamenti nodali { s}.

Gli autovettori possono quindi essere identificati come i modi di deforma-zione principali dell'elemento e gli autovalori come rigidezze modali relativea tali modi di deformazione. Essi sono inoltre proporzionali all'energia poten-ziale elastica relativa all'i-esimo modo di deformazione:

5.174

Si ricorda che il sistema di equazioni di rigidezza comprende le equazioni diequilibrio del sistema, equazioni a cui corrispondono i gradi di libertà di motorigido. In queste condizioni l'energia potenziale elastica assorbita dal sistema ènulla e quindi il corrispondente autovalore sarà nullo e l'autovettore ad esso asso-ciato rappresenterà uno dei moti rigidi del sistema.

Se il numero di autovalori nulli è minore del numero di gradi di libertà di motorigido del sistema allora il campo di spostamenti assunto non è adeguato allarappresentazione 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 meccani-smi nascosti, cioè modi di deformazione associati ad una energia di deforma-zione nulla.

Se inoltre gli autovalori non variano al variare dell'orientamento dell'elementorispetto al sistema di riferimento globale, il campo degli spostamenti assunto èinvariante.

Infine modi simili di deformazione devono corrispondere ad autovalori multipli.

k[ ] F[ ] a{ } f{ }=

F[ ]T k[ ] F[ ] a{ } F[ ]T f{ }=

L[ ] a{ } F[ ]T f{ }=

liai f{ }iT f{ }=

ai

f{ }iT f{ }

li-----------------------=

f{ }ili

Ui

f{ }iT k[ ] f{ }i 2Ui li= =


ESEMPIO 5.3

Calcolare gli autovalori della matrice di rigidezza dell'elemento asta a duenodi. Assumendo la matrice di rigidezza è data da:

5.175

e gli autovalori della matrice di rigidezza si ricavano da:

5.176

Il polinomio caratteristico è:

5.177

e gli autovalori:

5.178

ed i corrispondenti autovettori normalizzati sono:

5.179

L'elemento asta possiede un grado di libertà di moto rigido, lo spostamentoassiale; di conseguenza c'è un autovalore nullo a cui corrisponde un autovet-tore che rappresenta appunto un moto rigido. Il secondo autovettore rappre-senta invece il modo di deformazione proprio di questo elemento(allungamento/accorciamento) (fig. 5.28).

EA l§ 1=

k[ ] 1 1–

1– 1=

det 1 1–

1– 1l 00 l

– 0=

l2 2l– 0=

l1 0= l2 2=

f{ }1T 1

2------- 1 1{ }= f{ }2

T 1

2------- 1 1–{ }=

Fig. 5.28 – Modi di deformazione dell’elemento asta a 2 nodi.

182


183

ESEMPIO 5.4

Calcolare gli autovalori della matrice di rigidezza dell'elemento asta a trenodi. Assumendo , la matrice di rigidezza, è data da:

5.180

che integrata analiticamente o numericamente in modo completo con ilmetodo di Gauss con due punti di campionamento (si ha infatti ;

) da:

5.181

mentre sottointegrata, con un punto di campionamento, da:

5.182

Nel caso di integrazione analitica o numerica completa il polinomio caratteri-stico della matrice di rigidezza è:

5.183

gli autovalori risultano quindi:

5.184

ed i corrispondenti autovettori normalizzati sono:

5.185

EA l§ 1=

k[ ]

x 12---–Ë ¯

Ê � 2x2 1

4---–Ë ¯

Ê � 2x x 12---–Ë ¯

Ê �–

x2 14---–Ë ¯

Ê � x 12---+Ë ¯

Ê � 22x x 1

2---+Ë ¯

Ê �–

2x x 12---–Ë ¯

Ê �– 2x x 12---+Ë ¯

Ê �– 4x2

xd1–

1

Ú=

2m 1– 2=m 3 2§= 2Þ

k[ ]

73--- 1

3---

83--- –

13--- 7

3--- 8

3--- –

83---– 8

3--- –

163------

=

k[ ]1 1– 0 1– 1 0 0 0 0

=

l l2 10l– 16+( ) 0=

l1 0=

l2 2=

l3 8=

f{ }1T 1

3------- 1 1 1{ }=

f{ }2T 1

2------- 1 1– 0{ }=

f{ }3T 1

6------- 1 1 2–{ }=


Nel caso di sottointegrazione numerica il polinomio caratteristico dellamatrice di rigidezza è:

5.186

gli autovalori risultano quindi:

5.187

i primi due autovettori, corrispondenti agli autovalori nulli, hanno la terzacomponente, corrispondente al nodo centrale, indeterminata, mentre le duecomponenti, corrispondenti ai nodi di estremità, sono uguali; il terzo auto-vettore normalizzato è:

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 dideformazione 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 trenodi, in cui la matrice di rigidezza è integrata esattamente:

l3 2l2– 0=

l1 0= l2 0= l3 2=

f{ }3T 1

2------- 1 1– 0{ }=

Fig. 5.29 – Modi di deformazione dell’elemento asta a 3 nodi.

184

La figura 5.30 illustra i modi di deformazione relativi ad un elemento traveinflessa, calcolato con integrazione esatta; si notano i due modi relativi ai gradi dilibertà di moto rigido ed i due modi di deformazione, uno di flessione ed uno ditaglio.

La figura 5.31 illustra i modi di deformazione relativi ad un elemento isoparame-trico a quattro nodi in condizioni di stato di deformazione piano, in cui lamatrice di rigidezza è stata calcolata con integrazione esatta utilizzando 2 x 2punti campione; si notano i tre modi relativi ai gradi di libertà di moto rigido edi cinque modi di deformazione, due di flessione, uno di taglio, uno di esten-sione/compressione ed uno di estensione uniforme.


Fig. 5.30 – Modi di deformazione dell’elemento trave inflessa.

Fig. 5.31 – Modi di deformazione dell’elemento a 4 nodi (integrazione esatta).

185

La figura 5.32 illustra invece i modi di deformazione relativi ad un elemento iso-parametrico a quattro nodi in condizioni di stato di deformazione piano, in cuila matrice di rigidezza è stata calcolata con sottointegrazione utilizzando 1 puntocampione; si notano i tre modi relativi ai gradi di libertà di moto rigido, duemodi spuri e tre modi di deformazione, uno di taglio, uno di estensione/com-pressione ed uno di estensione uniforme.


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 isoparame-trico a quattro nodi in condizioni di stato di deformazione piano, in cui lamatrice di rigidezza è stata calcolata con integrazione selettiva utilizzando 2punti campione per le componenti normali della deformazione ed 1 punto cam-pione per la componente di taglio.

Fig. 5.33 – Modi di deformazione dell’elemento a 4 nodi (integrazione selettiva).

186


5.5.5 Patch test

Il patch test, ideato da Irons, è un metodo per verificare se nell'elemento, omeglio 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, inmodo tale che almeno un nodo sia compreso all'interno della struttura così for-mata come ad esempio nella figura 5.34.

Se un elemento supera il patch test è garantita la convergenza al risultato esattoall'aumentare della suddivisione in elementi.

Nella regione così definita si assuma un campo di spostamenti tale da produrreuno stato di deformazione costante.

Il problema così definito viene ora risolto numericamente applicando su tutti inodi esterni spostamenti congruenti con il campo imposto. Se il risultato del cal-colo fornisce, nella regione definita, uno stato di tensione costante ed uguale aquello teorico allora l'elemento in questione ha superato il patch test.

Fig. 5.34 – Patch test.

187

È da notare che il superare o meno il patch test può dipendere fortemente dallageometria dell'elemento utilizzata; ad esempio l'elemento a quattro nodi conmodi incompatibili supera il patch test se ha forma di un parallelogramma, men-tre non lo supera con forme diverse. Di conseguenza un solo patch test non èsufficiente per garantire la convergenza dell'analisi.


188

ESEMPIO 5.5

Nel caso della figura 5.34 sia definito il seguente campo di spostamenti:

5.189

e le coordinate dei nodi siano:

siano inoltre , e .

Gli spostamenti teorici dei nodi interni sono:

e le deformazioni teoriche:

5.190

Le condizioni al contorno per la soluzione numerica sono ricavate dalla leggedegli spostamenti definita ed applicate ai nodi 1, 2, 3, 4:

NODO 1 2 3 4 5 6 7 8

0 40 40 0 10 30 25 12

0 0 40 40 10 15 30 28

NODO 5 6 7 8

1.50 3.75 4.00 2.60

1.50 3.00 4.25 3.40

NODO 1 2 3 4

0.00 4.00 6.00 2.00

0.00 2.00 6.00 4.00

u110------ x y

2---+Ë ¯

Ê �=

v110------ x

2--- y+Ë ¯

Ê �=

X

Y

E 100.000= n 0.3= h 1=

U

V

exx eyy g xy 0.1= = =

U

V

6. SOLIDI ASSIALSIMMETRICI

6.1 I

NTRODUZIONE

Un solido di rivoluzione è un corpo tridimensionale generato per rotazione diuna sezione piana attorno ad un asse (fig.

6.1

). Un solido di rivoluzione è assial-simmetrico se le proprietà geometriche e del materiale sono indipendenti dallaanomalia

.

Se il carico è anch’esso assialsimmetrico il problema si riduce alcaso bidimensionale; ogni punto della sezione è caratterizzato dagli spostamenti

u

(radiale) e

w

(assiale); l’analisi è essenzialmente coincidente con quella del pro-blema piano con una variazione, dovuta alla presenza delle componenti circonfe-renziali della deformazione e della tensione.

J

Fig. 6.1 – Solido assialsimmetrico.

189

SOLIDI

ASSIALSIMMETRICI

Nel caso di carico non assialsimmetrico, il problema può ancora essere ricon-dotto ad un caso bidimensionale, scomponendo il carico in serie di Fourier. Datal’ortogonalità della serie di Fourier il problema è risolto ricorrendo al principiodi sovrapposizione degli effetti, sommando le

n

soluzioni ottenute applicando lecomponenti simmetriche e antisimmetriche del carico, corrispondenti agli

n

ter-mini dell’espressione in serie di Fourier del carico.

Un tipico elemento finito che descrive un solido di rivoluzione è caratterizzatoda

circonferenze nodali

, piuttosto che da

punti nodali

.

6.2 C

ARICO

ASSIALSIMMETRICO

Fig. 6.2 – Elemento triangolare assialsimmetrico.

190

La figura 6.2 illustra un elemento triangolare assialsimmetrico. Nel caso di mate-riale 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 com-ponente circonferenziale, v, nulla per simmetria:

6.1

Il vettore spostamento {u} è espresso in funzione degli spostamenti nodali come:

6.2

u{ }T u w{ }=

uvÓ þ

Ì ýÏ ¸ n1 n2 n3 0 0 0

0 0 0 n1 n2 n3

s{ }=

SOLIDI

ASSIALSIMMETRICI

dove le funzioni di forma dell’elemento triangolare, funzione delle coordinatenodali dell’elemento sono date dalla

4.111

:

6.3

Le deformazioni dell'elemento sono definite dalla

3.94

:

6.4

con:

n11

2A------ x2 y3 x3 y2–( ) y2 y3–( )x x3 x2–( )y+ +[ ]=

n21

2A------ x3 y1 x1 y3–( ) y3 y1–( )x x1 x3–( )y+ +[ ]=

n31

2A------ x1y2 x2 y1–( ) y1 y2–( )x x2 x1–( )y+ +[ ]=

e{ }T err ezz eJJ g rz{ }=

6.5

err r∂∂u=

ezz z∂∂w=

eJJ2p r u+( ) 2pr–

2pr-------------------------------------- u

r--= =

g rz z∂∂u

r∂∂w+=

191


6.6

e derivando:

6.7

Il vettore delle tensioni è:

6.8

e la relazione tensioni deformazioni:

e{ }

r∂∂

0

0 ∂∂z-----

1r--- 0

z∂∂

∂∂r-----

u{ } ∂[ ] n[ ] s{ }= =

e{ } b[ ] s{ }=

s{ }T srr szz sJJ trz{ }=

SOLIDI ASSIALSIMMETRICI

6.9

La matrice di rigidezza è data da:

6.10

l'elemento di volume dV è espresso da: dV = r dr dJ dz = r dr dA. Si ha:

6.11

dove A è l’area della sezione retta dell’elemento assialsimmetrico. L’integrale è disolito valutato numericamente in quanto i termini di [b] sono funzioni logarit-miche.

Nel caso di elementi isoparametrici la matrice di rigidezza si ricava in modo ana-logo a quanto visto in 5.2. Ad esempio per elementi quadrangolari assialsimme-trici (fig. 6.3).

srr

szz

sJJ

trzÓ þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

E1 n+( ) 1 2n–( )

--------------------------------------

1 n– n n 0n 1 n– n 0n n 1 n– 0

0 0 01 2n–

2---------------

err

ezz

eJJ

g rzÓ þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=

k[ ] b[ ]T E[ ] b[ ] VdVÚ=

k[ ] b[ ]T E[ ] b[ ]r Jd Ad0

2p

ÚAÚ 2p b[ ]T E[ ] b[ ]r AdAÚ= =

Fig. 6.3 – Elemento assialsimmetrico a 4 nodi.

192

gli spostamenti u e w e le funzioni di forma ni sono date dalle eq. 5.38 e 5.22:


193

6.12

6.13

Le deformazioni 6.5, sostituendo agli spostamenti u e v le espressioni delle fun-zioni di forma, sono espresse come:

6.14

e in forma matriciale, esplicitando il solo termine i-esimo:

6.15

e la matrice di rigidezza [k]:

6.16

u ni x h,( )uii 1=

4

Â=

w ni x h,( )wii 1=

4

Â=

n114--- 1 x–( ) 1 h–( )= n2

14--- 1 x+( ) 1 h–( )=

n314--- 1 x+( ) 1 h+( )= n4

14--- 1 x–( ) 1 h+( )=

err r∂∂ni ui

i 1=

n

Â=

ezz z∂∂ni wi

i 1=

n

Â=

eJJni

r---- ui

i 1=

n

Â=

g rz z∂∂niui r∂

∂ni wi+Ë ¯Ê �

i 1=

n

Â=

e{ }

º º r∂

∂ni 0 º º

º º 0 z∂

∂ni º º

º º ni

r---- 0 º º

º º z∂

∂ni r∂

∂ni º º

ººui

wi

ººÓ þ

Ô ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

k[ ] b[ ]T E[ ] b[ ]r Jd det J[ ] xd hdp–

p

Ú1–

1

Ú1–

1

Ú=

2p b[ ]T E[ ] b[ ]r det J[ ] xd hd1–

1

Ú1–

1

Ú=


194

Adottando lo schema di integrazione numerica di Gauss:

6.17

dove m è l'ordine di integrazione, wi, wj sono i pesi di integrazione, [b]ij, rij edet [ J ]ij sono calcolati nei punti di integrazione xi, hj. Alcuni termini sotto inte-grale sono funzione diretta di 1/r. Nel caso di integrazione numerica con ilmetodo di Gauss questi termini rimangono finiti poiché i punti di integrazionenon comprendono il caso r = 0. Tuttavia si può avere una forma indeterminatanel calcolo delle deformazioni circonferenziali e delle tensioni perpunti 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 , è quella disostituire, per il calcolo dei valori di tensione per punti giacenti sull’asse z, la rigacorrispondente a 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 scompo-nendo i carichi e gli spostamenti nodali in serie di Fourier lungo la circonferenza.

Nel caso di carico non assialsimmetrico la componente tangenziale dello sposta-mento, v non può più essere considerata nulla, per cui:

6.18

Le componenti radiale circonferenziale e assiale del carico, F, sono:

6.19

Queste possono essere espresse in serie di Fourier come:

6.20

k[ ] 2p wiwj b[ ]ijT E[ ] b[ ]ij rij det J[ ]ij

j 1=

m

Âi 1=

m

Â=

eJJ u r§=

err eJJ=

eJJ

u{ }T u w v{ }=

f{ }T fr fz fc{ }=

fr frs )( h hJ( )cos

h 0=

m

Â frA )( h hJ( )sin

h 0=

m

Â+=

fz fzs )( h hJ( )cos

h 0=

m

Â f zA )( h hJ( )sin

h 0=

m

Â+=

fc fcs )( h hJ( )sin

h 0=

m

Â f cA )( h hJ( )cos

h 0=

m

Â+=


Fig. 6.4 – Componenti simmetriche e antisimmetriche.

dove ciascun termine della serie è anche chiamata armonica ed h è l’ordine di cia-scuna armonica. L’esponente s ed A denotano le componenti simmetriche edantisimmetriche del carico (fig. 6.4), dove simmetria e antisimmetria sono rife-rite rispetto ad un piano di riferimento . La sommatoria è estesa agli mtermini necessari per descrivere correttamente la distribuzione di carico. Inmodo analogo gli spostamenti sono espressi in serie di Fourier:

J 0=

6.21

u r z J, ,( ) uhs hJ( )cos

h 0=

m

Â uhA hJ( )sin

h 0=

m

Â+=

w r z J, ,( ) whs hJ( )cos

h 0=

m

Â whA hJ( )sin

h 0=

m

Â+=

v r z J, ,( ) vhs hJ( )sin

h 0=

m

Â vhA hJ( )cos

h 0=

m

Â+=

195


dove e sono i coefficienti di Fourier per l’armonica diordine h.

In forma matriciale la 6.21 è:

6.22

ovvero:

6.23

con:

6.24

e:

6.25

Per l’armonica h-esima gli spostamenti di un punto dell’elemento sono espressiin funzione degli spostamenti nodali attraverso le funzioni di forma [n]:

6.26

e gli spostamenti complessivi:

uhs uh

A whs wh

A vhs, , , , vh

A

u{ } qhs[ ] uh

s{ }h 0=

m

Â q hA[ ] uh

A{ }h 0=

m

Â+=

u{ } q hs[ ] q h

A[ ]uh

s{ }

uhA{ }Ó þ

Ì ýÏ ¸

Ë ¯Á �Ê �

h 0=

m

Â=

uhs{ }T uh

s whs vh

s{ }= uhA{ }T uh

A whA vh

A{ }=

qhs[ ]

hJ( )cos 0 00 hJ( )cos 00 0 hJ( )sin

=

q hA[ ]

hJ( )sin 0 00 hJ( )sin 00 0 hJ( )cos

=

uhs{ } n[ ] sh

s{ }= uhA{ } n[ ] s h

A{ }=

6.27u{ } qhs[ ] n[ ] sh

s{ }h 0=

m

Â q A[ ] n[ ] s hA{ }

h 0=

m

Â+=

196

ovvero:

6.28u{ } qhs[ ] n[ ] qh

A[ ] n[ ]s h

s{ }

s hA{ }Ó þ

Ì ýÏ ¸

Ë ¯Á �Ê �

h 0=

m

Â=


Le deformazioni dell'elemento sono definite da:

6.29


6.30

con:

6.31

e, sostituendo la 6.27:

6.32

ovvero:

err r∂∂u=

g rJ1r---

J∂∂u v–Ë ¯

Ê �r∂

∂v+=

eJ Jur--=

g zJ z∂∂v 1

r---

J∂∂w

+=

ezz z∂∂w=

g rz z∂∂u

r∂∂w+=

e{ }

∂∂r----- 0 0

1r--- 0

1r--- ∂ ∂J-------

0 ∂∂z----- 0

1r--- ∂ ∂J------- 0 ∂

∂r-----

1r--- –

0 1r--- ∂ ∂J-------

∂∂z-----

∂∂z-----

∂∂r----- 0

u{ } ∂[ ] u{ }= =

e{ }T err eJJ ezz g rJ g zJ g rz{ }=

e{ } ∂[ ] qhs[ ] n[ ] s h

s{ }h 0=

m

Â ∂[ ] qhA[ ] n[ ] s h

A{ }h 0=

m

Â+=

6.33e{ } bhs[ ] s h

s{ }h 0=

m

Â bhA[ ] s h

A{ }h 0=

m

Â+=

197

6.34

Il vettore delle tensioni è:

e{ } bhs[ ] bh

A[ ]sh

s{ }

s hA{ }Ó þ

Ì ýÏ ¸

Ë ¯Á �Ê �

h 0=

m

Â bh[ ] sh{ }h 0=

m

Â= =


6.35

e la relazione tensioni deformazioni:

6.36

e, sostituendo per {e} la 6.33:

s{ }T srr sJJ szz trJ tzJ trz{ }=

srr

sJ J

szz

trJ

tzJ

trzÓ þÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

E1 n+( ) 1 2n–( )

--------------------------------------

1 n– n n 0 0 0n 1 n– n 0 0 0n n 1 n– 0 0 0

0 0 01 2n–

2--------------- 0 0

0 0 0 01 2n–

2--------------- 0

0 0 0 0 01 2n–

2---------------

err

eJ J

ezz

erJ

ezJ

g rzÓ þÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

6.37s{ } E[ ] bhs[ ] sh

s{ }h 0=

m

Â bhA[ ] sh

A{ }h 0=

m

Â+Ë ¯Á �Ê �

=

198

ovvero:

6.38

La matrice di rigidezza è data da:

6.39

avendo indicato con [khl ] il contributo alla rigidezza delle armoniche diordine h e l:

6.40

6.41

l'elemento di volume dV è espresso da: dV = r dr dJ dz = r dr dA. Si ha:

s{ } E[ ] bhs[ ] bh

A[ ]sn

s{ }

snA{ }Ó þ

Ì ýÏ ¸

Ë ¯Á �Ê �

h 0=

m

Â E[ ] bh[ ] sh{ }h 0=

m

Â= =

k[ ] khl[ ]l 0=

m

Âh 0=

m

Â=

khl[ ] bh[ ]T E[ ] bl[ ] VdVÚ=

khl[ ]bh

s[ ]T

bhA[ ]T

E[ ] bls[ ] b l

A[ ] VdVÚ=


6.42

dove A è l’area della sezione retta dell’elemento assialsimmetrico. La 6.42 si puòelaborare come:

6.43

La matrice [k] è quindi:

6.44

Sono quindi da calcolare i seguenti integrali:

khl[ ]bh

s[ ]T

bhA[ ]T

E [ ] bls[ ] bl

A[ ] r Jd Ad0

2p

ÚAÚ=

khl[ ]bh

s[ ]T E[ ] bls[ ] bh

s[ ]T E[ ] blA[ ]

bhA[ ]T E[ ] bl

s[ ] bhA[ ]T E[ ] bl

A[ ]r Jd Ad

0

2p

ÚAÚ=

k[ ]bh

s[ ]T E[ ] bls[ ] bh

s[ ]T E[ ] blA[ ]

bhA[ ]T E[ ] bl

s[ ] bhA[ ]T E[ ] bl

A[ ]r Jd Ad

l 0=

m

Âh 0=

m

Â0

2p

ÚAÚ=

6.45

bhs[ ]T E[ ] bl

s[ ] r Jd Adl 0=

m

Âh 0=

m

Â0

2p

ÚAÚ

bhA[ ]T E[ ] bl

A[ ] r Jd Adl 0=

m

Âh 0=

m

Â0

2p

ÚAÚ

199

6.46

Tutti i termini delle matrici in 6.46 contengono i seguenti integrali, il cui risul-tato è sempre nullo quando integrato tra 0 e 2p:

6.47

I termini simmetrici e antisimmetrici dell’espressione in serie sono quindi disac-coppiati, 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 adue tipi:

bhs[ ]T E[ ] bl

A[ ] r Jd Adl 0=

m

Âh 0=

m

Â0

2p

ÚAÚ

bhA[ ]T E[ ] bl

s[ ] r Jd Adl 0=

m

Âh 0=

m

Â0

2p

ÚAÚ

hJ( )sin hJ( ) Jdcos0

2p

Ú 0= h l, 0 1 2 º m, , , ,=


200

6.48

La matrice [k] sarà quindi composta solo da sottomatrici sulla diagonale princi-pale. Le singole armoniche possono quindi essere analizzate separatamente ed ilrisultato finale può essere ottenuto sovrapponendo i risultati relativi alle m armo-niche, come indicato in 6.27, 6.33, 6.37.

6.49

hJ( )cos hJ( ) Jdcos0

2p

Ú0

p2pÓ

ÔÌÔÏ

=

h lπl n 0π=

l n 0= =

hJ( )sin hJ( ) Jd( )sin0

2p

Ú0

p2pÓ

ÔÌÔÏ

=

h lπl n 0π=

l n 0= =

k[ ]

k00ss[ ] 0 0 0 0 0 0 0 0 0 0 0

0 k11ss[ ] 0 0 0 0 0 0 0 0 0 0

0 0 º 0 0 0 0 0 0 0 0 0

0 0 0 kiiss[ ] 0 0 0 0 0 0 0 0

0 0 0 0 º 0 0 0 0 0 0 0

0 0 0 0 0 kmmss[ ] 0 0 0 0 0 0

0 0 0 0 0 0 k00AA[ ] 0 0 0 0 0

0 0 0 0 0 0 0 k11AA[ ] 0 0 0 0

0 0 0 0 0 0 0 0 º 0 0 0

0 0 0 0 0 0 0 0 0 kiiAA[ ] 0 0

0 0 0 0 0 0 0 0 0 0 º 0

0 0 0 0 0 0 0 0 0 0 0 kmmAA[ ]

=

7. PIASTRE INFLESSE

7.1 R

ICHIAMI

DELLA

TEORIA

DELLE

PIASTRE

INFLESSE

Questo capitolo si riferisce solamente alle piastre piane in campo lineare caricateperpendicolarmente alla loro superficie media; se si suppone (fig.

7.1

) che ilpiano medio coincide con il piano di riferimento

xy

si avrà:

7.1

Verranno inoltre considerate solo piastre sottili (una delle due dimensioni è tra-scurabile rispetto alle altre due) e quindi è possibile fare le seguenti ipotesi:

1. La componente di tensione

s

zz

normale al piano medio può essere tra-scurata rispetto alle altre componenti

s

xx

,

s

yy

,

t

xy

(stato di tensionepiano)

2. Lo spostamento trasversale

w

non varia lungo lo spessore. (questa ipo-tesi è in contrasto con la prima)

3. Si ipotizza altresì, così come fatto per la teoria delle travi inflesse, chesezioni piane perpendicolari al piano medio restino piane dopo defor-mazione. Nell’ipotesi di Kirchoff si ammette inoltre che il contributodel taglio trasversale

t

yz ,

t

zx

sia nullo ed allora segmenti inizialmenteperpendicolari alla superficie media restano tali dopo la deformazione.Nell’ipotesi di Mindlin si considera invece il contributo del taglio tra-sversale e quindi segmenti inizialmente perpendicolari alla superficiemedia in generale non restano tali dopo deformazione

7.1.1 Sistemi di riferimento

L’elemento piastra è caratterizzato, nel suo piano medio, da uno spostamentoverticale

w

, e da due rotazioni

a

x

,

a

y

. È possibile definire spostamenti e rotazionisecondo due sistemi di riferimento, il primo di tipo convenzionale ed il secondopiù utile per la caratterizzazione dell’elemento piastra (fig.

7.1a

e fig.

7.1b

). Si ha:

sxx syy txy 0= = = per z 0=

7.2ax a y–= ay ax=

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à dilunghezza al posto delle tensioni, utilizzare cioè le caratteristiche delle sollecita-zioni (in tal modo si effettua una preintegrazione lungo lo spessore della piastra).

Fig. 7.2 – Tensioni e carichi per unità di lunghezza.

202

Detto h lo spessore della piastra (fig. 7.2) si ha:

7.3

Mx sx z zdh 2§–

h 2§

Ú= My sy z zdh 2§–

h 2§

Ú= Mxy txy z zdh 2§–

h 2§

Ú=

Fx txz zdh 2§–

h 2§

Ú= Fx tyz zdh 2§–

h 2§

Ú=

PIASTRE

INFLESSE

Le tensioni espresse in funzione dei carichi per unità di lunghezza sono:

7.4

ed i valori massimi, in corrispondenza delle superficie di intradosso ed estra-dosso, valgono in valore assoluto:

7.5

Normalmente gli sforzi di taglio

t

xz

e

t

yz

sono trascurabili; essi raggiungono ilvalor massimo per

z

= 0:

7.6

7.2 I

POTESI

DI

K

IRCHOFF

7.2.1 Spostamenti e deformazioni

Nell’ipotesi di Kirchoff le sezioni inizialmente perpendicolari alla superficiemedia restano tali dopo deformazione (fig.

7.3

). Utilizzando il sistema di riferi-mento illustrato in figura

7.1b

, si ha:

7.7

7.8

sx

12Mx

h3------------- z= sy

12My

h3------------- z= txy

12Mxy

h3--------------- z=

sx

6Mx

h2----------= sy

6My

h2----------= txy

6Mxy

h2------------=

txz 1.5 Fx

h-----= tyz 1.5

Fy

h-----=

ax x∂∂w= ay y∂

∂w=

u zx∂

∂w– zax–= = v z

y∂∂w

– zay–= =

Fig. 7.3 – Spostamenti e rotazioni: ipotesi di Kirchoff.

203

PIASTRE INFLESSE

Le deformazioni sono espresse da:

7.9

avendo definito le curvature kx, ky, kxy come:

exx x∂∂u z

x∂∂ax– z

x 2

2

∂∂ w

– z kx–= = = =

eyy y∂∂v z

y∂∂ay– z

y2

2

∂∂ w

– z ky–= = = =

g xy y∂∂u

x∂∂v+ z

y∂∂ax

x∂∂ay+Ë ¯

Ê �– z2xy

2

∂∂ w

– z kxy–= = = =

7.10kx x∂∂ax

x 2

2

∂∂ w= = ky y∂

∂ay

y 2

2

∂∂ w= = kxy

∂ax

∂y---------

x∂∂ay+Ë ¯

Ê � 2xy

2

∂∂ w

= =

7.2.2 Relazioni tensioni-deformazioni

Per una piastra isotropa si ha, nell’ipotesi di Kirchoff:

7.11

in funzione delle curvature e dei carichi per unità di lunghezza si ha:

7.12

e per il termine della temperatura:

sxx

syy

txyÓ þÔ ÔÌ ýÔ ÔÏ ¸

E1 n2–--------------

1 n 0n 1 0

0 01 n–

2------------

exx

eyy

g xyÓ þÔ ÔÌ ýÔ ÔÏ ¸ ax

*T

ay*T

0Ó þÔ ÔÌ ýÔ ÔÏ ¸

–

Ë ¯Á �Á �Á �Ê �

=

exx

eyy

g xyÓ þÔ ÔÌ ýÔ ÔÏ ¸

z

cxx

cyy

cxyÓ þÔ ÔÌ ýÔ ÔÏ ¸

–=

sxx

syy


12zh3

--------

Mx

Mxy

MxyÓ þÔ ÔÌ ýÔ ÔÏ ¸

=

Fig. 7.4 – Gradiente di temperatura.

T2zh-----T0=

204

PIASTRE INFLESSE

7.13

Si ha quindi:

7.14

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 mediaprima della deformazione rimane piana ma non più perpendicolare allasuperficie media (fig. 7.5). Utilizzando il sistema di riferimento illustrato infigura 7.1b, si ha:

7.15

ax*T

ay*T


2zT0

h------------–

ax*

ay*


z

kxo

kyo


–= =kxo ax

*2T0

h---------=

kyo ay*2T0

h---------=

Mx

My


Eh3

12 1 n2–( )-------------------------

1 n 0n 1 0

0 01 n–

2------------

kx

ky

kxyÓ þÔ ÔÌ ýÔ ÔÏ ¸ kxo

kyo


–

Ë ¯Á �Á �Á �Ê �

=

ax x∂∂w g zx–= ay y∂

∂w g zy–=

Fig. 7.5 – Spostamenti e rotazioni: ipotesi di Mindlin.

205

PIASTRE INFLESSE

206

7.16

Le deformazioni sono espresse da:

7.17

avendo definito le curvature kx, ky, kxy, kyz, kzx come:

7.18

7.3.2 Relazioni tensioni-deformazioni

Per una piastra isotropa si ha, nell’ipotesi di Mindlin:

7.19

in funzione delle curvature e dei carichi per unità di lunghezza si ha:

7.20

Si ha quindi:

u zax–= v zay–=

exx x∂∂u z

x∂∂ax– zkx–= = =

eyy y∂∂v z

y∂∂ay– zky–= = =

g xy zy∂

∂ax

x∂∂ay+Ë ¯

Ê �– zkxy–= =

g yz y∂∂w ay– kyz–= =

g zx x∂∂w ax– kzx–= =

kx x∂∂ax= ky y∂

∂ay=

kxy y∂∂ax

x∂∂ay+Ë ¯

Ê �= kyz y∂∂w ay–= kzx x∂

∂w ax–=

sxx

syy


E1 n2–--------------

1 n 0n 1 0

0 01 n–

2------------

exx

eyy

g xyÓ þÔ ÔÌ ýÔ ÔÏ ¸ ax

*T

ay*T


–

Ë ¯Á �Á �Á �Ê �

=

exx

eyy

g xyÓ þÔ ÔÌ ýÔ ÔÏ ¸

z

kx

ky

kxyÓ þÔ ÔÌ ýÔ ÔÏ ¸

–=

g yz

g zxÓ þÌ ýÏ ¸ kyz

kzxÓ þÌ ýÏ ¸

=

sxx

syy


12zh3

--------

Mx

My


=

tyz

tzxÓ þÌ ýÏ ¸ 1.5

h-------

Fy

FxÓ þÌ ýÏ ¸

=

PIASTRE INFLESSE

7.21

7.4 FORMULAZIONE DI RIGIDEZZA NELL’IPOTESI DI KIRCHOFF

Nell’ipotesi in cui si trascura il contributo del taglio alla deformazione (ipotesi diKirchoff ), appare chiaro che se viene assegnato il campo degli spostamenti w, ènecessario calcolarne le derivate seconde per ricavare le deformazioni, e che inol-tre si dovrà formulare un elemento in classe C1, ovvero con continuità sulle deri-vate prime dello spostamento, dovendo assicurare la congruenza (anche lungo lospessore dell’elemento).

La cosa non è impossibile, ma diventa immediatamente piuttosto complicata.Sia ammetta a tale scopo di avere un elemento a quattro nodi rappresentato infigura 7.6 nello spazio naturale x, h.

Mx

My

Mxy

Fy

FxÓ þÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÏ ¸

Eh3

12 1 n2–( )--------------------------

1 n 0 0 0 n 1 0 0 0

0 0 1 n–

2------------ 0 0

0 0 06 1 n–( )

h2 c-------------------- 0

0 0 0 06 1 n–( )

h2 c--------------------

kx

ky

kxy

kyz

kzxÓ þÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÏ ¸ kxo

kyo

000Ó þ

Ô ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÏ ¸

–

Ë ¯Á �Á �Á �Á �Á �Á �Á �Ê �

=

Fig. 7.6 – Elemento piastra quadrangolare.

Dovendo interpolare in classe C1 una funzione f (x, h) si definirà innanzituttouna funzione interpolatrice sul lato 1-2:

7.22f12 G1 x( ) f1 G2 x( ) x∂∂f1 G3 x( ) f2 G4 x( ) x∂

∂f2+ + +=

207

PIASTRE INFLESSE

essendo Gi(x) le funzioni di forma cubiche già viste per le travi inflesse, fi i valoridella funzione f nel nodo i ed il valore che la derivata di f secondo lavariabile h assume nel nodo i. Sul lato 1-2 occorre anche definire indipendente-mente la derivata rispetto a h:

fi∂ h∂§

7.23

Analogamente sul lato opposto, 4-3 si ha:

7.24

7.25

Interpolando secondo h:

7.26

h∂∂f12 G1 x( ) h∂

∂f1 G2 x( ) x h∂

2

∂∂ f1 G3 x( ) h∂

∂f2 G4 x( ) x h∂

2

∂∂ f2+ + +=

f43 G1 x( )f4 G2 x( ) x∂∂f4 G3 x( )f3 G4 x( ) x∂

∂f3+ + +=

h∂∂f43 G1 x( ) h∂

∂f4 G2 x( ) x h∂

2

∂∂ f4 G3 x( ) h∂

∂f3 G4 x( ) x h∂

2

∂∂ f3+ + +=

f x h,( ) G1 h( ) f12 G2 h( ) h∂∂f12 G3 h( ) f43 G4 h( ) h∂

∂f43+ + +=

208

Sostituendo le 7.22, 7.23, 7.24, 7.25 nella 7.26 si ottiene la funzione di interpola-zione bicubica nella quale appaiono come coefficienti i valori, calcolati nei nodi,della funzione, di tutte le sue derivate prime e delle derivate seconde miste; intotale 16 coefficienti.

Per il calcolo delle rotazioni si ha:

7.27

I valori delle derivate di x, h rispetto a x e y si possono ricavare numericamenteinvertendo lo jacobiano, così come si è già visto nel caso degli elementi isopara-metrici membranali.

Le deformazioni si ricavano dalla 7.10:

ax x∂∂w

x∂∂w

x∂∂x◊ ∂w

∂h------- ∂h

∂x------◊+= =

ay y∂∂w

x∂∂w

y∂∂x◊ ∂w

∂h------- ∂h

∂y------◊+= =

PIASTRE INFLESSE

7.28

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 (derivata dire-zionale); questo rende difficile costruire funzioni di forma congruenti. Se siusano come gradi di libertà , vi sarà discontinuità di

.

Si escogitano allora elementi con sottocampi aventi differenti funzioni di forma.Una delle soluzioni possibili si ottiene rinunciando a lavorare in classe C1; si pos-sono definire funzioni interpolatrici indipendenti per lo spostamento e per lerotazioni, lavorando così in classe C0. Occorre però che rotazioni e spostamentosiano effettivamente indipendenti; bisogna quindi abbandonare l’ipotesi di Kir-choff e questo può essere fatto tenendo conto dell’effetto del taglio sulla defor-mata della piastra.

Nei paragrafi seguenti si tratterà prima il problema, più semplice, relativo allatrave per estendere successivamente il tutto al caso bidimensionale.

cx x2

2

∂∂ w

x∂∂x

Ë ¯Ê �

2

x∂∂w

x 2

2

∂∂ x

2 x∂h

2

∂∂ w ∂x

∂x------ ∂h

∂x------ ∂w

∂h-------

x 2

2

∂∂ h

h2

2

∂∂ w

x∂∂h

Ë ¯Ê �

2+ + + +=

cy x2

2

∂∂ w

y∂∂x

Ë ¯Ê �

2

x∂∂w

y 2

2

∂∂ x

2 x∂h

2

∂∂ w ∂x

∂y------ ∂h

∂y------ ∂w

∂h-------

y 2

2

∂∂ h

h2

2

∂∂ w

y∂∂h

Ë ¯Ê �

2+ + + +=

cxy x2

2

∂∂ w ∂x

∂x------ ∂x

∂y------

x∂∂w

x∂y

2

∂∂ x

2 x h∂

2

∂∂ w

x∂∂x

y∂∂h

y∂∂x

x∂∂h

+Ë ¯Ê � ∂w

∂h-------

x∂y

2

∂∂ h

h2

2

∂∂ w ∂h∂h

∂x∂y--------------+ + + +=

w∂ n∂§

w w∂ x∂§ w∂ y∂§, ,∂2w ∂x∂y§

7.5 LA TRAVE DI TIMOSHENKO

Nella letteratura anglosassone si chiama trave di Timoshenko la trave inflessanella quale la rotazione della sezione è svincolata da quella della normale all’asseneutro.

Si ha:

7.29

Lo spostamento v e la rotazione az sono interpolati indipendentemente cosìl’elemento può rappresentare la deformazione dovuta al taglio. È tuttavia danotare che questo elemento, se lineare, può simulare un momento costante manon uno variabile linearmente.

exx x∂∂u y

x∂∂az–= =

g xy x∂∂v az–=

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.

210

Lo spostamento v e la rotazione az sono interpolati indipendentemente come:

7.30

con

v n1 n2v1

v2Ó þÌ ýÏ ¸

=

azn1 n2

az1

az2Ó þÌ ýÏ ¸

=

PIASTRE INFLESSE

7.31

Definendo il vettore degli spostamenti come:

7.32

le deformazioni exx e gxy sono espresse come:

7.33

con

7.34

La matrice di rigidezza dell’elemento può essere calcolata come somma di unaparte flessionale [kF ] ed una di taglio [kT ], essendo i due fenomeni disaccop-piati.

n1 1 xl--–= n2

xl--=

s{ }T v1 az1 v2 az2{ }=

exx yx∂

∂az– y 0x∂

∂n1 0x∂

∂n2 s{ }– bF[ ] s{ }= = =

x∂∂n1 1

l---–=

x∂∂n2 1

l---=

7.35k[ ] kF[ ] kT[ ]+=

con:

7.36

kF[ ] bF[ ]T E[ ] bF[ ] VdVÚ=

kT[ ] bT[ ]T E[ ] bT[ ] VdVÚ=

211

Si ha:

7.37kF[ ] E Jz

0 0 0 0

0x∂

∂n1Ë ¯Ê �

2

0x∂

∂n1

x∂∂n2◊

0 0 0 0

0x∂

∂n1

x∂∂n2◊ 0

x∂∂n1

Ë ¯Ê �

2

xdlÚ=

PIASTRE INFLESSE

212

e:

7.38

Integrando analiticamente (o numericamente con due punti di Gauss) si ha:

7.39

7.40

Se si sottostima la componente di rigidezza dovuta al taglio, sottointegrando lamatrice di rigidezza [kT] con un solo punto di Gauss si ottiene:

kT[ ] GAc

--------

x∂∂n1

Ë ¯Ê �

2

n1 x∂∂n1–

x∂∂n1

x∂∂n2◊ n2 x∂

∂n1–

n1 x∂∂n1–

x∂∂n1

Ë ¯Ê �

2

n1 x∂∂n2– n1n2

x∂∂n1

x∂∂n2◊ n1 x∂

∂n2– x∂

∂n2Ë ¯Ê �

2

n2 x∂∂n2–

n2 x∂∂n1–

x∂∂n1

x∂∂n2◊ n2 x∂

∂n2– x∂

∂n2Ë ¯Ê �

2

xdlÚ=

kF[ ]E Jz

l--------

0 0 0 0 0 1 0 1 –

0 0 0 0 0 1– 0 1

=

kT[ ] GAcl

--------

1 l2--- 1–

l2---

l2---

l 2

3---- l

2---–

l 2

6----

1– l2---– 1

l2--- –

l2---

l 2

6---- l

2---–

l 2

3----

=

PIASTRE INFLESSE

7.41

I vettori dei carichi nodali equivalenti ad un carico distribuito qv e quello dovutoad una differenza di temperatura 2DT fra intradosso ed estradosso sono rispetti-vamente:

7.42

7.5.2 Elemento parabolico

La figura 7.9 rappresenta l’elemento trave di Timoshenko parabolico. Le leggidello spostamento v e della rotazione az sono in tal caso date da funzioni disecondo grado, così come nel caso dell’elemento asta parabolico visto nel cap. 5,eq. 5.65 e 5.66:

7.43

con:

kT[ ] GAcl

--------

1 l2--- 1–

l2---

l2---

l 2

4---- l

2---–

l 2

4----

1– l2---– 1

l2--- –

l2---

l 2

4---- l

2---–

l 2

4----

=

fe{ }qT qv

l2--- 0

l2--- 0

Ó þÌ ýÏ ¸

=

fe{ }DTT

2E Jz a*DT

h--------------------------- 0 1 0 1–{ }=

v n1 n2 n3

v1

v2

v3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=

7.44

n112---x 1 x–( )–= n2

12---x 1 x–( )= n3 1 x2–=

x∂∂n1 1

2---– x+=

x∂∂n2 1

2--- x+=

x∂∂n3 2x–=

213

PIASTRE INFLESSE


214

Supponendo il nodo centrale in mezzeria (x3 = (x1 + x2)/2), lo Jacobiano ed ilsuo determinante valgono, eq. 5.68:

7.45

Ordinando il vettore degli spostamenti come:

7.46

le deformazioni exx e gxy sono date da:

7.47

7.48

J[ ] l2---= J[ ] 1– 2

l---= det J[ ] l

2---=

s{ }T v1 az1 v2 az2 v3 az3{ }=

exx yx∂

∂a– y 0

x∂∂n1 0

x∂∂n2 0

x∂∂n3 s{ }–= =

y2l--- 0 x∂

∂n1 0 x∂∂n2 0 x∂

∂n3s{ }–=

y2l--- 0 1

2---– x+Ë ¯

Ê � 0 12--- x+Ë ¯

Ê � 0 2x–( ) s{ }–=

bF[ ] s{ }=

g xy x∂∂v a–

x∂∂n1 n– 1 x∂

∂n2 n– 2 x∂∂n3 n– 3 s{ }= =

2l--- 1

2---– x+Ë ¯

Ê � 12---x 1 x–( ) 2

l--- 1

2--- x+Ë ¯

Ê � 12---x 1 x+( )–

4l---x– 1 x2–( )– s{ }=

bT[ ] s{ }=

PIASTRE INFLESSE

La matrice di rigidezza si calcola come somma dei due contributi flessionale e ditaglio così come per l’elemento lineare:

7.49

k[ ] kF[ ] kT[ ]+=

kF[ ] bF[ ]T E[ ] bF[ ] Ad( ) xdVÚ=

kT[ ] bT[ ]T E[ ] bT[ ] Ad( ) xdVÚ=

Integrando analiticamente i due termini si ottiene per la parte flessionale (stessitermini dell’elemento asta parabolica 5.71):

7.50kF[ ]2E Jz

l------------

0 0 0 0 0 0

076--- 0

16--- 0 4

3---–

0 0 0 0 0 0

016--- 0

76--- 0 4

3---–

0 0 0 0 0 0

0 43---– 0 4

3---– 0

83---

=

215

e per il contributo dovuto al taglio:

7.51

mentre una sottointegrazione (con due punti di Gauss) fornisce:

kT[ ] GAl2c---------

143l 2-------- 1

l--- 2

3l 2-------- 1

3l----– 16

3l 2--------–

43l----

1l--- 4

15------ 1

3l---- 1

15------– 4

3l----–

215------

23l 2-------- 1

3l---- 14

3l 2-------- 1

l---– 16

3l 2--------– 4

3l----–

13l----– 1

15------– 1

l---–

415------ 4

3l---- 2

15------

163l 2--------– 4

3l----– 16

3l 2--------–

43l---- 32

3l 2-------- 0

43l---- 2

15------ 4

3l----–

215------ 0

1615------

=

PIASTRE INFLESSE

7.52

I vettori dei carichi nodali equivalenti ad un carico uniformemente distribuito qve ad una distribuzione di temperatura 2DT tra intradosso ed estradosso sonorispettivamente:

7.53

kT[ ] GAl2c---------

143l 2-------- 1

l--- 2

3l 2-------- 1

3l----– 16

3l 2--------–

43l----

1l--- 2

9--- 1

3l---- 1

9---– 4

3l----–

29---

23l 2-------- 1

3l---- 14

3l 2-------- 1

l---– 16

3l 2--------– 4

3l----–

13l----– 1

9---– 1

l---–

29--- 4

3l---- 2

9---

163l 2--------– 4

3l----– 16

3l 2--------–

43l---- 32

3l 2-------- 0

43l---- 2

9--- 4

3l----–

29--- 0

89---

=

fe{ }qT

qv l2

------ 13--- 0

13--- 0

43--- 0

Ó þÌ ýÏ ¸

=

fe{ }DTT

2E Jz a*DT

h--------------------------- 0 1 0 1– 0 0{ }=

ESEMPIO 7.1

Calcolare lo spostamento dell’estremo libero della mensola illustrata in figura7.10, modellando la trave sia con un solo elemento lineare che con un soloelemento parabolico.

Fig. 7.10 – Mensola.

216

PIASTRE INFLESSE

La soluzione esatta fornisce per la freccia e la rotazione dell’estremolibero:

7.54

La formulazione di rigidezza nel caso dell’elemento lineare è:

7.55

e, imponendo le condizioni al contorno si ha:

v* a*

v* Pl 3

3E Jz------------– c Pl

GA--------–= a* Pl 2

2E Jz------------–=

E Jz

l--------

0 0 0 00 1 0 1–

0 0 0 00 1– 0 1

GAcl

--------

1l2--- 1–

l2---

l2--- l 2

3----- l

2---–

l 2

6-----

1– l2---– 1 l

2---–

l2--- l 2

6----- l

2---–

l 2

3-----

+

V1

A1

V2

A2Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸ Fv1

Az1

Fv2

Az2Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=

V1 0 Az1 0 Fv2 P Mz2 0=;–=;=;=

7.56E Jz

l--------

0 00 1

GAcl

--------1

l2--- –

l2---–

l 2

3-----

+ V2

A2Ó þÌ ýÏ ¸ P–

0Ó þÌ ýÏ ¸

=

217

Risolvendo il sistema si ricava:

7.57

Ponendo , con raggio giratore d’inerzia della sezione, si ha:

7.58

e nel caso di una trave snella si ha:

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.

V24Pl cGA

------------–3E Jz

c l 2GA+

12E Jz c l 2GA+---------------------------------------◊=

Az22P cGA---------- 1 4

3E Jzc l 2GA+

12E Jz c l 2GA+----------------------------------------–

Ë ¯Á �Ê �

=

Jz r2A= r

V24Pl cGA

-------------–3E c r l§( )2 G+12E c r l§( )2 G+-----------------------------------------=

Az22P cGA---------- 1 4

3E c r l§( )2 G+12E c r l§( )2 G+------------------------------------------–Ë ¯

Ê �=

r l§ 1«

V24Pl cGA

------------–= Az26Pl cGA

------------–=

PIASTRE INFLESSE

218

D’altronde in un caso del genere, essendo e espressi da una legge line-are, ed essendo il taglio costante su tutta la lunghezza della trave, si ha che:

7.60

ma se lo spostamento è lineare allora e allora anche ;

ma all’estremo sinistro, quindi dovunque.

La componente dovuta al taglio prevale quindi su quella flessionale; se si sot-tostima il contributo alla rigidezza dovuto al taglio, utilizzando la matrice dirigidezza sottointegrata, il sistema risolutivo 7.56 diventa:

7.61

Si ricava:

7.62

risposta che approssima meglio il risultato esatto della trave incastrata, for-nendo anzi per ciò che riguarda la rotazione il risultato esatto; con la sottoin-tegrazione della parte di matrice di rigidezza dovuta al taglio (integrazioneselettiva) si evita, almeno in questo caso, il fenomeno del locking.Utilizzando un solo elemento parabolico e imponendo le condizioni al con-torno ( ) la formulazione di rigidezzacon la matrice dovuta al taglio integrata esattamente da:

7.63

Risolvendo il sistema si ricava:

7.64

v az

g xy x∂∂v az– cost= =

vx∂

∂v cost= az cost=

az 0= az 0=

kT[ ]

E Jz

l--------

0 00 1

GAcl

--------1 l

2---–

l2---–

l 2

4-----

+ V2

A2Ó þÌ ýÏ ¸ P–

0Ó þÌ ýÏ ¸

=

V2Pl 3

4E Jz------------– Pl c

GA----------–= Az2

Pl 2

2E Jz------------–=

V1 0 Az1 0 Fv2 P Mz2 0=;–=;=;=

2E Jz

l------------

0 0 0 0

076--- 0 4

3---–

0 0 0 0

0 43---– 0

83---

GAl2c---------

143l 2-------- 1

l---– 16

3l 2--------– 4

3l----–

1l---–

415------ 4

3l---- 2

15------

163l 2--------–

43l---- 32

3l 2-------- 0

43l----–

215------ 0

1615------

+

V2

Az2

V3

Az3Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸ P–

000Ó þ


=

V2Pl

4GA-----------–

240E 2 Jz2 c2 84l 2GAE Jz

c l 4A2G 2+ +

E Jz 60E Jzc l 2GA+( )

----------------------------------------------------------------------------------------------=

Az2Pl 2

2E Jz------------–=

PIASTRE INFLESSE

219

Ponendo dove è il raggio giratore d’inerzia della sezione, si ha:

7.65

Se si sottostima la componente di rigidezza dovuta al taglio il sistema risolu-tivo 7.2 diventa:

7.66

Si ricava:

7.67

risposta che coincide con il risultato esatto della trave incastrata.Utilizzando i seguenti dati:

, , , ,, , A = 100 mm2, e facendo variare il rapporto

tra la lunghezza e l’altezza della trave tra 1 e 10, si ottengono (mante-nendo costante il valore dell’area ) i risultati illustrati nella figura 7.11 siaper l’elemento lineare che per quello parabolico.

Jz r2A= r

V2Pl

4GA-----------–

240E 2 r l§( )4c2 84l 2GAE r l§( )2c l 4A2G 2+ +E r l§( )2 60E r l§( )2c l 2GA+( )

------------------------------------------------------------------------------------------------------------------=

Az2Pl 2

2E Jz-----------–=

2E Jz

l------------

0 0 0 0

076--- 0 4

3---–

0 0 0 0

0 43---– 0

83---

GAl2c---------

143l 2-------- 1

l---– 16

3l 2--------– 4

3l----–

1l---–

29--- 4

3l---- 2

9---

163l 2--------–

43l---- 32

3l 2-------- 0

43l----–

29--- 0

89---

+

V2

Az2

V3

Az3Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸ P–

000Ó þ


=

V2Pl 3

3E Jz------------– Pl c

GA---------–= Az2

Pl 2

2E Jz------------–=

E 210000 N/mm2

= n 0.3= G 80770 N/mm2

= J 8.3 104 mm¥

4=

P 100= l 1000 mm= l h§l h

A

PIASTRE INFLESSE

Fig. 7.11 – Trave incastrata con carico concentrato.

ESEMPIO 7.2


Fig. 7.12 – Mensola.

220


7.68

v* a*

v* Ml 2

2E Jz------------= a* Ml

E Jz--------=

PIASTRE INFLESSE

Utilizzando un solo elemento lineare integrato esattamente e risolvendo ilsistema si ricava:

7.69

mentre un solo elemento lineare con integrazione selettiva fornisce, al pari diun solo elemento parabolico integrato sia esattamente sia in modo selettivo:

7.70

Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nellafigura 7.13.

V26Ml 2c

12E Jz c GAl 2+---------------------------------------= Az2

12Ml c12E Jz c GAl 2+---------------------------------------=

V2Ml 2

2E Jz-------------= Az2

MlE Jz---------=

Fig. 7.13 – Trave incastrata con momento.

ESEMPIO 7.3


Fig. 7.14 – Mensola con carico distribuito.

221

PIASTRE INFLESSE

222


7.71

Utilizzando un solo elemento lineare integrato esattamente e risolvendo ilsistema si ricava:

7.72

mentre un solo elemento lineare con integrazione selettiva fornisce:

7.73

Utilizzando un solo elemento parabolico integrato esattamente e risolvendo ilsistema si ricava:

7.74

mentre un solo elemento parabolico con integrazione selettiva fornisce:

7.75

Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nellafigura 7.15.

v* a*

v* ql 4

8E Jz------------

ql 2

2GA----------- c+= a* ql 3

6E Jz------------=

V22ql 2c

GA---------------

3E Jz c GAl 2+

12E Jz c GAl 2+---------------------------------------◊=

Az2ql 3

4E Jz------------ q GAl 5

4E Jz 12E Jz c GAl 2+( )--------------------------------------------------------–=

V2ql 4

8E Jz------------

ql 2

2GA----------- c+= Az2

ql 3

4E Jz------------=

V2 ql 2360 EJz c( )2 96E Jz GAl 2c GAl 2( )2+ +

12E Jz GA 60E Jz c GAl 2+( )------------------------------------------------------------------------------------------------= Az2

ql 3

6E Jz------------=

V2ql 4

8E Jz------------

ql 2

2GA-----------c+= Az 2

ql 3

6E Jz------------=

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 alproblema iperstatico di una trave caricata in mezzeria e incastrata agliestremi. Calcolare lo spostamento del punto centrale modellando la trave siacon elementi lineari che con elementi parabolici.

Fig. 7.16 – Trave incastrata con carico in mezzeria.

223

PIASTRE INFLESSE


7.76

Utilizzando due elementi lineari con integrazione esatta o selettiva e risol-vendo il sistema si ricava:

7.77

mentre un solo elemento parabolico (a parità quindi di gradi di libertàrispetto alla soluzione precedente) con integrazione esatta o selettiva forni-sce:

7.78

Si nota che manca totalmente la parte dovuta alla flessione pura. Nel caso diuna trave a sezione rettangolare di altezza h e larghezza b si hanno i seguentirapporti rispetto alla soluzione esatta, rispettivamente per gli elementi linearie parabolici:

7.79

7.80

e si ha quindi il risultato esatto solo nel caso trave a snellezza nulla e di dueelementi lineari.Utilizzando i dati dell’esempio 7.1 si ottengono i risultati illustrati nellafigura 7.17.

v* a*

v* Pl 3

24E Jz---------------– Pl

2GA-----------– c= a* 0=

V2Pl

2GA-----------– c= Az2 0=

V2316------ Pl

GA--------◊– c= A z 2 0=

v*V2------ 1

GE c------- l

h--Ë ¯

Ê � 2+=

v*V2------

83--- 1

GE c------- l

h--Ë ¯

Ê � 2+Ë ¯

Ê �=

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 elementiTPE o TPS.

Fig. 7.18 – Trave incastrata con carico concentrato in mezzeria.

225

7.6 MODELLI DISCRETI (LEGAME DISCRETO TRA FLESSIONE E TAGLIO)

Si è visto negli esempi precedenti che nel caso di travi sottili, dove l’effetto deltaglio può essere trascurato rispetto al fenomeno flessionale, l’elemento a trenodi, con ordine di interpolazione parabolico, fornisce risultati migliori dell’ele-mento lineare, a meno che non si adottino tecniche particolari nel caso di questiultimi (ad esempio l’integrazione selettiva).

Una spiegazione di ciò può essere ottenuta considerando l’energia potenziale ela-stica della trave nella sua formulazione di trave di Timoshenko; separando laparte flessionale da quella a taglio si ha, dalle 7.35 e 7.36 o 7.49:

7.81

e dividendo per 1/2EJz:

PE Jz2

--------x∂

∂azË ¯Ê �

2xd

0

l

Ú GmA2c------------

x∂∂v az–Ë ¯

Ê �2

xd0

l

Ú+=

PIASTRE INFLESSE

7.82

Il termine GA/EJz che moltiplica la quota parte di energia potenziale dovuta altaglio assume valori molto grandi nel caso di una trave snella puòessere interpretato come un fattore di penalty. È quindi importante che ilmodello utilizzato (il tipo di elemento e le funzioni di forma scelte) sia in gradodi simulare piccole deformazioni a taglio per grandi valori di GA/EJz in mododa annullare il contributo di energia elastica dovuto al taglio. Se invece, ilsecondo integrale a secondo membro non è piccolo, allora nasce il problema dellocking, perché il termine GA/EJz funziona da moltiplicatore dell’energia ditaglio. Questo è quello che succede nel caso dell’elemento lineare a due nodiintegrato esattamente che quindi non può e non deve essere usato nel caso ditravi snelle. Queste conclusioni sono applicabili anche alle piastre inflesse.

Un modo di superare questo inconveniente è quello di sottointegrare la energiadovuta al taglio, rendendo così singolare la matrice di rigidezza al taglio. Esisteun ulteriore metodo, detto Kirchoff discreto, per superare questo problema. Que-sto metodo consiste nell’imporre l’uguaglianza della rotazione totale dellasezione e derivata dello spostamento in un numero discreto di punti.

Si consideri la trave illustrata in figura 7.19.

Px∂

∂azË ¯Ê �

2

xd0

l

Ú GAE Jz c------------

x∂∂v az–Ë ¯

Ê �2

xd0

l

Ú+=

l h§ �Þ( )


226

Le funzioni di interpolazione siano le 7.44. Si possono eliminare le due variabiliv3 e a3 relative al nodo centrale e imporre quindi due condizioni di Kirchoff.Scelti due punti si impone in essi che:

7.83

sia uguale a:

x q±=

x∂∂v 2

l---

∂n1

∂x--------

∂n2

∂x--------

∂n3

∂x--------◊ ◊

v1

v2

v3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=

PIASTRE INFLESSE

227

7.84

Scrivendo le due equazioni e separando le variabili v3 e a3, si ottiene:

7.85

ovvero:

7.86

Se la matrice [a] è invertibile:

7.87

e, mediante la 7.87 si potrà definire una matrice [t] che permette di eseguire latrasformazione:

7.88

ovvero:

7.89

dove è il vettore delle variabili ridotto.

L’equazione dei lavori virtuali diventa allora:

7.90

a n1 n2 n3

a1

a2

a3Ó þÔ ÔÌ ýÔ ÔÏ ¸

=

x∂∂n3 l

2---n3–

x∂∂n3 l

2---n3–

v3

a3Ó þÌ ýÏ ¸ x∂

∂n1 l2---n1 x∂

∂n2 l2---n2

x∂∂n1 l

2---n1 x∂

∂n2 l2---n2

v1

a1

v2

a2Ó þÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÏ ¸

=

a[ ]v3

a3Ó þÌ ýÏ ¸

b[ ]

v1

a1

v2


=

v3

a3Ó þÌ ýÏ ¸

a[ ] 1– b[ ]

v1

a1

v2


=

v1 a1 v2 a2 v3 a3{ }T t[ ] v1 a1 v2 a2{ }T=

s{ } t[ ] sr{ }=

sr{ }

dsr{ } t[ ]T kF[ ] t[ ] sr{ } dsr{ } t[ ]T f{ }=

PIASTRE INFLESSE

228

Calcolando la 7.85 in :

7.91

che produce:

7.92

la matrice [t] diventa:

7.93

e la matrice [kF] è la 7.50. La matrice di rigidezza [k] coniugata al vettore ridottodelle variabili nodali è:

7.94

Il coefficiente k11 della matrice [k] è:

7.95

mentre quello della matrice di rigidezza dell’elemento hermitiano (a due nodi):

7.96

Questi due coefficienti sono uguali se:

7.97

che comporta:

x q±=

2q–l2--- 1 q2–( )–

2q l2--- 1 q2–( )–

v3

a3Ó þÌ ýÏ ¸

12--- q–

l4--- q2 q–( ) q– 1

2---–

l4--- q2 q+( )

12--- q+

l4--- q2 q+( ) q 1

2---–

l4--- q2 q–( )

v1

a1

v2


=

v3

a3Ó þÌ ýÏ ¸

12--- l

8--- 1

2--- l

8---–

1l 1 q2–( )---------------------– q2

2 1 q2–( )----------------------–

1l 1 q2–( )--------------------- q2

2 1 q2–( )----------------------–

v1

a1

v2


=

t[ ]

1 0 0 00 1 0 012--- l

8--- 1

2--- l

8---–

1l 1 q2–( )---------------------– q2

2 1 q2–( )----------------------–

1l 1 q2–( )--------------------- q2

2 1 q2–( )----------------------–

0 0 1 00 0 0 1

=

k[ ] t[ ]T kF[ ] t[ ]=

k1116EJ

2l3 1 q2–( )2-----------------------------=

k1112EJ

l 3------------=

36 1 q2–( )2 16=

PIASTRE INFLESSE

7.98

È un fatto notevole che questi siano proprio i valori delle ascisse dei punti diintegrazione di Gauss. Si può verificare che in tali ascisse anche tutti gli altricoefficienti, e non solo l’elemento k11, coincidono con quelli della matrice her-mitiana.

Se ne deduce che, accettando di trascurare gli effetti del taglio, i risultati migliorisi ottengono imponendo le condizioni di Kirchoff proprio nei due punti di inte-grazione 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 Mindlinconsidera la rotazione del piano medio separata dalla rotazione della sezione. Sihanno le relazioni viste in 7.3.1:

7.99

Lo spostamento w e le rotazioni ax e ay sono descritti indipendentemente come:

7.100

ovvero:

q 1

3-------±=

exx x∂∂u z

x∂∂ax– zkx–= = =

eyy y∂∂v z

y∂∂ay– zky–= = =

g xy zy∂

∂ax

x∂∂ay+Ë ¯

Ê �– zkxy–= =

g yz y∂∂w ay– k yz–= =

g zx x∂∂w ax– kzx–= =

w ni x h,( )wii 1=

n

Â= ax ni x h,( )axi

i 1=

n

Â= ay ni x h,( )ayi

i 1=

n

Â=

7.101u{ } n[ ] s{ }=

229

La matrice di rigidezza [k] si può ottenere come somma della parte flessionale edi quella dovuta al taglio, poiché i due fenomeni sono disaccoppiati:

7.102

È conveniente riferirsi a momenti per unità di lunghezza e curvature, cioè allecaratteristiche di sollecitazione (momenti flettenti e forze di taglio). Le deforma-zioni sono espresse da:

k[ ] kF[ ] kT[ ]+ bF[ ]T EF[ ] bF[ ] VdVÚ bT[ ]T ET[ ] bT[ ] Vd

VÚ+= =

PIASTRE INFLESSE

230

7.103


7.104

7.105

kx x∂∂ax

x∂∂niaxi

i 1=

n

Â= =

ky y∂∂ay

y∂∂niayi

i 1=

n

Â= =

kxy y∂∂ax

x∂∂ay+

y∂∂niaxi x∂

∂niayi+Ë ¯

Ê �

i 1=

n

Â= =

kyz y∂∂w ay+

y∂∂niwi niayi

+Ë ¯Ê �

i 1=

n

Â= =

kzx x∂∂w ax+

x∂∂niwi niaxi

+Ë ¯Ê �

i 1=

n

Â= =

kx

ky

kxy

kyz

kzxÓ þÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÏ ¸

º º 0 x∂

∂ni 0 º º

º º 0 0 y∂

∂ni º º

º º 0 y∂

∂ni x∂

∂ni º º

º º y∂

∂ni 0 ni– º º

º º x∂

∂ni ni – 0 º º

ººwi

axi

ayi

ººÓ þ

Ô ÔÔ ÔÔ ÔÔ ÔÔ ÔÌ ýÔ ÔÔ ÔÔ ÔÔ ÔÔ ÔÏ ¸

=

k{ }g{ }Ó þ

Ì ýÏ ¸ º º kF i,[ ] º º

º º kT i,[ ] º º

ººsi{ }

ººÓ þ


bF[ ]

bT[ ]s{ }= =

PIASTRE INFLESSE

231

e le derivate si calcolano come in 5.51.

La relazione carichi curvature è:

7.106

e le matrici di rigidezza sono normalmente calcolate mediante integrazionenumerica:

7.107

7.7.2 Carico nodale equivalente a carico superficiale distribuito

Sia qv il carico per unità di superficie agente perpendicolarmente alla superficiemedia dell’elemento piastra. Il vettore del carico nodale equivalente è:

7.108

con:

7.109

e [n] è la matrice delle funzioni di forma dell’elemento 7.101. Il calcolo delcarico nodale equivalente si esegue numericamente come:

7.110

x∂∂ni

y∂∂ni,

EF[ ] Eh3

12 1 n2–( )-------------------------

1 n 0n 1 0

0 01 n2–

2--------------

=

ET[ ] Eh2 1 n+( ) c------------------------ 1 0

0 1=

kF[ ] wiwj bF[ ]ijT EF[ ] bF[ ]

ijhij det J ][ ij

j 1=

m

Âi 1=

m

Â=

kT[ ] wiwj bT[ ]ijT ET[ ] bT[ ]

ijhij det J ][ ij

j 1=

m

Âi 1=

m

Â=

fe{ }q n[ ]T q{ } AdAÚ=

fe{ }qT Fzi

MxiMyi

{ }=

q{ }T q 0 0{ }=

fe{ }q nij[ ] qij{ }wiwj det J ][ ijj 1=

m

Âi 1=

m

Â=

PIASTRE INFLESSE

ESEMPIO 7.5

È data una piastra circolare non forata, incastrata al bordo esterno, di raggioesterno e spessore costante . La piastra realizzata inacciaio, con modulo elastico e coefficiente di Poisson

, è soggetta ad un carico uniformemente distribuito pari a .Calcolare la deformata ed i momenti flettenti radiale e circonferenziale. Lasoluzione classica del problema fornisce i seguenti risultati:

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 esaminatedue suddivisioni, rispettivamente con 5 e 9 elementi e 72 e 120 gradi dilibertà (fig. 7.20)

R 50 mm= h 2mm=E 200000 N/mm

2=

n 0.3=p 1 N/mm

2=

w12 1 n2–( )p

64Eh3---------------------------- R2 r2–( )2 4h2

1 n–( )----------------- R2 r2–( )+=

mrp

16------ R2 1 n+( ) r2 3 n+( )–[ ]=

mcp

16------ R2 1 n+( ) r2 1 3n+( )–[ ]=

Fig. 7.20 – Schematizzazione di piastra circolare: a. 5 elementi, b. 9 elementi.

La tabella 7.1 riassume i risultati ottenuti dal calcolo numerico confrontaticon quelli analitici.

Tab. 7.1 – Valori massimi della freccia e dei momenti flettenti

soluzione

ANALITICO 0.669 203.1

FEM 5 EL. 0.659 205.8

FEM 9 EL. 0.667 203.4

w mr mc( )

232

PIASTRE INFLESSE

La figura 7.21 illustra la deformata calcolata con le due suddivisioni, mentrein 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 sulla precisione del calcolo; la figura 7.23 riporta i risultati ottenuti per variabile da 0.01 a 0.4.

h R§h R§

233

Fig. 7.23 – Andamento della freccia in funzione del rapporto h/R.

Elementi finiti - Parte III

Documents

Transcript of Elementi finiti - Parte III