Appunti di Progettazione agli Elementi Finiti di Strutture Metalliche + Progetto ANSYS [PEFSM - ex PASM]

Part III

Metodi numerici per la risoluzione

di sistemi di equazioni lineari

Abbiamo visto che il metodo agli elementi finiti permette di ricondurre la risoluzione

di problemi statici alla soluzione di un sistema di equazioni lineari del tipo =

Ax b.

Poiché tipicamente si ha a che fare con sistemi aventi un elevato numero di gradi di

libertà, il numero di equazioni componenti il sistema è generalmente molto elevato. Di

conseguenza, l’approccio che prevede l’ottenimento della soluzione tramite inversione

della matrice dei coefficienti non è praticabile per il notevole costo computazionale

che verrebbe richiesto in termini di memoria.

Per tale motivo, vengono usualmente adottate metodologie differenti. In partico-

lare, possono essere utilizzati metodi diretti (che forniscono la soluzione esatta del

sistema) o metodi iterativi (che forniscono una soluzione approssimata con un grado

di approssimazione che viene opportunamente scelto).

In questa sezione, vedremo solo alcuni dei possibili metodi adottabili. Nello speci-

fico, ci occuperemo del metodo diretto basato sulla fattorizzazione di Cholesky della

matrice e dei metodi iterativi di Jacobi e Gauss-Seidel.

A

1 Metodi diretti

1.1 Fattorizzazione di Cholesky

Sia una matrice quadrata, reale, simmetrica e definita positiva, tale matrice può

A

essere decomposta come: T

=

A LL

1 T

con matrice triangolare inferiore con elementi diagonali positivi. è la trasposta

L L

di L.

1.1.1 Algoritmo di Cholesky-Banachiewicz

L’algoritmo di Cholesky-Banachiewicz permette di calcolare direttamente gli elementi

della matrice triangolare inferiore Esso inizia con il calcolo del termine superiore

L.

sinistro della matrice e procede a calcolare la matrice riga per riga:

∀ i = 1, ..., n

∀ j = 1, ..., (i − 1)

j−1 !

1 X

− l l

a

l =

i,j i,j i,k j,k

l j,j k=1

v i−1

u X

u 2

= a − l .

l t i,i

i,i i,k

k=1

dove a ed l sono gli elementi delle matrici ed rispettivamente.

A L,

i,j i,j =

1.2 Soluzione del sistema di equazioni lineari Ax b

Dato il sistema di equazioni lineari = con matrice simmetrica e definita

Ax b, A

positiva, è possibile risolvere il sistema sfruttando la fattorizzazione di Cholesky. Una

volta calcolata la matrice si procede con la risoluzione dei due sistemi lineari:

L, T

= = (1)

Ly b; L x y

T matrici triangolari, la risoluzione dei

Si noti che essendo e la sua trasposta

L L

due sistemi la si può effettuare rapidamente in accordo al seguente algoritmo (che

non prevede l’inversione diretta delle matrici):

!

i−1

1 X

−

y b

= l y per i = 1, ..., n (2)

i i,k k

i l i,i k=1 !

n

1 X

= −

y l x per i = n, ..., 1 (3)

x i i k,i k

l i,i k=i+1 2

1.2.1 Esempio di applicazione dell’algoritmo di Cholesky-Banachiewicz

Risolvere, sfruttando la fattorizzazione di Cholesky, il seguente sistema lineare =

Ax b,

dove:   



  x 3

4 −1 −1  



 1  



  

 

  

  

; =

= ; =

x

A b 9

−1 6 1 x



 2



  

 



  

 

−1 1 7 −6

x

   

 

 

3

Calcolo matrice di fattorizzazione L.

√ a = 2

=

l 1,1

1,1 1 a = −1/2

=

l 2,1 2,1

l 1,1 √

q 2

a − l = 23/2

l =

2,2 2,2 2,1

1

l = a = −1/2

3,1 3,1

l 1,1 √

1

= (a − l l ) = 3/ 2 23

l 3,2 3,2 3,1 2,1

l 2,2

q p

2 2

l = a − l − l = 3 17/23

3,3

3,3 3,1 3,2

Calcolo vettore intermedio y.

b 3

1

=

y =

1 l 2

1,1

1 39

√

= (b − l y ) =

y 2 2 2,1 1

l 2 23

2,2

1 50

√

= (b − l y − l y ) = −

y 3 3 3,1 1 3,2 2

l 391

3,3

Calcolo soluzione x. 3

y 49

3

x = =

3 l 51

3,3 31

