Relazione parte 3 Metodi Numerici

Alma Mater Studiorum – Università di Bologna

Scuola di ingegneria e architettura

Sede di Bologna

Corso di laurea magistrale in ingegneria civile

Relazione esercitazioni

Metodi numerici per l'ingegneria civile M

Professore: Custodi Alberto
Gruppo progetto:

• Donini Giovanni - Matricola: 896077
• Ferrari Mattia - Matricola: 952362
• Diana Francesco - Matricola: 9195061

Esercizio 01 - Calcolo pi greco

Per la determinazione di pi greco si è proceduto al calcolo dei perimetri dei poligoni inscritti ad una circonferenza di raggio unitario r=1, aumentando progressivamente il numero di lati attraverso un ciclo for in Matlab.

Una volta calcolato il perimetro di un poligono oggetto dello studio, si è imposta l'eguaglianza tra lo stesso e il perimetro della generica circonferenza, pari a 2rπ. Nel caso studio, essendo la circonferenza di raggio unitario, si ha 2π.

Per il calcolo dei perimetri dei poligoni inscritti si è utilizzata la seguente formula, ponendola successivamente uguale a 2rπ, estraendo poi π approssimato per difetto: π = Lnew · n dove Lnew è la lunghezza dell’n-esimo lato del nuovo poligono (inscritto nella circonferenza) oggetto di studio ed è calcolato con la formula: 2√(1 - (Lold/2)²). Lold è la lunghezza del lato del poligono inscritto relativo alla iterazione precedente.

Una volta calcolato il π approssimato, si aumenteranno progressivamente i lati dei poligoni fino a cercare di arrivare a convergenza tra i due π approssimati. Come già detto, il programma lavora variando il numero dei lati per poligoni inscritti con n.

Un incremento del numero di lati ad ogni iterazione pari a 2, con n volta mandato in esecuzione, si ottiene un'approssimazione aumentando il numero di lati attraverso un ciclo for, pi greco come all'interno di differenza fra il pi greco di macchina (che è memorizzato Matlab) e il valore approssimato calcolato. Si è optato per l'utilizzo del comando vpa. Questo perché il comando vpa fondamentalmente opera un aumento di cifre significative sul calcolo aritmetico. Il comando è stato utilizzato perché nel caso oggetto di studio arrecava una migliore approssimazione legata di cifre decimali.

Aumentando il numero di lati ci si aspetterebbe un progressivo aumento della precisione, cosa che però non accade. Questo dipende dal numero di cifre decimali con le quali stiamo lavorando. Esso è un problema legato alla precisione di macchina. Aumentando il numero di iterazioni si può notare che l’errore, che per un i=20 iniziale, mostrava un andamento solo decrescente (e quindi un miglioramento dell’approssimazione), per i=50 ha una diminuzione fino ad un valore di circa -6,7 · 10^-17, per poi aumentare fino alla fine del ciclo.

Questo significa che l’aumento del numero dei lati del poligono non porta necessariamente ad un miglioramento della precisione dei risultati. Il motivo è da attribuirsi all’approssimazione delle cifre decimali dopo la virgola, che dopo un certo valore vengono arrotondate secondo le direttive di macchina.

Vengono ora riportati i grafici dell’errore:

• Figura 01.1. Andamento dell’errore partendo da un quadrato.
• Figura 01.2. Andamento dell’errore partendo da un triangolo.

Si confrontano ora i risultati ottenuti su Matlab con quelli ottenuti su Excel. I risultati ottenuti su Excel sono inizialmente analoghi a quelli ottenuti su Matlab; dalla decima iterazione però, i risultati incominciano a discostarsi da quelli ottenuti su Matlab e ad essere piuttosto imprecisi.

Arrivati alla 13-esima iterazione, si comincia a notare un frequente cambio di segno nell’errore, indice che ormai il grado di precisione non è più affidabile. Infine, superate circa le 20 iterazioni, il programma non è più in grado di calcolare la lunghezza del lato del nuovo poligono e di conseguenza vengono stampati una serie di zeri.

Nel caso di Excel, anche se Excel è in grado di archiviare numeri da 1.79769313486232 E308 a 2.2250738585072 E-308, è possibile farlo solo entro 15 cifre di precisione. Questa limitazione è il risultato diretto dell'applicazione rigorosa della specifica IEEE 754 e non è una limitazione di Excel.

Quindi si può dedurre che l’implementazione di un programma di questo genere, che richiede una certa precisione nelle cifre significative dopo la virgola, non può essere affidata ad un programma come Excel. È quindi da preferire in tal caso un’analisi numerica su Matlab.

Se si utilizza come poligono iniziale un quadrato, l’errore ha un andamento che inizialmente decresce con il numero di iterazioni e dopo un certo valore inizia a crescere. Cambiando il poligono iniziale da quadrato a triangolo si può notare che l’andamento dell’errore è simile a quello del quadrato.

Si può allora dedurre che la formula utilizzata nel calcolatore per l’approssimazione del pi greco vada bene sia nel caso del quadrato che del triangolo inscritto.

Esercitazione 1 - Strutture reticolari

Struttura reticolare A

Si analizza la struttura in Figura 1.1, assumendo i seguenti dati:

• h = 3 m
• L = 4 m
• E = 30 x 10⁶ kN/m²
• Aste di parete: A = 0.02 m²
• Aste corrente superiore e inferiore: A = 0.045 m
• F = 1500 kN

Figura 1.1. Struttura reticolare A: dimensioni e carichi assegnati.

Di seguito si riportano gli spostamenti, le reazioni vincolari, le sollecitazioni nelle aste (indicando con gli sforzi positivi quelli di trazione) e la matrice di rigidezza per la reticolare A calcolati utilizzando il programma “Ret2d_incl”:

Spostamenti nodali

Reazioni vincolari

Sforzi nelle aste

Figura 1.2. Deformata della reticolare A.

Figura 1.3. Diagramma degli sforzi normali della reticolare A (verde = positivo = trazione; rosso = negativo = compressione).

Struttura reticolare B

La seconda struttura reticolare presenta la medesima geometria e le medesime condizioni di carico della precedente, tuttavia viene modificata la numerazione dei nodi costituenti, come indicato in Figura 1.4.

Figura 1.4. Struttura reticolare B.

Spostamenti nodali

Reazioni vincolari

Sforzi nelle aste

Figura 1.5. Deformata della reticolare B.

Figura 1.6. Diagramma degli sforzi normali della reticolare B (verde = positivo = trazione; rosso = negativo = compressione).

Matrice di rigidezza

La matrice di rigidezza globale (K) compare all’interno di un’equazione lineare K a = f. Essa lega in una relazione di equilibrio le forze esterne (f) con la variabile configurazione a (spostamenti nodali).

Noto il vettore degli spostamenti nodali si ricava il vettore delle forze che produce quegli spostamenti. Nella pratica noto f si cerca la configurazione a che ne consegue. Per quanto concerne una struttura priva di vincoli esterni, il determinante della matrice di rigidezza è nullo e quindi la soluzione non è unica. Questo è dovuto alla presenza di moti di corpo rigido; per ovviare a tale problema occorre impedire tali moti, ovvero introdurre delle condizioni al contorno per le variabili di configurazione a, ossia occorre vincolare la struttura in modo tale da impedire tali moti.

Imponendo quindi un vincolo agli spostamenti, si otterrà una matrice definita positiva, simmetrica, con determinante diverso da zero e quindi invertibile. Risolvendo il sistema così ottenuto si otterranno i valori degli spostamenti relativi agli altri nodi e la reazione vincolare.

La struttura analizzata negli esercizi 1 e 2 è costituita da 8 nodi e, poiché si sta lavorando nel piano, la dimensione della matrice di rigidezza è 16x16 in quanto, nel piano, un punto possiede due gradi di libertà.

Per entrambe le disposizioni dei nodi, la matrice K risulta essere simmetrica. Dalla Figura 1.7 è possibile effettuare il confronto della disposizione dei valori all’interno delle matrici di rigidezza.

Figura 1.7. Distribuzione dei valori di K della reticolare A (sinistra) e della reticolare B (destra).

Confrontando la matrice di rigidezza ottenuta con la numerazione del caso A e quella del caso B si può osservare come nel primo caso questa presenti termini non nulli nelle vicinanze della diagonale principale, determinando una matrice approssimabile con una diagonale a bande. Viceversa, con la seconda numerazione, la matrice presenta dei valori non nulli anche alle estremità. L’utilizzo di una matrice diagonale a bande rende più agevole l’operazione di fattorizzazione, per tanto la numerazione della struttura A ottimizza considerevolmente la matrice delle rigidezze. In particolare è sempre bene evitare che un’asta colleghi due nodi con una numerazione molto differente, in quanto tali elementi andrebbero ad occupare i posti più esterni della matrice.