1

x (y − l x ) =

= 2 3,2 3

2 l 17

2,2

1 50

x = (y − l x − l x ) = −

1 1 2,1 2 3,1 3

l 51

1,1

2 Metodi iterativi

Ci occuperemo dei seguenti metodi iterativi:

• Metodo di Jacobi;

• Metodo di Gauss-Seidel.

2.1 Metodo di Jacobi

Si supponga di dover risolvere il sistema di equazioni lineari = dove

Kx b,



 k k . . . k

k 11 12 13 1n 

 k k k . . . k 

 21 22 23 2n 

 



=

K k k k . . . k

 

31 32 33 3n 

 .. ..

.. ..

..

 

.

. .

. . 



 



 k k k . . . k

n1 n2 n3 nn

La matrice può essere scritta come la somma di due matrici e

K A B:

= + (4)

K A B

in cui è una matrice triangolare inferiore e è una matrice triangolare strettamente

A B

4

superiore 

 

 0 k

k 0 0 ... 0 k . . . k

11 12 13 1n



 

 k 0 . . . 0

k 0 0 k . . . k





 

21 22 23 2n



 

 

 

 ; =

= B

A 0 0 0 . . . k

k k k . . . 0 

 

 31 32 33 3n





 

.. ..

..

.. .. .. .. ...

... 



 

k

. . .

.

. . . 



 

in





 





 

k k k . . . k 0 0 0 ... 0

n1 n2 n3 nn

L’equazione = è, pertanto, riscrivibile come:

Kx b = − . (5)

Ax b Bx

Se per il vettore incognito si assume una soluzione di primo tentativo:

x 

 x 

 1 

 

 

 

 x 

 

 2 

 

 



(0) ,

=

x x 3

.. 

 

 . 

 

 

 

 



 

x

 

 

n (1)

è possibile calcolare la soluzione all’iterazione successiva :

x

(0)

(1) −1

= (b − ) . (6)

A Bx

x

(1) (2)

Noto il vettore è possibile ricavare :

x x (1)

(2) −1

= (b − ) . (7)

A Bx

x

e, in generale, tale schema può essere implementato iterativamente

(k)

(k+1) −1

= (b − ) . (8)

x A Bx

Il processo di iterazione sarà arrestato quando la differenza tra la soluzione all’iterazione

k + 1 e quella precedente all’iterazione k è sufficientemente piccola, inferiore a un

5

valore stabilito. La facilità di procedere in questo modo risiede nel fatto che si utiliz-

zano matrici trangolari inferiori e superiori, per cui non è in realtà necessario effet-

tuare l’inversione diretta della matrice Infatti, il generico termine della soluzione

A.

all’iterazione k + 1 può scriversi come

n

1 (k)

(k+1) X

(b − k x ) i = 1, 2, 3, ..., n. (9)

x = i ij j

i k ii =

j=1j6 i

Il metodo di Jacobi, consiste dunque nell’applicare iterativamente l’equazione (9),

(0)

partendo da una soluzione arbitraria di primo tentativo . Il processo di iterazione

x

si arresta quando (k+1) (k)

− k < ε (10)

kx x

essendo ε un valore arbitrariamente piccolo. In alternativa, è possibile definire l’errore

ε in termini percentuali e, quindi, arrestare il processo iterativo quando

(k+1) (k)

kx − k

x · 100 < ε (%). (11)

(k)

kx k

2.1.1 Esempio di applicazione del metodo di Jacobi

Si consideri il seguente sistema di equazioni lineari = dove

Kx b, 









 3 1 0 x

1 

 

 1



 

 







 



 



; =

; =

= x

b

K x

1

1 5 2

  2

  



 

   

 

2

0 2 7 x

 

 

 

 

3 (k+1) −

e si supponga di volerlo risolvere con il metodo di Jacobi, da arrestare quando kx

(k) (k)

k/kx k < 0.01.

x Si assuma come soluzione di primo tentativo la seguente:

 

0 

 

 

 



(0) =

x 0 

 

 1 

 



6

e si applichi iterativamente l’equazione (9). Si avrà: 1

1

(1)

x (1 − 1 · 0 − 0 · 1) = (12)

=

1 3 3 1

1

(1) =

x (1 − 1 · 0 − 2 · 1) = − (13)

2 5 5

2

1

(1) =

x (2 − 0 · 1 − 2 · 0) = . (14)

3 7 7

Iterando una seconda volta:

1 2

1

(2) (1 − 1 · − (15)

− 0) =

=

x 1 3 5 5