Struttura reticolare B con spostamento imposto

Imponendo ora di applicare uno spostamento in direzione y, con verso positivo e modulo pari a 0.5 m nel nodo numero 5 (quello centrale), risulta necessario imporre un vincolo nel nodo per effettuare le operazioni. In questo caso il numero di vincoli passa a 4 e la matrice degli spostamenti Uimp non avrà più tutti gli elementi nulli, ma nella posizione (4,1), in corrispondenza del vincolo aggiunto, presenterà il valore di spostamento richiesto.

La geometria della struttura, le forze applicate e la rigidità delle aste risultano le stesse della reticolare B. Si è risolta la reticolare in maniera analoga ai casi precedenti ottenendo i seguenti risultati:

Spostamenti nodali

Reazioni vincolari

Sforzi nelle aste

Figura 1.8. Deformata della reticolare B con spostamento imposto.

Figura 1.9. Diagramma degli sforzi normali della reticolare B con spostamento imposto (verde = positivo = trazione; rosso = negativo = compressione).

Struttura reticolare C

Si analizza ora la struttura in Figura 1.10 che si differenzia dalla struttura precedente dall’aggiunta di due aste al modello, rispettivamente l’asta 8-5 e l’asta 6-5. Essa sarà sottoposta alle medesime condizioni di carico delle precedenti.

Figura 1.10. Struttura reticolare C.

Spostamenti nodali

Reazioni vincolari

Sforzi nelle aste

Figura 1.11. Deformata della reticolare C.

Figura 1.12. Diagramma degli sforzi normali della reticolare C (verde = positivo = trazione; rosso = negativo = compressione).

Struttura reticolare C con spostamento imposto

Come nel caso della reticolare B, si procede imponendo alla struttura uno spostamento del nodo 5 in direzione y pari a 0.5 m, ottenendo i seguenti risultati:

Spostamenti nodali

Reazioni vincolari

Sforzi nelle aste

Figura 1.13. Deformata della reticolare C con spostamento imposto.

Figura 1.14. Diagramma degli sforzi normali della reticolare C con spostamento imposto (verde = positivo = trazione; rosso = negativo = compressione).

Come è stato osservato nel confronto della reticolare B, tra la struttura a cui è stato applicato lo spostamento nel nodo 5 e la struttura senza imposizioni, anche nella struttura C si possono effettuare le stesse considerazioni. Lo spostamento imposto risulta essere uno spostamento notevole se confrontato con la dimensione della struttura, infatti, le deformazioni tendono a variare quasi tutte di un ordine di grandezza.

Ovviamente, siccome lo spostamento è applicato in direzione opposta alle forze esterne, si osserva anche un cambio di direzione degli spostamenti riferiti ai nodi. Stessa osservazione viene fatta per gli sforzi nelle aste. Anche in questo caso si osserva un'inversione di segno nelle varie aste che ne comporta una variazione del comportamento: da trazione a compressione e viceversa.

Inoltre, confrontando i risultati ottenuti con quelli del caso della reticolare B, si nota come la struttura sia molto più rigida quando vengono aggiunti due diagonali nella parte centrale: in questo caso gli spostamenti nodali risultano essere più esigui, facendo deformare meno la struttura. Per quanto riguarda le sollecitazioni invece, a meno delle aste più esterne, che presentano lo stesso valore di sollecitazione, le aste interne presentano delle sollecitazioni maggiori nel caso B, ossia quando si modellano solamente due diagonali.

Struttura reticolare A con spostamenti noti

In questo caso si risolve la reticolare A, andando a rimuovere i carichi applicati lungo la corrente inferiore e imponendo nei nodi gli spostamenti ottenuti nella medesima. Per effettuare tali operazioni è necessario imporre dei vincoli fittizi ad ogni nodo, perciò il numero complessivo delle componenti di spostamento vincolate cresce a 16 e la matrice degli spostamenti Uimp risulterà non nulla, in quanto in corrispondenza di ogni vincolo fittizio aggiunto, presenterà il valore di spostamento richiesto. La geometria della struttura e la rigidità delle aste risulteranno invece le stesse del primo studio effettuato sulla reticolare A.

Risolvendola si sono ottenuti i seguenti risultati:

Spostamenti nodali

Nodo Direzione Spostamento [m]

Reazioni vincolari

Nodo Direzione Reazione vincolare [kN]