1 2

1 2

(2) − 2 · (16)

1 − 1 · =

=

x 2 5 3 7 105

1

2 12

1

(2)

x − 2 · − . (17)

2 − 0 · =

=

3 7 7 5 35

e cosı̀ via. Il processò iterativo verrà arrestato quando

v 2

u (k+1) (k)

3

P − x

x

u i i

i=1

u < 0.01 (18)

u 2

(k)

t 3

P x i

i=1

che nel nostro caso avviene dopo 7 iterazioni, producendo la seguente soluzione

(7)

x = 0.325627 (19)

1

(7)

x = 0.021933 (20)

2

(7) = 0.279109 (21)

x 3

con un errore v 2

u (7) (6)

3

P x − x

u i i

i=1

u = 0.007777 < 0.01. (22)

u 2

(6)

t 3

P x i

i=1

La soluzione esatta del sistema di partenza è, invece:

= 0.3255814 (23)

x 1 = 0.0232558 (24)

x 2

x = 0.2790698 (25)

3 7

2.2 Metodo di Gauss-Seidel

Il metodo iterativo di Gauss-Seidel è una variante del metodo di Jacobi, usata per

accelerare la convergenza della soluzione. La formula iterativa del metodo di Gauss-

Seidel è analoga a quella di Jacobi, ma si differenzia perché nel calcolo del termine

(k+1) (k+1)

vengono utilizzati anche i termini di già calcolati. Di conseguenza

x x

i

l’equazione (9) si modifica come segue

i−1 n

1

(k+1) (k+1) (k)

X X

x = (b − k x − k x ) i = 1, 2, 3, ..., n. (26)

i ij ij

i j j

k ii j=1 j=i+1

La velocità di convergenza è aumentata per il fatto che nell’equazione (26) sono

presenti dei valori aggiornati delle incognite. Il metodo di Gauss-Seidel è sempre

convergente quando la matrice dei coefficienti è simmetrica e definita positiva o è

K nj=1,j6

P

| > |k |).

a diagonale dominante in senso stretto (cioè vale la condizione |k ii ij

= i

Tuttavia, non è infrequente che il metodo di Gauss-Seidel converga anche quando tali

condizioni non sono soddisfatte.

Anche in questo caso, il processo iterativo verrà arrestato quando sarà soddisfatta

la condizione (10) (o la (11)).

2.2.1 Esempio di applicazione del metodo di Gauss-Seidel

Si consideri il sistema di equazioni = preso in esame nel paragrafo precedente

Kx b,

 

 



 x

1

3 1 0  

 

1

 

 

 

 



 



 

; =

; =

= .

x

b

K 1

1 5 2 x



 2

  



 



 



 

2

0 2 7 x





 

 

 

3

Partendo dalla stessa soluzione di primo tentativo



 0

 

 



 



(0) =

x 0 

 

 1 

 



8

e applicando iterativamente la formula (26), si ha:

1

1

(1)

x (1 − 1 · 0 − 0 · 1) = = 0.3333 (27)

=

1 3 3

1 4

1

(1) (1 − 1 · − 2 · 1) = − = −0.266667 (28)

x =

2 5 3 15

1 4 38

1

(1) (2 − 0 · + 2 · ) = = 0.361905 . (29)

x =

3 7 3 15 105

Iterando una seconda volta:

1 19

4

(2) =

x − 0) =

(1 − 1 · − = 0.422222 (30)

1 3 15 45

38

19 46

1

(2) − 2 · = −0.0292063 (31)

1 − 1 · = −

=

x 2 5 45 105 1575

1 46

19 3242

(2) =

x − 2 · − = 0.294059 . (32)

2 − 0 · =

3 7 45 1575 11025

e cosı̀ via. Il processò iterativo verrà arrestato quando

v 2

u (k+1) (k)

3

P − x

x

u i i

i=1

u < 0.01 (33)

u 2

(k)

t 3

P x i

i=1

che nel nostro caso avviene dopo 5 iterazioni (quindi prima rispetto al metodo di

Jacobi), producendo la seguente soluzione

(5)

x = 0.326154 (34)

1

(5)

x = 0.022945 (35)

2

(5)

x = 0.279159 (36)

3

con un errore v 2

u (5) (4)

3

P − x

x

u i i

i=1

u = 0.00688838 < 0.01. (37)

u 2

(4)

t 3

P x i

i=1 9

Part V

Elemento asta

Gli elementi asta sono in grado di supportare esclusivamente azioni assiali e, quindi,

esplicano una reazione solo in tale direzione. Ricordiamo che un sistema è schematiz-

zabile con elementi di tipo asta quando gli elementi sono sollecitati solo assialmente

e la sezione geometrica ha dimensioni molto più piccole rispetto alla lunghezza degli

elementi stessi.

L’elemento asta è, pertanto, utilizzabile nell’analisi di sistemi strutturali quali,

ad esempio, le travature reticolari. Ovviamente, affichè una generica travatura possa

schematizzarsi con elementi asta è necessario che i vari elementi della stessa siano

connessi tra loro mediante perni o cerniere, cosı̀ che non vengano trasmessi momenti,

ma solo forze.

Inoltre, è possibile definire un elemento 1D analogo all’elemento asta a sforzo

normale, che però sia in grado di supportare momenti torcenti e, quindi, di deformarsi

solo a torsione. In questo capitolo, ci occuperemo anche dell’elemento soggetto a

pura torsione, per la cui formulazione si potranno adottare le stesse funzioni di forma

impiegate per quello soggetto a solo sforzo assiale. Di entrambi daremo le matrici di

rigidezza e delle masse (potendo poi calcolare la matrice di dissipazione con

K M C

il metodo di Rayleigh).

1 Elemento a sforzo normale

Con riferiemnto alla Figura 1, dove è rappresentato un elemento asta, si assuma l’asse

x del sistema di riferimento locale coincidente con l’asse dell’asta e orientato dal nodo

1 verso il nodo 2.

In tale riferimento locale, gli spostamenti significativi sono quelli lungo la direzione

x, u e u . In accordo a quanto visto in precedenza, lo spostamento u (x) di un gener-

1 2 1

Figure 1: Elemento asta. Gradi di libertà nel riferimento locale e globale.

ico punto interno all’elemento può ottenersi come interpolazione degli spostamenti

nodali tramite apportune funzioni di forma: 

 



[ ] u 1 = N (x)u + N (x)u , (1)

u(x) = =

N(x)d N (x) N (x)

e 1 1 2 2

1 2 

 u 2

dove N (x) e N (x) sono funzioni polinomiali che, nel caso specifico, si assumeranno

1 2

lineari.

Se N e N sono funzioni lineari, lo sarà ovviamente anche lo spostamento u (x):

1 2 u(x) = α + α x . (2)

0 1

Le costanti α e α possono determinarsi imponendo le condizioni al contorno sui

0 1

nodi, e cioè u(x = 0) = u (3)

1

u(x = L) = u . (4)

2

2

da cui, α = u e α = (u − u ) /l , essendo l la lunghezza dell’elemento. Sos-

e e

0 1 1 2 1

tituendo nell’equazione precedente: ( ) ( )

u − u x x

2 1

u(x) = u + x = 1 − u + u . (5)

1 1 2

l l l

e e e

Confrontando le relazioni (1) e (5), si ottiene per la matrice delle funzioni di forma

la seguente espressione ]

[ x

x , . (6)

= 1 −

N(x) l l

e e

E’ facile verificare che le due funzioni di forma trovate soddisfano le proprietà

fondamentali delle funzioni di forma.

1.1 Matrice di rigidezza K e

Dalle equazioni di congruenza e costitutiva, si ricavache la deformazione e la tensione

all’interno dell’elemento sono costanti. Infatti, avremo

∂u u − u

2 1

ε = = , (7)

xx ∂x l e

u − u

2 1 . (8)

σ = E ε = E

xx xx l e ∂

Essendo ε = = = , dove ≡ , la matrice di compatibilità

LU LNd Bd L B

xx e e ∂x

può scriversi come [ ]

] [

∂N ∂N 1

1

∂N 1 2

= , , , (9)

= −

=

B ∂x ∂x ∂x l l e

2

La matrice di rigidezza per definizione è:

k e ∫ T

= dV ,

k B c B

e V

dove in questo caso ≡ E (essendo E il modulo di Young). Assumendo costante la

c

sezione dell’elemento, si ha e, pertanto,

V=A dx

 



 ∫

∫ [ ] 1

1

1 l

l −

− e

e  



 2 2

l l l dx,

= E A dx = E A

1 1

k −

e e e

e l le 1

1 1

−

e 0

0 2 2

l l l

e e e

3

ovvero: 

 1 −1

E A 



= . (10)

k e l −1 1

e

1.1.1 Matrice del cambiamento di riferimento T

Un elemento asta, in generale, sarà parte di un sistema con più aste aventi differenti

orientazioni. Poiché le matrici di rigi