ESERCITAZIONE Con Confronto Matlab-SAP2000

1

ESERCITAZIONE con confronto Matlab/SAP2000 Corso di Meccanica Computazionale delle Strutture Professore: Prof. Ing. M. Gioffrè Studente: Alessandro Cancelli

Università degli studi di Perugina. Facoltà di Ingegneria.

2

1 Introduzione

La presente relazione tecnica descrive le procedure eseguite nella risoluzione di un telaio piano costituito da profilati in acciaio e sottoposto all’azione di forze sia di natura statica sia dinamica. L’analisi, effettuata mediante il programma Matlab R2010b, ha riguardato lo studio del problema statico e del problema dinamico; i risultati ottenuti sono stati confrontati con quelli conseguiti risolvendo la struttura con il programma di calcolo SAP 2000 v14. In prima istanza è stata eseguita l’analisi statica sottoponendo la struttura a carichi statici. Una volta determinata la matrice di rigidezza della struttura e definiti i carichi, sono stati valutati gli spostamenti dei nodi non vincolati e determinate le reazioni vincolari. Quindi è stato possibile plottare la deformata della struttura considerando gli spostamenti dei nodi non vincolati. Sono stati poi calcolati gli spostamenti a livello locale dei singoli nodi nonché le sollecitazioni in tutte le aste che compongono il telaio. Infine sono stati plottati i diagrammi delle sollecitazioni di sforzo normale, taglio e momento flettente asta per asta. Il problema dinamico è stato affrontato analizzando il comportamento della struttura sia in oscillazioni libere, con smorzamento nullo, sia in oscillazioni forzate, causate da una forzante di tipo sinusoidale. Lo studio della struttura in oscillazioni libere è stato eseguito considerando quali condizioni iniziali gli spostamenti dei nodi ottenuti dall’analisi statica. L’analisi in oscillazioni forzate con smorzamento nullo è stata affrontata introducendo una forzane armonica di assegnata frequenza ed ampiezza. Tale analisi è stata eseguita mediante i Metodi di integrazione diretta: Metodo alle differenze centrali, Metodo di Hubolt, Metodo di Wilson-, Metodo di Newmark ed il Metodo della Sovrapposizione Modale, mettendo in evidenza le differenze riscontrate fra le suddette procedure di calcolo. Al termine si andrà inoltre ad analizzare una platea di fondazione in c.a. mediante l’utilizzo di elementi shell, cioè elementi finiti piani, all’interno del codice di calcolo SAP2000. In particolare verranno utilizzati elementi tipo membrane, che sono adatti a simulare il comportamento a piastra della platea di fondazione in esame, sia di tipo thin che di tipo thick e si discuterà poi su quale sia la scelta migliore fra i due.

3

2 Analisi Statica 2.1 Definizione della geometria della struttura La struttura oggetto di studio è un telaio piano in acciaio [ E  2.1  10 8 KN m 2 ] avente la seguente conformazione e soggetto ai seguenti carichi:

La struttura è caratterizzata esternamente da tre vincoli a incastro, mentre internamente sono presenti anche delle sconnessioni a cerniera.

2.2 Discretizzazione Definita la geometria del sistema si procede alla sua discretizzazione, individuando i nodi della struttura e i conseguenti elementi finiti presenti fra questi. Si utilizzano due tipi di elementi finiti frame: - gli elementi tipo ‘trave’ che nel piano hanno 6 gradi di libertà, 3 per ogni nodo, che sono traslazione lungo l’asse dell’elemento, traslazione in direzione ortogonale all’asse e la rotazione nel piano; - gli elementi di tipo ‘biella’ che nel piano hanno 2 gradi di libertà, 1 per ogni nodo, ed è solamente la traslazione lungo l’asse dell’elemento. Quindi si è suddivisa la struttura in 8 nodi, individuando 13 aste; tali nodi e aste vanno opportunamente numerate in modo da consentire l’acquisizione dei dati geometrici in ambiente Matlab. Inoltre un’opportuna numerazione diventa molto importante per curare due aspetti fondamentali: - Consente di individuare il sistema di riferimento dell’elemento finito, fissandone il nodo iniziale e quello finale e consentendo quindi una corretta interpretazione dei risultati sia in termini di spostamenti che in termini di sollecitazioni. È quindi opportuno che l’orientamento del sistema di riferimento degli elementi finiti sia omogeneo in tutta la struttura (ad esempio sempre dal basso verso l’alto e sempre da sinistra verso destra)

4 -

Il secondo aspetto riguarda l’ampiezza della banda della matrice di rigidezza della struttura Ks, più precisamente un’opportuna numerazione consente di far addensare i valori intorno alla diagonale principale di tale matrice; ciò semplifica notevolmente le elaborazioni svolte dal calcolatore.

Per quanto riguarda la struttura in esame il risultato della discretizzazione può essere osservato nella figura precedente. All’interno del programma Matlab l’inserimento delle coordinate dei nodi avviene attraverso la definizione della matrice Mxy (=matrice delle coordinate x e y) nella quale sono inserite le coordinate x e y dei vari nodi in due righe separate, seguendo l’ordine della loro numerazione:

Con il comando di riga 11 si è definita la grandezza n che rappresenta il numero di nodi della struttura.

2.3 Collegamenti Una volta definita la discretizzazione della struttura, cioè il numero di nodi e la loro posizione, è necessario andare a definire come questi nodi sono collegati dai vari elementi computazionali; per fare ciò in ambiente Matlab si è definita una matrice M_coll (=matrice dei collegamenti) in cui si pone: 1 se i due nodi sono collegati da una ‘trave’, 2 se i due nodi sono collegati da una ‘biella’, 0 se i due nodi non sono collegati.

5

2.4 Attribuzione sezioni e Matrice riassuntiva Ora che abbiamo definito sia i nodi che le aste della struttura possiamo passare all’assegnazione delle sezioni strutturali alle aste stesse. Per la nostra struttura si decide di utilizzare un HEB 240 per quello che riguarda i pilastri, un IPE 240 per quello che riguarda le travi e un profilato angolare accoppiato 2Lx100x100x10 per quel che riguarda i controventi. Per definire la caratteristiche delle varie sezioni in Matlab si sono create le 3 variabili ‘IPE240_valori’, ‘HEB240_valori’, ‘doppia_L_valori’ mediante lettura delle caratteristiche da un file Excel (preventivamente realizzato) in cui sono racchiuse le caratteristiche di molteplici tipologie di sezioni; di seguito si riporta la parte di codice che realizza ciò:

A questo punto nel programma si definisce una matrice riassuntiva degli elementi chiamata M_ele_riass in cui si elencano gli elementi (ordinando gli elementi a partire dai pilastri, poi le travi e poi i controventi) e ad ognuno di essi si associa (nell’ordine): -tipo elemento, -identificativo dell’elemento, -identificativo del primo nodo dell’elemento, -coordinata x del primo nodo xi, -coordinata y del primo nodo yi, -identificativo del secondo nodo dell’elemento, -coordinata x del secondo nodo xj, -coordinata y del secondo nodo yj, -lunghezza dell’elemento l, -angolo di inclinazione  [cioè l’angolo che il sistema di riferimento locale dell’elemento (che da adesso verrà indicato con SRL) forma con il sistema di riferimento globale fissato per la struttura(che da adesso indicheremo con SRG)], -modulo elastico E, -area della sezione A, -momento di inerzia intorno a y Iy; -momento di inerzia intorno ad x Ix. Di seguito si riporta la parte di codice Matlab necessaria a creare tale matrice (compresa la parte necessaria al calcolo di le):

6

7

8

2.5 Calcolo delle rigidezze delle singole aste. Come già accennato in precedenza la struttura presa in esame comprende al suo interno sia elementi tipo trave che elementi tipo biella le cui matrici di rigidezza verranno trattate in maniera separata. Per quanto riguarda la matrice di rigidezza degli elementi trave dobbiamo ricordarci che tali elementi (nel piano) offrono una rigidezza a 6 gradi di libertà, 3 per ogni nodo, e perciò la loro matrice di rigidezza sarà una 6x6. Tale matrice si ottiene a partire dall’elemento nel sistema geometricamente determinato a nodi bloccati, da qui si sbloccano uno alla volta i vari gradi di libertà e ogni volta si assegna al grado di libertà sbloccato uno spostamento unitario con il quale si vanno a determinare tutte le forze che insorgono nell’elemento a causa di tale spostamento; proprio tali forze sono le componenti di rigidezza della trave per tale spostamento unitario. Quindi il generico elemento K ij della matrice di

rigidezza dell’elemento è la forza in direzione dell’i-esimo grado di libertà provocata da uno spostamento unitario in direzione del j-esimo grado di libertà; inoltre per il Teorema di Betti risulta che K ij  K ji quindi tali matrici risulteranno simmetriche. Quindi per l’elemento trave la matrice di rigidezza è la seguente:  EA  l   0   0  [K ' ]   EA   l  0    0 

-

0

0

12EJ l3 6EJ l2

6EJ l2 4EJ l

0

0

12EJ l3 6EJ l2

-

-

6EJ l2 2EJ l

EA l 0 0 EA l 0 0

0 12EJ l3 6EJ - 2 l

-

0 12EJ l3 6EJ - 2 l

 0  6EJ  - 2  l  2EJ  l   0   6EJ  - 2 l  4EJ   l 

A differenza della trave, l’elemento biella (nel piano) ha 2 gradi di libertà per i quali offre rigidezza (cioè le due traslazioni assiali ai due nodi) quindi la sua matrice di rigidezza sarà una 2x2, e avrà la seguente forma:

K   EAl  - 11 '



-1   1

In ambiente Matlab si è realizzata un apposita parte di codice che a seconda del particolare tipo di elemento computazionale crea la relativa matrice di rigidezza; per fare ciò si utilizza il comando switch che ci permette di passare da una caso ad un altro a seconda del valore che assume una determinata variabile, più precisamente in tale caso quando la variabile hh è ‘pilastro’ o ‘trave’ all’elemento viene associata la matrice di rigidezza dell’elemento computazionale trave, mentre se tale variabile è ‘controvento’ all’elemento viene assegnata la matrice di rigidezza dell’elemento computazionale biella. Di seguito si riporta la parte di codice relativa a quanto appena descritto:

9

Tuttavia tali matrici di rigidezza così calcolate sono riferite al SRL dell’elemento, per riportarle nel SRG è necessario utilizzare la matrice di trasposizione (che non è altro che una matrice di cambiamento di base) L che è definita a partire dai coseni direttori dell’elemento che sono n x  cos  n y  cos (



-  )  sin 

n z  1 , dove  è l’angolo fra l’asse x del SRG e l’asse x '

2 (cioè quello lungo l’asse dell’elemento) del SRL.

10

Più precisamente la matrice di trasposizione per un elemento trave ha la seguente forma:  nx  - n y  0 L    0   0  0

ny

0

0

0

nx

0

0

0

0

1

0

0

0

0

nx

ny

0

0

- ny

nx

0

0

0

0

0  cos    0 - sin  0  0  0  0   0  0 1  0

sin 

0

0

0

cos 

0

0

0

0

1

0

0

0

0

cos 

sin 

0

0

- sin 

cos 

0

0

0

0

0  0 0  0  0 1

Invece la matrice di trasposizione per l’elemento biella ha la seguente forma: n x

ny

0

0

0

nx

L  

0  cos   n y   0

sin 

0

0

cos 

0   sin  

Per calcolare tali matrici in ambiente Matlab è stato previsto un apposito ciclo sempre utilizzando il comando switch come visto sopra per le matrici di rigidezza locali; di seguito si riporta tale parte di codice:

11

Ora che abbiamo a disposizione le matrici di trasposizione possiamo passare le matrici di rigidezza del singolo elemento dal SRL al SRG mediante la seguante relazione:

K   L  K  L T

glob

loc

Nell’ambiente Matlab tale relazione è implementata in appositi cicli che calcolano le matrici di rigidezza dei vari elementi nel SRG e li salvano in delle matrici tridimensionali; di seguito si riporta la parte di codice relativa:

2.6 Ripristino della congruenza e determinazione della matrice delle rigidezze della struttura è necessario ora ripristinare la congruenza, poiché finora sono state sempre prese in considerazione le singole aste in modo separato. Si procede quindi alla determinazione della matrice topologica [A]; essa è costituita da 0 e da 1 e serve a correlare gli spostamenti nodali dei singoli elementi agli spostamenti globali dei nodi della struttura secondo la relazione: u  [A]  u s

Dove: - u è il vettore degli spostamenti nodali delle singole aste ed ha dimensioni [(d 1  1)] dove d 1  (m1  n 1  p1 )  (m 2  n 1  p 2 ) con: m 1  numero delle travi, m 2  numero delle bielle, n 1  numero di nodi x elemento, p1  gdl del singolo nodo della trave, p 2  gdl del singolo nodo della biella - u s è il vettore degli spostamenti globali dei nodi della struttura ed ha dimensioni [(d 2  1)] dove d 2  n  p con: n = numero di nodi della struttura, p = numero di gdl del singolo nodo della struttura

12 - [A] è la matrice topologica che quindi avrà dimensioni [(d 1  d 2 )] dove d 1 e d 2 hanno il significato visto precedentemente Si può quindi ricavare la matrice di rigidezza della struttura [K s ] tramite la seguente razione: [K s ]  [A ]T  [K ]  [A ]

Dove con [K ] si è indicata la matrice a banda contenete tutte le matrici di rigidezza dei singoli elementi nel sistema di riferimento globale, di dimensioni [d 1  d 1 ] , cioè:

K 1  [K ]    0 

  

0     K m 

In ambiente Matlab per fare ciò si è prima creata la matrice [K ] , si è poi creata la matrice topologica mediante un apposito ciclo che assegna gli uno e gli zero della matrice a in base all’identificativo dei nodi iniziali e finali dell’elemento, in modo da garantire la congruenza fra spostamenti locali dell’elemento e spostamenti globali della struttura; si riportano di seguito le parti di codice relative:

13

14

2.7 Introduzione dei vincoli I vincoli vengono inizialmente introdotti mediante il vettore ‘tipo_vincolo’ in cui si inserisce 0 se il nodo non è vincolato, 1 se il nodo è incastrato, 2 se il nodo è incernierato (sempre all’esterno); dopodichè tali condizioni di vincolo vengono salvate all’interno della variabile num_punti, in cui sono contenuti (in ordine): - identificativo del nodo, -coordinata x del nodo, - coordinata y del nodo, - condizione di vincolo del nodo. Grazie a questa grandezza è possibile definire quali gradi di libertà sono vincolati attraverso un apposito ciclo che tramite la condizione di vincolo e l’identificativo del punto vincolato calcola il numero dei gdl vincolati (ad esempio se la condizione di vincolo è 1 e l’identificativo del punto è 2 i numeri dei gdl vincolati sarà 4,5 e 6) e li inserisce tutti in un apposita variabile chiamata ‘gdlr’ (=gradi di libertà vincolati); dopodichè vengono calcolati i gradi di libertà non vincolati togliendo ai gdl totali della struttura quelli vincolati già calcolati. Di seguito si riportano le varie parti di codice con cui si è realizzato quanto esposto poc’anzi:

15

2.8 Inserimento dei carichi Per quanto riguarda i carichi essi vengono definiti mediante l’ausilio di due variabili, che sono: - F che è il vettore delle forzanti applicate nei nodi della struttura - q che è il vettore dei carichi distribuiti applicati nei vari elementi della struttura Una volta definite tali variabili è necessario ricordarsi che per eseguire il calcolo mediante il metodo degli spostamenti è necessario togliere al vettore dei carichi esterni il vettore delle forze del sistema geometricamente determinato a nodi bloccati; per fare ciò si sono impostati due cicli che trasformano gli eventuali carichi distribuiti presenti su travi e pilastri nelle equivalenti forze nodali che sono le reazioni vincolari dell’elemento calcolate con lo schema incastro-incastro, come riportato nella figura sottostante:

Una volta determinate tali forze per ogni elemento esse vanno sottratte al vettore dei carichi esterni così da trovare il vettore dei carichi necessario alla risoluzione della struttura. Quanto sopra riportato è stato implementato in Matlab tramite una serie di istruzioni che vengono di seguito riportate:

16

17

18

2.9 Determinazione degli spostamenti ai nodi non vincolati e reazioni vincolari Ora che abbiamo definito sia la matrice di rigidezza che i carichi agenti sul sistema ci ricordiamo che per una generica struttura vale la relazione:

Fs   K s    s Però in tale relazione sorge il problema che K s  è una matrice singolare, quindi non possiamo invertirla, poiché ancora non abbiamo usato la definizione dei vincoli per eliminare i loro gdl dal calcolo. Per fare ciò si divide il vettore spostamenti in due sottovettori  r e  f che rappresentano rispettivamente gli spostamenti dei nodi vincolati (che sono già noti essendo nulli nel caso di vincoli perfetti) e gli spostamenti dei nodi non vincolati (che sono incogniti), e poi si suddivide la matrice di rigidezza nelle seguenti quattro sottomatrici:

Ks rr 

-

che è la sottomatrice ottenuta eliminando da K s  le righe e le colonne relative ai gdl non vincolati

-

Ks ff  che è la sottomatrice ottenuta eliminando da K s  le righe e le colonne relative ai gdl vincolati

Ks rf 

-

che è la sottomatrice ottenuta eliminando da K s  le righe relative ai gdl non vincolati e le colonne relative ai gdl vincolati

-

Ks fr  che è la sottomatrice ottenuta eliminando da K s  le righe relative ai gdl vincolati e le colonne relative ai gdl non vincolati

Quindi riarrangiando l’espressione sopra riportata secondo questa nuova suddivisione si ottiene che: Fr  Ks rr   Ff  Ks fr

Ks rf   r    Ks ff   f 

Che è l’espressione matriciale di due equazioni, se lo esplicitiamo otteniamo che : Fr   Ks rr    r   Ks rf    f   Ff   Ks fr    r   Ks ff    f 

Dove Fr = forze nei nodi vincolati (che sono le reazioni vincolari) e Ff = forze nei nodi non vincolati; ora ricordandoci che gli spostamenti dei gdl vincolati sono tutti nulli, cioè  r   0 , avremo che gli spostamenti nei nodi non vincolati sono dati dalla relazione:

 f   Ks ff 1  Ff  In ambiente Matlab si è prima provveduto alla suddivisione di K s  nelle quattro sottomatrici e di F nei due sottovettori mediante le variabili gdlr e gdlf prima definite:

19

Una volta fatto ciò è semplice determinare gli spostamenti dei nodi non vincolati  f  mediante la semplice relazione sopra esposta; dopodichè si possono anche facilmente calcolare le reazioni vincolari tramite la relazione

Fr   Ks rf    f  Si riporta qui di seguito la relativa parte di codice Matlab:

2.10 Calcolo delle sollecitazioni interne Per calcolare le sollecitazioni su ogni singolo elemento è necessario per prima cosa conoscere gli spostamenti dei nodi del singolo elemento nel SRG così da poter usare la relazione:

 

Fi'  K i'   i' Per determinare tali spostamenti del singolo elemento per prima cosa si crea un vettore contenete gli spostamenti dei gdl liberi nelle posizioni dei gdl liberi e gli spostamenti dei gdl vincolati (che ricordiamo essere 0) nelle posizioni dei gdl vincolati che viene poi moltiplicato per la matrice topologica [A] così da ottenere un vettore in cui sono racchiusi gli spostamenti dei singoli elementi nel sistema di riferimento globale; infine si estraggono di volta in volta gli spostamenti di un elemento da tale vettore, li si moltiplica per la matrice di trasposizione [L] dell’elemento così da riportarli nel SRL e poi si salvano in un vettore in cui sono contenuti tutti questi spostamenti (elemento per elemento) nel SRL. Tutto ciò è implementato mediante le seguenti istruzioni Matlab:

20

Una volta noti tali spostamenti possiamo andare a determinare le sollecitazioni sui vari elementi mediante la relazione sopra descritta, ricordandoci però che per il calcolo degli spostamenti si erano tolte le forze del sistema a nodi bloccati che ora dovranno essere risommate alle sollecitazioni così calcolate per ottenere i valori finali delle sollecitazioni stesse; di seguito si riporta la parte di codice Matlab che realizza il calcolo di tali sollecitazioni:

21

Ora non ci rimane che andare a cambiare i segni delle sollecitazioni, poiché quelle appena calcolate hanno il segno riferito al sistema di riferimento locale dell’elemento mentre, per una pura semplicità di esposizione,a noi fa più comodo riferirlo alla convenzioni dei segni di scienza delle costruzioni; di seguito si riporta la parte di codice che effettua questo cambiamento di segno:

22

2.11 Confronto grafico-numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di spostamento Ora andremo ad esporre i risultati si qui ottenuti con il codice Matlab e li confronteremo con quelli ottenuti utilizzando il codice di calcolo SAP2000 sia in maniera grafica che in maniera numerica. Per prima cosa si riporta la parte di codice che in Matlab disegna lo schema della struttura comprensivo di identificativi grafici quali: - dei quadrati rossi ove siano presenti vincoli esterni (che siano essi cerniere o incastri), - dei pallini gialli in corrispondenza delle varie sconnessioni interne a cerniera, - dei puntini blu per identificare la presenza di un nodo generico; vengono inoltre stampati a video gli identificativi numerici dei nodi e degli elementi.

23

24

25

Il risultato che si ottiene da questa parte di codice è la seguente figura che rappresenta lo schema della struttura in esame:

26

Adesso si va a riportare la parte di codice che a partire dagli spostamenti dei nodi dei vari elementi nel sistema di riferimento globale prima calcola le costanti della linea elastica elemento per elemento (mediante gli spostamenti ortogonali all’asse dell’elemento e alle rotazioni, nei nodi iniziali e finali) e poi con tali costanti va a disegnare la deformata di ogni elemento così da disegnare la deformata globale della struttura; tale deformata ottenuta tramite Matlab verrà poi confrontata con quella ottenuta dal codice di calcolo SAP2000.

27

28

29

30

31

32 La deformata che si ottiene mediante tale parte di codice Matlab è la seguente:

Invece la deformata che si ottiene mediante il codice SAP2000 è la seguente:

33 Come si vede le due deformate sono uguali, però per un raffronto più preciso di seguito si riporta una tabella con i valori numerici degli spostamenti nodali calcolati con il Matlab, una contenente quelli calcolati con il SAP2000 e una con il confronto fra i due calcoli in termini percentuali:

Joint OutputCase Text Text 1 COMBO 2 COMBO 3 COMBO 4 COMBO 5 COMBO 6 COMBO 7 COMBO 8 COMBO

TABLE: CaseType Text Combination Combination Combination Combination Combination Combination Combination Combination

Joint Text 1 2 3 4 5 6 7 8

TABLE: Joint Displacements Matlab CaseType U1 U2 U3 R1 R2 Text m m m Radians Radians 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,002063 -0,000497 0 0 0 0,001894 -0,000974 0 0 0 0,001578 -0,000378 0 0 0 0,008058 -0,000781 0 0 0 0,006866 -0,001325 0 0 0

OutputCase Text COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO n° nodo / 4

5

6

7

8

spostamento / U1 U2 R2 U1 U2 R2 U1 U2 R2 U1 U2 R2 U1 U2 R2 R

Joint Displacements SAP2000 U1 U2 U3 R1 R2 R3 m m m Radians Radians Radians 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0,002063 -0,000497 0 0 0 -0,002367 0,001894 -0,000974 0 0 0 -0,002304 0,001578 -0,000378 0 0 0 0,004233 0,008058 -0,000781 0 0 0 -0,008587 0,006866 -0,001325 0 0 0 0,007164

Matlab / 0,002063 -0,000497 -0,002367 0,001894 -0,000974 -0,002304 0,001578 -0,000378 0,004233 0,008058 -0,000781 -0,008587 0,006866 -0,001325 0,007164

SAP2000 / 0,002063 -0,000497 -0,002367 0,001894 -0,000974 -0,002304 0,001578 -0,000378 0,004233 0,008058 -0,000781 -0,008587 0,006866 -0,001325 0,007164

diff. Percentuale (%) -0,016 0,030 -0,017 -0,022 0,004 0,002 -0,026 0,065 0,008 -0,002 -0,031 0,002 0,005 -0,004 0,005

R3 Radians 0 0 0 -0,002367 -0,002304 0,004233 -0,008587 0,007164

34

2.12 Confronto grafico-numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di sollecitazioni interne Si va ora a riportare la parte di codice Matlab che disegna i diagrammi delle sollecitazioni interne (che sono già state calcolate al punto 2.10) mediante il comando fill che permette di riempire una data area con un determinato colore, in particolare per i nostri diagrammi quando la sollecitazione interna (sia che sia N, T o M) è positiva il diagramma viene riempito con il colore blu, invece quando essa è negativa viene riempita con il colore rosso; la definizione di positivo e negativo fa riferimento alla convenzione di scienza delle costruzioni cioè:

Si consideri che le fibre di riferimento sono sempre poste alla destra dell’elemento che viene guardato a partire dal nodo i fino al nodo j. Partiamo dal codice per il plot del diagramma dello sforzo normale:

35

36

37 Ora andiamo a riportare il diagramma così ottenuto e lo confrontiamo con quello ottenuto mediante il codice SAP2000: MATLAB

SAP2000

38 Adesso verrà riportata la parte di codice che permette di effettuare il plot del diagramma del taglio:

39

Ora andiamo a riportare il diagramma così ottenuto e lo confrontiamo con quello ottenuto mediante il codice SAP2000:

40 MATLAB

SAP2000

41 Infine si riporta la parte di codice che permette di effettuare il plot del diagramma del momento:

42

43

44

45

Ora andiamo a riportare il diagramma così ottenuto e lo confrontiamo con quello ottenuto mediante il codice SAP2000: MATLAB

46 SAP2000

Come già fatto per la deformata, per avere un raffronto più preciso fra i risultati, si riporta una tabella contenente i valori delle sollecitazioni, elemento per elemento, ottenute mediante il SAP200, una contenente le sollecitazioni ottenute con il Matlab e una con le loro differenze percentuali: TABLE: Element Forces – Frames SAP2000 Frame Station OutputCase CaseType P V2 V3 Text m Text Text KN KN KN 1 0,000 COMBO Combination -221,199 -8,751 0,000 1 5,000 COMBO Combination -221,199 -8,751 0,000 2 0,000 COMBO Combination -433,606 -8,775 0,000 2 5,000 COMBO Combination -433,606 -8,775 0,000 3 0,000 COMBO Combination -168,177 27,604 0,000 3 5,000 COMBO Combination -168,177 27,604 0,000 4 0,000 COMBO Combination -158,261 -70,555 0,000 4 4,000 COMBO Combination -158,261 -70,555 0,000 5 0,000 COMBO Combination -195,381 65,133 0,000 5 4,000 COMBO Combination -195,381 65,133 0,000 6 0,000 COMBO Combination -18,494 117,913 0,000 6 7,500 COMBO Combination -18,494 -125,837 0,000 7 0,000 COMBO Combination -43,244 129,857 0,000 7 6,000 COMBO Combination -43,244 -125,143 0,000

T KN-m 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000

M2 KN-m 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000

M3 KN-m 10,683 -33,074 11,041 -32,833 -48,992 89,026 104,344 -177,876 -74,298 186,233 -137,419 -167,131 -125,666 -111,526

47 8 8 9 9 10 10 11 11 12 12 13 13

0,000 7,500 0,000 9,014 0,000 9,014 0,000 7,810 0,000 7,810 0,000 8,500

Frame Station Text m 1 0,000 1 5,000 2 0,000 2 5,000 3 0,000 3 5,000 4 0,000 4 4,000 5 0,000 5 4,000 6 0,000 6 7,500 7 0,000 7 6,000 8 0,000 8 7,500 9 0,000 9 9,014 10 0,000 10 9,014 11 0,000 11 7,810 12 0,000 12 7,810 13 0,000 13 8,500

COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO COMBO

Combination -130,555 -158,261 0,000 0,000 0,000 -177,876 Combination -130,555 160,489 0,000 0,000 0,000 -186,233 Combination 18,826 0,000 0,000 0,000 0,000 0,000 Combination 18,826 0,000 0,000 0,000 0,000 0,000 Combination -36,206 0,000 0,000 0,000 0,000 0,000 Combination -36,206 0,000 0,000 0,000 0,000 0,000 Combination 20,359 0,000 0,000 0,000 0,000 0,000 Combination 20,359 0,000 0,000 0,000 0,000 0,000 Combination -43,598 0,000 0,000 0,000 0,000 0,000 Combination -43,598 0,000 0,000 0,000 0,000 0,000 Combination 74,145 0,000 0,000 0,000 0,000 0,000 Combination 74,145 0,000 0,000 0,000 0,000 0,000

TABLE: Element Forces - Frames Matlab OutputCase CaseType P V2 V3 Text Text KN KN KN COMBO -221,199 -8,751 0,000 COMBO -221,199 -8,751 0,000 COMBO -433,606 -8,775 0,000 COMBO -433,606 -8,775 0,000 COMBO -168,177 27,604 0,000 COMBO -168,177 27,604 0,000 COMBO -158,261 -70,555 0,000 COMBO -158,261 -70,555 0,000 COMBO -195,381 65,133 0,000 COMBO -195,381 65,133 0,000 COMBO -18,494 117,913 0,000 COMBO -18,494 -125,837 0,000 COMBO -43,244 129,857 0,000 COMBO -43,244 -125,143 0,000 COMBO -130,555 158,261 0,000 COMBO -130,555 -160,489 0,000 COMBO 18,826 0,000 0,000 COMBO 18,826 0,000 0,000 COMBO -36,206 0,000 0,000 COMBO -36,206 0,000 0,000 COMBO 20,359 0,000 0,000 COMBO 20,359 0,000 0,000 COMBO -43,598 0,000 0,000 COMBO -43,598 0,000 0,000 COMBO 74,145 0,000 0,000 COMBO 74,145 0,000 0,000

T KN-m 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000

M2 KN-m 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000

M3 KN-m 10,682 -33,074 11,041 -32,8326 -48,992 89,026 104,344 -177,876 -74,298 186,233 -137,419 -167,131 -125,666 -111,526 -177,876 -186,233 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000 0,000

48

CONFRONTO SOLLECITAZIONI MATLAB-SAP2000 ELEMENTO STAZIONE P V  M / / (%) (%) (%) 0,000 0,000090 -0,004571 0,000936 1 5,000 0,000090 -0,004571 0,000000 0,000 -0,000046 0,003419 -0,000906 2 0,003419 -0,000305 5,000 -0,000046 0,000 0,000119 0,001449 0,000000 3 5,000 0,000119 0,001449 0,000000 0,000 0,000190 0,000000 0,000000 4 4,000 0,000190 0,000000 0,000000 0,000 -0,000102 0,000461 -0,000135 5 4,000 -0,000102 0,000461 0,000000 0,000 0,002163 -0,000339 0,000000 6 7,500 0,002163 0,000318 0,000000 0,000 0,000925 0,000308 0,000000 7 6,000 0,000925 -0,000320 0,000000 0,000 0,000000 0,000190 0,000000 8 -0,000187 0,000000 7,500 0,000000 0,000 0,000000 0,000000 0,000000 9 9,014 0,000000 0,000000 0,000000 0,000 0,000276 0,000000 0,000000 10 9,014 0,000276 0,000000 0,000000 0,000 0,000982 0,000000 0,000000 11 7,810 0,000982 0,000000 0,000000 0,000 -0,000688 0,000000 0,000000 12 7,810 -0,000688 0,000000 0,000000 0,000 -0,000405 0,000000 0,000000 13 8,500 -0,000405 0,000000 0,000000

49

2.13 Confronto numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di reazioni vincolari Di seguito si riportano tre tabelle: - una contenente i valori delle reazioni vincolari calcolate mediante il Matlab (come visto al termine del paragrafo 2.9); - un’altra tabella contenente i valori delle reazioni vincolari calcolate mediante il SAP2000; - infine un’ultima contenente le differenze percentuali fra le due metodologie di calcolo

Joint OutputCase Text Text 1 COMBO 2 COMBO 3 COMBO

Joint Text 1 2 3

OutputCase Text COMBO COMBO COMBO

TABLE: Joint Reactions SAP2000 CaseType F1 F2 F3 Text KN KN KN Combination -6,913 210,756 0,000 Combination -36,990 440,656 0,000 Combination -61,097 196,088 0,000

M1 KN-m 0,000 0,000 0,000

M2 KN-m 0,000 0,000 0,000

TABLE: Joint Reactions SAP2000 CaseType F1 F2 F3 Text KN KN KN -6,913 210,756 0,000 -36,990 440,656 0,000 -61,097 196,088 0,000

M1 KN-m 0,000 0,000 0,000

M2 M3 KN-m KN-m 0,000 -10,682 0,000 -11,041 0,000 48,992

CONFRONTO REAZIONI MATLAB-SAP2000 NODO OutputCase F1 F2 M3 Text / (%) (%) (%) 1 COMBO 0,002992 -0,000009 0,000474 2 COMBO -0,001252 -0,000048 -0,000507 3 0,000117 -0,000064 COMBO 0,000419

M3 KN-m -10,683 -11,041 48,992

50

3 Analisi Dinamica 3.1 Costruzione della matrice delle masse della struttura L’equazione di equilibrio dinamico per un sistema elastico discreto, che è proprio il caso della struttura oggetto di studio, è la seguente: [M s ]  u  [C]  u  [K s ]  u  F( t ) Visto che si sta trattando il caso di un analisi di una struttura per cui il comportamento dello smorzamento non è ben noto ma solo approssimabile (ad esempio con una formulazione alla Rayleigh in cui [C]    [M s ]    K s  ) si è deciso di svolgere tutte le analisi considerando lo smorzamento nullo, cioè [C] = 0, così da determinare i massimi della risposta e non il reale decorso della risposta nel tempo; quindi a causa di questa ipotesi l’equazione della dinamica diventa: [M s ]  u  [K s ]  u  F( t ) Siccome la matrice di rigidezza della struttura è già nota, in quanto determinata nell’analisi statica, non ci rimane che andare a determinare la matrice di massa della struttura. Per fare ciò in generale è necessario definite per prima cosa la matrice di massa dei singoli elementi nel loro sistema di riferimento locale e poi riportarla nel sistema di riferimento globale tramite la relazione [M ele_glob ]  L  [M ele_loc ]  L T

Fatto ciò si deve andare a comporre la matrice delle masse della struttura mediante la matrice topologica [A] ( già definita per l’assemblaggio della [K s ] ) mediante la relazione:

M s   AT  [M ele _ glob ]  A Quanto sopra esposto vale nel caso generale in cui la matrice delle masse dell’elemento non sia diagonale ma sia del tipo congruente, cioè ottenuta a partire dalle funzioni di forma dell’elemento e dove quindi esistono anche termini fuori diagonale sia per gli elementi trave che per gli elementi biella; però nel nostro caso si utilizzerà una formulazione a masse concentrale (lumped mass) ove la massa dell’elemento è divisa equamente fra i due nodi dell’elemento per entrambi i gradi di libertà traslazionali del nodo stesso mentre a quello rotazionale non viene assegnata nessuna massa. Quindi in questa formulazione la matrice di massa dell’elemento verrà definita direttamente nel sistema di riferimento globale, poiché sarebbe inutile definire la matrice delle masse nel SRL e poi passarla nel SRG in quanto essa è diagonale e quindi nel passaggio di sistema di riferimento non verrebbe modificata. Quindi la forma della matrice di massa dell’elemento trave e di quella dell’elemento biella (nel SRG) sono le seguenti:



Elemento Trave M ele_glob



1  0   A  l 0   2 0  0 0

0

0

0

0

1

0

0

0

0

0

0

0

0

0

1

0

0

0

0

1

0

0

0

0

0  0 0  0  0 0

51



Elemento Biella M ele_glob



1    A  l 0  2 0  0

0

0

1

0

0

1

0

0

0  0 0  1

Tale formulazione per la matrice delle masse è stata scelta per ottenere un miglior riscontro con i risultati del programma di calcolo SAP2000 che utilizza proprio tale formulazione per la matrice delle masse. Una volta costruita la matrice delle masse però ci dobbiamo ricordare che alla risposta della struttura partecipano solamente i gdl non vincolati e perciò dovremo andare a ridurre il problema iniziale selezionando solamente le righe e le colonne della matrice di massa relative a gdl non vincolati, riducendo quindi il sistema iniziale al seguente sistema: [M s _ ff ]  u  [K s _ ff ]  u  F( t ) ff  Di seguito si riporta la parte di codice Matlab necessaria a costruire la matrice di massa della struttura (sia quella completa che quella ridotta):

52

53

3.2 Analisi in oscillazioni libere con smorzamento nullo Quello che ci proponiamo di fare adesso è andare a studiare il nostro sistema nel caso in cui non siano presenti forzanti esterne, ma le condizioni iniziali in termini di spostamento e velocità sono diverse da 0 (o entrambe o almeno una delle due), in pratica ci proponiamo di andare a risolvere il seguente sistema di equazioni differenziali: [M s _ ff ]  u  [K s _ ff ]  u  0  u (0)  u 0 u (0)  u 0 

Per risolvere tale problema si ipotizza che la soluzione sia del tipo u ( t )    f ( t ) dove f(t) è una funzione del tempo valida per tutti i gradi di libertà e   è un vettore di numeri di dimensioni pari al numero di gradi di libertà del sistema, in pratica si pensa che la soluzione abbia lo stesso andamento per ogni massa ma sia proporzionale al valore dello i relativo alla massa stessa. Ora derivando due volte questa espressione e sostituendola nell’equazione iniziale si ottiene il seguente sistema lineare omogeneo:

   M   K    0 2

s _ ff

s _ ff

Tale sistema ammette una soluzione (che non sia quella banale di    0) se e solo se è verificata la seguente relazione:





 



det   2  M s _ ff  K s _ ff  0

La risoluzione di questo problema è detta problema agli autovalori e si risolve tramite l’equazione di grado N ( = numero di gdl del sistema) che si ottiene dalla condizione sopraesposta e che ha per incognite gli N autovalori  2 ; a tali N autovalori sono associate N forme modali   della struttura (che si ottengono reinserendo uno alla volta i valori di  2 nell’equazione iniziale). Una volti noti gli autovalori si possono calcolare sia le pulsazioni proprie del sistema  i   i sia i periodi propri del sistema Ti  2    i . È importante notare che a causa della matrice delle masse utilizzata sulla diagonale della stessa ci sono dei termini nulli, tali termini corrispondono ad autovalori con valore numerico non finito (cioè  ) e quindi ad autovettori nulli, di conseguenza si avrà che il sistema è definito da un numero di 2 forme modali non nulle minore rispetto al numero di gdl del sistema, precisamente si avranno N 3 forme modali; tale numero deriva dal fatto che si deve togliere una forma modale per ogni gdl rotazionale (a cui non è assegnata alcuna massa) e quindi, essendo i gdl rotazionali un terzo di quelli totali, ne risulta che i restanti sono due terzi. In ambiente Matlab il calcolo degli autovettori e degli autovalori di effettua mediante il comando eig, mediante il quale si ottengono la matrice diagonale [aut_val] contenente gli autovalori del sistema e la matrice [aut_vett] contenente gli autovettori del sistema; una volte note tali quantità è facile calcolare le pulsazioni e i periodi propri del sistema. Di seguito si riporta la parte di codice che realizza tali operazioni: 2

54

Gli autovettori che si sono trovati sono definiti a meno di una costante moltiplicativa che può essere qualunque (poiché per gli autovettori non è importante il valore assoluto delle componenti ma è importante il rapporto fra le componenti stesse), allora normalmente tale costante viene scelta in maniera che siano soddisfatte le due seguenti relazioni:

 T  M s _ ff     I  T  K s _ ff      2  Per fare ciò basta moltiplicare l’i-esimo autovetture i per il coefficiente  i 





1 , dove m i è mi

l’i-esimo della matrice M s _ ff     M s _ ff    che è una matrice diagonale. In Matlab tale T

procedura è stata implementata nel seguente ciclo:

55

Ora che abbiamo a disposizione sia le pulsazioni naturali che le forme modali si può andare a determinare la risposta del problema in oscillazioni libere come combinazione lineare delle risposte associate ai singoli modi di vibrare e cioè: u ( t )   j11  j   j  sen ( j  t   j ) N

Dove  j e  j sono due parametri che dipendono dalle condizioni iniziali e N 1 è il numero di modi

2 N ). Il parametro  j prende il nome di contributo 3 modale , in quanto rappresenta il contributo della j-esima forma modale sulla risposta complessiva del sistema. Per analizzare le oscillazioni libere si assumo due diversi set di condizioni al contorno: -nel primo si assumeranno velocità nulle e spostamenti proporzionali alla prima forma modale; - nel secondo si assumeranno sempre velocità nulle ma spostamenti pari a quelli ottenuti nel calcolo della deformata statica. propri del sistema (come già visto è pari a

Set 1

u(0)    1  u (0)  0

Set 2

u(0)  u s  u (0)  0

Dalla soluzione del problema con il primo set di condizioni al contorno ci si aspetta di ottenere degli spostamenti sinusoidali con periodo pari al periodo del primo modo di vibrare, ciò perché se si assumono spostamenti proporzionali ad una forma modale sarà solo quella forma a dare contributo, cioè il  di quella forma modale sarà ari a 1 mentre tutti gli altri saranno nulli. Dalla soluzione del secondo set di condizioni al contorno invece ci si aspetta di ottenere una forma generica data dalla sovrapposizione dei contributi dei vari modi. Di seguito si riporta la parte di codice che realizza il calcolo delle oscillazioni libere per entrambi i set di condizioni iniziali (mentre le elaborazioni grafiche vengono riportate nell’apposito paragrafo):

56

57

58

3.3 Analisi in oscillazioni forzate con smorzamento nullo Quello che ci proponiamo di fare adesso è andare a studiare il nostro sistema nel caso in cui siano presenti forzanti esterne con andamento nel tempo di tipo sinusoidale, e che le condizioni iniziali in termini di spostamento e velocità sono entrambe pari 0, in pratica ci proponiamo di andare a risolvere il seguente sistema di equazioni differenziali: [M s _ ff ]  u  [K s _ ff ]  u  F( t )  u (0)  0 u (0)  0 

59 L’analisi è stata eseguita impostando come forzanti esterne due forze, agenti contemporaneamente, di ampiezza e periodo differente rispettivamente pari a: - F0 _ 1  50 KN , - F0 _ 2  80 KN , - T1  0.1 s , - T2  0.2 s . Queste due forzanti sono applicate rispettivamente al nodo 4 in direzione x e al nodo 7 sempre in direzione x. Di seguito si riporta la parte di codice relativa all’inserimento di tali forze:

La risoluzione di tale problema può essere effettuata sia mediante metodi diretti sia mediante metodi indiretti, ora andremo ad analizzare nello specifico entrambe le soluzioni.

3.3.1 Metodi di integrazione diretta Tali metodi sono detti diretti poiché non è prevista alcuna trasformazione delle equazioni in altre forme prima dell’integrazione numerica. L’idea base di questi metodi è quella di andare a discretizzare il dominio del tempo in un numero finito di intervalli di ampiezza pari a t e di calcolare in tali punti i valori di spostamenti velocità e accelerazioni mediante espressioni di carattere algebrico (e non più differenziale) che si ottengono inserendo nell’equazione generale del moto del sistema al posto della derivata prima e della derivata seconda degli spostamenti una loro approssimazione espressa in termini finiti (cioè mediante il t ). Quindi nella nostra analisi confronteremo 4 diversi metodi di integrazione diretta nel dominio del tempo che sono: -

Metodo delle differenze centrali

-

Metodo di Houbolt

-

Metodo di Newmark

-

Metodo di Wilson  

Concettualmente tali metodi prevedono gli stessi passi da seguire ma differiscono fra loro essenzialmente per: -

espressioni delle variazioni di velocità e accelerazione in funzione degli spostamenti

-

istante temporale in cui si effettua l’equilibrio dinamico

60 Andiamo quindi ad analizzare nel dettaglio questi quattro metodi di integrazione, considerando che per ogni metodo è stato scelto un passo di integrazione t  1  10 5 s , la scelta di tale passo è stata fatta tenendo in considerazione che esso deve essere minore sia del passo massimo per il metodo T delle differenze centrali che è pari a t max  min (che è il passo massimo oltre il quale non si può



andare altrimenti la soluzione “esplode”) sia dell’ordine di grandezza dei termini trascurabili che sono stati aggiunti sulla diagonale della matrice delle masse affinché essa possa essere numericamente invertita (altrimenti una matrice con uno 0 sulla diagonale non è invertibile) e che nel nostro caso sono pari a 1  10 3 . Sempre per ogni metodo si è deciso di studiare un intervallo temporale che va da zero a due secondi, cioè 0  t  2 s Metodo delle differenze centrali

Sia assume che: 1 variazione di velocità  ( U t   t  U t  t ) 2  t  t  1  ( U t  t  2 U t  U t  t ) variazione di accelerazione U t 2 Dove U è assunta come funzione degli spostamenti agli istanti t  t , t, t  t . L’istante in cui si va a effettuare l’equilibrio dinamico è t e quindi di ha che: t  U

 t  [C]  U  t  [K ]  U t  R t [M s ]  U s  t si giunge ad un’equazione di soli spostamenti, che  t e di U Sostituendo le espressioni di U espressa in maniera compatta ha la seguente forma:

ˆ  U t  t  Rˆ t M Dove:

ˆ a Ma C M 0 1 Rˆ t  R t  (K  a 2 M) U t  (a 0 M  a 1C) U t  t All’inizio della procedura risolutiva è necessario una inizializzazione che consiste nel calcolo degli spostamenti all’istante  t , che è pari a:  0 a U  0 U  t  U 0   t  U 3

con a 0 , a 1 , a 2 e a 3 costanti di integrazione che valgono: 1 ; t 2 1 a1  ; 2t 2 a0 

a 2  2a 0

61

a3 

1 a2

Di seguito si riporta la parte di codice Matlab che esegue il calcolo con il metodo delle differenze centrali:

62

Metodo di Hubolt

In questo metodo le leggi di variazione di velocità e accelerazione sono definite all’istante t  t e hanno le seguenti espressioni: 1  (11  U t  t  18  U t  9  U t  t  2  U t  2 t ) 6  t 1  2  (2  U t  t  5  U t  4  U t  t  U t  2 t ) t

U t  t  U t  t

Quindi l’equazione di equilibrio dinamico viene scritta all’istante t  t , cioè:  t  t  [C]  U  t  t  [K ]  U t  t  R t  t [M s ]  U s

 t  t e di U  t  t si giunge ad un’equazione di soli spostamenti, che Sostituendo le espressioni di U espressa in maniera compatta ha la seguente forma:

ˆ  U t  t  Rˆ t  t K Dove

ˆ  K  a M  a C matrice di rigidezza equivalente K 0 1 Rˆ t  t  R t  t  M(a 2 U t  a 4 U t  t  a 6 U t  2 t )  C(a 3 U t  a 5 U t t  a 7 U t  2 t ) forza equivalente Con 11 ; 6t

a2 

5 ; t 2

a0 ; 2

a7 

a3 9

a0 

2 ; t 2

a1 

a5 

 a3 ; 2

a6 

a3 

3 ; t

a 4  2a 0 ;

Dall’espressione risolutiva si riscontra che lo spostamento U t  t è funzione degli spostamenti calcolati ai tre istanti precedenti; di conseguenza tali spostamenti vanno calcolati con una procedura

63 separata, e più precisamente nella nostra analisi tali spostamenti agli intervalli t, t  t e t  2  t vengono valutati mediante il Metodo delle differenze centrali e quindi si ha che

U 0 H  U 0 DC U t H  U t DC U 2 t H  U 2 t DC A questo punto è possibile valutare gli spostamenti a tutti gli intervalli successivi mediante l’espressione ricorsiva sopra esposta. Di seguito si riporta la parte di codice Matlab che esegue il calcolo con il metodo di Houbolt:

64 Metodo di Wilson  

Il metodo di integrazione nel dominio del tempo Wilson -  rappresenta un’estensione del “metodo di accelerazione lineare”. In tale metodo si assume lineare la variazione di accelerazione tra gli istanti t e t  t ; invece nel metodo di Wilson -  l’accelerazione varia linearmente tra gli istanti t e t    t . In base al valore assunto dal parametro  il metodo è condizionatamente oppure incondizionatamente stabile, in particolare se si assume   1,.37 il metodo diventa incondizionatamente stabile, pertanto nell’analisi si è definito   1,4 . La legge di variazione dell’accelerazione viene definita all’istante t   come:

 t   U  t    U  t  t  U  t  con 0      t . Integrando una volta tale espressione si può U   t ottenere l’espressione della velocità e integrando due volte si può ottenere quella degli spostamenti. Le leggi di variazione scelte sono quelle sopra esposte scritte all’istante t    t . Quindi si avrà che l’equazione di equilibrio dinamico va scritta all’istante t    t e avrà la seguente forma:

ˆ  U t t  Rˆ t t K Dove:

ˆ  Ka Ma C K 0 1

 t  2U  t )  C(a U t  2U  t a U  t ) Rˆ t t  R t   (R t  t  R t )  M(a 0 U t  a 2 U 1 3 Con le costanti di integrazione che assumono i seguenti valori: a0 

6 ; (t ) 2

a7 

t 2

a8 

a1 

3 ; t

a 2  2a 1 ;

a3 

t 2

;

a4 

a0



;

a5 

 a2



;

3 a6  1 ;



t 2 6

 0  0 , si possono determinare istante Una volta fissate le condizioni iniziali, che sono U 0  0 e U per istante i valori degli spostamenti nell’istante t    t , ma a noi interessa quello all’istante t  t e quindi dovremo ripassare a tale valore mediante le seguenti relazioni





 t  t  a  U t t  U t  a  U  t a U  t U 5 6 4

65



 t  t  U  t a  U  t  t  U  t U 7





 t a  U  t  t  2  U  t U t  t  U t  t  U 8



Di seguito si riporta l’implementazione in Matlab di tale metodo:

66

Metodo di Newmark

Il metodo di Newmark è un’estensione del metodo Wilson -  , pertanto rappresenta anche un ampliamento del metodo di accelerazione lineare. In questo metodo le leggi di variazione sono le seguenti:

 t  t  U  t  [(1   )  U  t    U  t  t ]  t U  t  t  U  t  [( 1   )  U  t    U  t  t ]  t 2  t  t  U U 2 Affinché tale metodo risulti incondizionatamente stabile si deve avere che δ ≥ 0,5 e che α ≥ 0,25 · (0,5+δ). Ponendo  

1 1 e   si ottiene il “metodo di accelerazione media costante” 2 4

Da ciò si deduce che l’equilibrio dinamico viene effettuato all’istante t  t ; inserendo l’espressione delle leggi di variazione in quella dell’equilibrio dinamico si ottiene l’equazione risolutiva che in forma compatta può essere espressa come:

ˆ  Ka Ma C K 0 1

 t a U  t )  C(a U t  a U  t a U  t ) Rˆ t  t  R t  t  M(a 0 U t  a 2 U 3 1 4 5

67 Con le costanti di integrazione che assumono i seguenti valori:

 ; t

a0 

1 ;  t 2

a4 

t    1; a 5  (  2)  2 

a1 

a2 

1 ;  t

a3 

a 6  t (1   )

1  1; 2

a 7  t

Di seguito si riporta l’implementazione in Matlab di tale metodo:

68

3.3.2 Metodi di integrazione indiretta Metodo della sovrapposizione modale

I metodi di integrazione nel dominio del tempo appena visti determinano la soluzione del problema in forma discreta ovvero in un numero finito di istanti temporali. Tutte le metodologie ad ogni step di calcolo necessitano di un numero di operazioni da implementare pari a:

 mn Dove: - n = numero di gradi di libertà, cioè le dimensioni del problema - m = metà larghezza di banda delle matrici coinvolte -  2 Come detto questo è il numero di operazioni per ogni step di calcolo cioè per ogni istante temporale, quindi se S è il numero totale di istanti temporali in cui si è discretizzato il problema si avrà che il numero complessivo di operazioni è pari a: S   m  n

Per ridurre il numero di operazioni è necessario diminuire la larghezza della banda delle matrici coinvolte. A tal riguardo la numerazione dei nodi è importante ma, ovviamente, la riduzione della larghezza della banda è limitata dalla discretizzazione della struttura. Da ciò si deduce che sotto un certo numero di operazioni non è possibile scendere a meno di non effettuare una manipolazione delle equazioni risolventi. In particolare il metodo della sovrapposizione modale consente di stringere la larghezza della banda delle matrici operando una trasformazione di coordinate del tipo:

U( t )  [P]  {x ( t )} In pratica si introduce un nuovo sistema di coordinate x(t) definite coordinate generalizzate; nel caso in esame come matrice [P] si sceglie la matrice delle forme modali   poiché ci permette di ottenere la diagonalizzazione sia della matrice di massa sia della matrice di rigidezza, cioè:

 T  [M]     [I]  T  [K ]      2 

69 Quindi grazie a questa trasformazione l’equazione di equilibrio dinamico (sempre in assenza di smorzamento) diventa la seguente: [I]  x ( t )  [ 2 ]  x ( t )     F( t ) T

Avendo diagonalizzato tutte le matrici coinvolte si è ottenuto un sistema disaccoppiato che presenta il grande vantaggio di poter essere risolto mediante la soluzione di N equazioni distinta ad un grado di libertà ciascuna. La soluzione in termini di spostamenti generalizzati x(t) può essere calcolata mediante vari metodi, nel nostro caso si è scelto di utilizzare il metodo di Houbolt, con la stessa metodologia vista in precedenza. Una volta determinate le soluzioni in termini di coordinate generalizzate si può calcolare le soluzione in termini di coordinate del sistema attraverso l’equazione di passaggio di coordinate e cioè:

U( t )  [ ]  {x ( t )} Un altro grande vantaggio di questo metodo è che ci permette di ottenere una soluzione accurata (dal punto di vista ingegneristico) anche non considerando tutte le N equazioni, ma considerandone solo un loro numero limitato. Tale numero significativo di forme modali da considerare dipende da vari fattori quali: -tipo di struttura; - tipo di carico, inteso come contenuto in frequenza del carico e distribuzione del carico fra le varie forme modali (comunque il linea generale si può dire che: - per 2 analisi sismica n  10 , - per analisi di esplosioni n  N , - per analisi di vibrazioni n = tutti i modi 3 compresi fra  inf e  sup ). Nel nostro caso si studierà sia la soluzione ottenuta con una sovrapposizione di tutte le forme modali sia la soluzione ottenuta considerando la metà delle forme modali a disposizione e poi si confronteranno i risultati per valutarne l’errore. Di seguito si riporta la parte di codice che esegue il calcolo della sovrapposizione a modi completi:

70

71

Di seguito si riporta la parte di codice che esegue il calcolo della sovrapposizione a modi troncati:

72

73

3.4 Risposta in frequenza della struttura Ora pensiamo di effettuare un’analisi del sistema condotta per frequenze circolari comprese tra 0 e 1000 Hz con passo di 0,01 Hz che ha messo in evidenza il seguente andamento della risposta massima (su un tempo di 0,1 secondi). Da questo è possibile desumere quale siano le frequenze che inducono risonanza. Per fare ciò si impone che le frequenze delle due forzanti precedentemente definite siano uguali e che varino nell’intervallo prestabilito, per ognuna di queste frequenze, mediante il metodo della sovrapposizione modale (risolto con il metodo delle differenze centrali), si calcola prima la risposta e poi il massimo della stessa per poi plottare i vari massimi in funzione delle  . Di seguito si riporta la parte di codice che realizza quanto descritto:

74

75

3.5 Fenomeno dei battimenti Andiamo ad analizzare la risposta della struttura quando viene sollecitata con una forzante armonica avente pulsazione  alla prima pulsazione propria della struttura. Quindi questa volta le forzanti che andremo a definire saranno sempre due ma avranno una stessa pulsazione circa uguale alla prima pulsazione propria della struttura. Una volta definite le forzanti si procederà al calcolo degli spostamenti mediante il metodo della sovrapposizione modale risolto con il metodo di Houbolt. Di seguito si riporta la parte di codice che realizza quanto descritto:

76

77

3.6 Fenomeno della risonanza Studiamo il fenomeno della risonanza adottando due forzanti con pulsazione uguale alla prima pulsazione propria del sistema. Il fenomeno della risonanza è dovuto al fatto che quando la pulsazione della forzante è pari a una delle pulsazioni proprie del sistema il fattore di 1 amplificazione di quel modo  j  tende ad infinito e quindi gli spostamenti di quel (1  ( /  j ) 2 ) modo tendono ad avere un aumento indefinito nel tempo. Quindi ora si andrà a risolvere il sistema sottoposto a tali forzanti mediante il metodo della sovrapposizione modale risolto mediante il metodo di Houbolt. Di seguito si riporta la parte di codice che realizza quanto descritto:

78

79

3.7 Confronto grafico-numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di periodi propri e modi propri di vibrare Per prima cosa si riporta un confronto numerico fra i periodi propri calcolati con il Matlab e quelli calcolati con il SAP2000, più precisamente si riportano due tabelle contenenti rispettivamente i dati ottenuti dal SAP e dal Matlab e poi si riporta una terza tabelle contenente le differenze percentuali fra i due sistemi di calcolo:

TABLE: Modal Periods And Frequencies SAP2000 OutputCase StepType StepNum Period CircFreq Text Text Unitless Sec rad/sec MODAL Mode 1 0,0497 126,41 Mode 2 0,0259 242,65 MODAL MODAL Mode 3 0,0111 568,62 MODAL Mode 4 0,0094 670,18 MODAL Mode 5 0,0089 704,53 MODAL Mode 6 0,0072 878,77 MODAL Mode 7 0,0071 886,33 MODAL Mode 8 0,0053 1184,50 MODAL Mode 9 0,0036 1726,60 MODAL Mode 10 0,0034 1828,90

80 TABLE: Modal Periods And Frequencies Matlab OutputCase StepType StepNum Period CircFreq Text Text Unitless Sec rad/sec MODAL Mode 1 0,0497 126,41 MODAL Mode 2 0,0259 242,66 MODAL Mode 3 0,0110 568,63 MODAL Mode 4 0,0094 670,18 MODAL Mode 5 0,0089 704,53 MODAL Mode 6 0,0071 878,78 MODAL Mode 7 0,0071 886,33 MODAL Mode 8 0,0053 1184,51 MODAL Mode 9 0,0036 1726,60 MODAL Mode 10 0,0034 1828,90

Confronto periodi propri Matlab-SAP n°modo diff. percentuale / (%) 1 -0,000350 2 0,002206 3 0,002052 4 -0,003655 5 -0,002652 6 0,001192 7 0,000244 8 -0,009028 9 -0,001616 10 0,014443

Ora si andrà a riportare la parte di codice che esegue il plot delle forme modali della struttura mediante il calcolo della linea elastica, utilizzando di volta in volta come condizioni al contorno i valori dell’autovettore preso in considerazione (con la stessa procedura già vista per il plot della deformata statica):

81

82

83

84

85

Ora si vanno a riportare le prime 5 forme modali che si ottengono con questo plot e anche quelle ottenute con il SAP in modo da avere un confronto grafico fra i due codici:

86 MATLAB

SAP2000

87 MATLAB

SAP2000

88 MATLAB

SAP2000

89 MATLAB

SAP2000

90 MATLAB

SAP2000

91

3.8 Plot delle oscillazioni libere mediante Matlab Di seguito si riporta la parte di codice che realizza il plot delle oscillazione libere:

92

Di seguito si riporta il risultato di tale plot sia per il primo set di condizioni al contorno sia per il secondo:

93

3.9 Plot delle oscillazioni forzate calcolate con il metodo delle differenze centrali mediante Matlab Di seguito si riporta sia la parte di codice che realizza il plot delle oscillazioni forzate calcolate con il metodo delle differenze centrali che il grafico ottenuto da tale codice :

94

3.9 Plot delle oscillazioni forzate calcolate con il metodo di Houbolt mediante Matlab Di seguito si riporta sia la parte di codice che realizza il plot delle oscillazioni forzate calcolate con il metodo di Houbolt sia il grafico ottenuto da tale codice :

95

3.10 Confronto grafico-numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di oscillazioni forzate calcolate con il metodo di Newmark Di seguito si riporta sia la parte di codice che realizza il plot delle oscillazioni forzate calcolate con il metodo di Newmark:

Ora si va a riportare il confronto grafico fra la risposta calcolata con il Matlab (mediante il codice sopra riportato) e quella calcolata con il SAP: MATLAB

96 SAP2000

Per il confronto numerico fra i due codici di calcolo, siccome il numero di intervalli temporali necessari all’analisi è sufficientemente elevato (200001 intervalli), era impensabile eseguire un confronto su ogni singolo punto perciò si è eseguito un confronto su una serie di punti ad intervalli maggiori dell’intervallo di analisi, più precisamente il confronto è stato effettuato con un’ampiezza dell’intervallo di 0,1 s. Le differenze percentuali ottenuti in questi punti, per comodità di lettura dei risultati, sono stati riportati in un istogramma dove in ascissa abbiamo il tempo e in ordinata le differenze percentuali fra i due codici. Di seguito si riporta la parte di codice che calcola e diagramma tali differenze percentuali:

97

Il risultato che si ottiene da questo confronto numerico, come detto, è il seguente istogramma delle u  u Matlab differenze percentuali, calcolate come   SAP  100 (%) u SAP

98

3.11 Confronto grafico-numerico fra i risultati ottenuti con Matlab e quelli ottenuti con il codice di calcolo SAP2000 in termini di oscillazioni forzate calcolate con il metodo di Wilson   Di seguito si riporta sia la parte di codice che realizza il plot delle oscillazioni forzate calcolate con il metodo di Wilson   :

99 Ora si va a riportare il confronto grafico fra la risposta calcolata con il Matlab (mediante il codice sopra riportato) e quella calcolata con il SAP: MATLAB

SAP2000

100 Anche in questo caso per il confronto numerico fra i due codici si calcolano le differenze solo in alcuni punti (intervallati con un intervallo temporale di 0,1 s) e poi si riportano tali differenze in un istogramma. Di seguito si riporta la parte di codice che calcola e diagramma tali differenze percentuali:

101

Il risultato che si ottiene da questo confronto numerico, come detto, è il seguente istogramma delle u  u Matlab differenze percentuali, calcolate come   SAP  100 (%) u SAP

3.12 Plot delle oscillazioni forzate calcolate con il metodo di sovrapposizione modale a modi completi mediante Matlab Di seguito si riporta sia la parte di codice che realizza il plot delle oscillazioni forzate calcolate con il metodo della sovrapposizione modale a modi completi sia il grafico ottenuto da tale codice :

102

3.12 Plot delle oscillazioni forzate calcolate con il metodo di sovrapposizione modale a modi troncati mediante Matlab

103

3.13 Confronto grafico-numerico fra i vari metodi di integrazione nel dominio del tempo implementati in Matlab Ora si andranno a riportare una serie di confronti fra i metodi di integrazione sopra elencati di natura sia solamente grafica (sovrapposizione grafica di tutti i grafici ottenuti con i vari metodi) che grafica-numerico cioè un istogramma delle differenze percentuali (calcolato come visto sia nel metodo di Newmark che in quello di Wilson   ); tale differenze percentuali saranno calcolate per: - differenza fra metodo delle differenze centrali e metodo di Houbolt, - differenza fra il metodo delle differenze centrali e quello della sovrapposizione modale a modi completi, - differenza fra il metodo della sovrapposizione modale a modi completi e il metodo della sovrapposizione modale a modi troncati. Si inizia riportando il codice che realizza la sovrapposizione grafica dei grafici di tutti i metodi e di seguito il relativo grafico:

104

105

106 Ora si va a riportare la parte di codice che calcola e plotta le differenze percentuali fra il metodo delle differenze centrali e quello di Houbolt:

107

Il risultato che si ottiene è il seguente:

Ora si va a riportare la parte di codice che calcola e plotta le differenze percentuali fra il metodo delle differenze centrali e quello della sovrapposizione modale a modi completi:

108


109 Ora si va a riportare la parte di codice che calcola e plotta le differenze percentuali fra il metodo della sovrapposizione modale a modi completi e quello della sovrapposizione modale a modi troncati:

110


3.14 Plot della risposta in frequenza della struttura calcolata mediante Matlab Di seguito si riporta la parte di codice che esegue il plot della risposta in frequenza della struttura:

111


3.14 Plot del fenomeno dei battimenti della struttura calcolato mediante Matlab Di seguito si riporta la parte di codice che esegue il plot dei battimenti della struttura:

112


3.14 Plot del fenomeno della risonanza della struttura calcolato mediante Matlab Di seguito si riporta la parte di codice che esegue il plot dei battimenti della struttura:

113


114

4. Modellazione di una piastra quadrata tramite elementi shell L’elemento Shell nel SAP ha una formulazione a tre o quattro nodi che combina il comportamento separato a membrana (lastra) e quello a piastra flettente. La formulazione a quattro nodi permette di ottenere risultati più accurati pertanto quando possibile è da preferire. Il comportamento a membrana (o lastra) usa una formulazione isoparametrica che comprende le componenti di rigidezza traslazionali nel piano e una componente di rigidezza rotazionale nella direzione normale al piano dell’elemento. Il comportamento a piastra flettente comprende due componenti di rigidezza rotazionali della piastra, fuori dal piano, e una componente di rigidezza traslazionale nella direzione normale al piano dell’elemento. Per default viene usata una formulazione a piastra spessa (Mindlin/Reissner) che comprende gli effetti della deformazione di taglio trasversale. E’ possibile anche scegliere una formulazione a piastra sottile (Kirchoff) che trascuri la deformazione di taglio trasversale. Per conseguire la maggiore accuratezza possibile nei risultati la posizione dei nodi dovrebbe essere scelta in modo tale da rispettare le seguenti condizioni geometriche:  il valore di ciascun angolo interno fra due lati concorrenti deve essere minore di 180°; i risultati migliori si ottengono per angoli compresi fra 45° e 135°.  il rapporto di forma (aspect ratio) di un elemento non dovrebbe essere troppo elevato. Per il quadrilatero si considera il rapporto fra la distanza più lunga fra i punti mediani dei lati opposti e la più corta di queste distanze. I risultati migliori si ottengono quando i rapporti di forma sono vicini all’unità o almeno minori di quattro. Il rapporto di forma non dovrebbe superare dieci.  per il quadrilatero, i quattro nodi non devono essere necessariamente complanari. L’angolo fra le normali ai vertici fornisce una misura del grado di torsione. La normale ad un vertice è perpendicolare ai due lati che si incontrano in quell’angolo. Questo angolo non dovrebbe superare i 45°. Normalmente queste condizioni vengono rispettate con un’adeguata scelta della maglia. L’elemento Shell attiva sempre tutti i sei gradi di libertà relativi a ciascun nodo ad esso connesso. Quando l’elemento viene usato solo come membrana (lastra), l’utente deve assicurarsi che i gradi di libertà per la traslazione normale e per le rotazioni flettenti abbiano vincoli esterni. Quando l’elemento è usato solo come piastra, l’utente deve assicurarsi che i gradi di libertà per le traslazioni nel piano e per la rotazione intono alla normale abbiano vincoli esterni. Ciascun elemento Shell ha un proprio sistema di coordinate locale usato per definire le proprietà del materiale, i carichi e l’output. Gli assi di questo sistema locale sono indicati con i numeri 1,2 e 3. I primi due assi giacciono nel piano dell’elemento con orientamento specificato dall’utente; il terzo asse è normale.

115

4.1 Modellazione in SAP2000 Per analizzare lo stato di sollecitazione della struttura in esame è stato utilizzato il programma Sap2000, che permette di effettuare sia analisi statiche che dinamiche. Nel nostro lavoro ci limitiamo a studiare il comportamento della struttura in campo statico.

4.1.1 Definizione della struttura La seguente relazione ha per oggetto l’analisi delle sollecitazioni in una platea in calcestruzzo con le seguenti caratteristiche: -

Lunghezza

= 2m

-

Larghezza = 2 m

-

Spessore = 0.25 m

Il materiale scelto ha le seguenti caratteristiche: -

E = 33000000 KN/m2 ϒ = 25 KN/m3 ν = 0,3

116 Schema della struttura nel piano xy:

Schema della struttura 3d:

Il modello ad elementi finiti della platea è stato ricostruito mediante degli elementi bidimensionali a quattro nodi, gli shell.

117

Sono state adottate le seguenti meshature: -

4x4 8x8 16x16 32x32

Il modello è stato vincolato ai 4 lati adottando due schemi limite: a) cerniera b) incastro

4.1.2 Carichi: La piastra è soggetta a una forza concentrata applicata lungo l’asse z in corrispondenza del baricentro della platea pari a - 400 KN e ad un carico uniformemente distribuito sul piano xy della piastra pari a - 50 KN/m2. Nell’analisi non si è tenuto conto del peso proprio degli elementi (DEAD).

Fig.- Carichi applicati sulla struttura

118

4.1.3 Discretizzazione e scelta del tipo di elemento La struttura è stata modellata sia con elementi bidimensionali plate-thin che plate-thick in modo da tener conto anche delle deformazioni a taglio. Inizialmente la superficie totale è stata suddivisa in sottoaree 0,5x0,5 m ( 4x4 mesh) raffittite progressivamente in modo da avere una modellazione migliore. La mesh più fitta è costituita da elementi rettangolari 6,25x6,25 cm.

Mesh 4x4

Mesh 8x8

Mesh 16x16

Mesh 32x32

119

4.1.4 Risultati in termini di spostamento I risultati ottenuti con il Sap2000 sono stati poi confrontati con i valori ricavati analiticamente, di seguito riportati. In particolare è stato determinato l’abbassamento massimo della platea, che si ha in corrispondenza del baricentro della stessa, sia nel caso di carico concentrato che nel caso di carico uniformemente distribuito per entrambe le tipologie di vincolo. L’analisi è stata effettuata sia nel caso platethin che nel caso platethick. Il calcolo manuale degli spostamenti è stato eseguito utilizzando le equazioni riportate nel testo “Timoshenko e Woinowsky-Krieger 1959”. Introduciamo la notazione:

dove D è chiamata rigidezza flessionale della piastra. Vincolo cerniera, carico uniformemente distribuito (Tab.8 pag 120, con

, “Timoshenko e Woinowsky-Krieger 1959”)

Vincolo cerniera, carico concentrato (Tab.23 pag 143, con


Vincolo incastro, carico uniformemente distribuito (Tab.37 pag 206, con


120 Vincolo incastro, carico concentrato (Tab.37 pag 206, con


Di seguito sono riportate le tabelle con il confronto degli spostamenti. PLATE THICK ‐ CARICO CONCENTRATO

PLATE THICK ‐ CARICO UNIFORME

U3

U3

SAP2000 CALCOLO ERRORE

VINCOLO MESH SAP2000 CALCOLO ERRORE

VINCOLO MESH

4 x 4 0,000079 8 x 8 0,000079 Cerniera 0,000068 16 x 16 0,000078

13,92 13,92 12,82

4 x 4 0,000515 8 x 8 0,000509 Cerniera 0,000393 16 x 16 0,000516

23,69 22,79 23,84

32 x 32 0,000078

12,82

32 x 32 0,000530

25,85

4 x 4 0,000031 8 x 8 0,000028 Incastro 0,000021 16 x 16 0,000027

32,26 25,00 22,22

4 x 4 0,000283 8 x 8 0,000287 Incastro 0,000189 16 x 16 0,000297

33,22 34,15 36,36

32 x 32 0,000027

22,22

32 x 32 0,000312

39,42

PLATE THIN ‐ CARICO UNIFORME

PLATE THIN ‐ CARICO CONCENTRATO

U3

U3

VINCOLO MESH SAP2000 CALCOLO ERRORE

VINCOLO

MESH SAP2000 CALCOLO ERRORE

4 x 4 0,000069 8 x 8 0,000069 Cerniera 0,000068 16 x 16 0,000069

1,45 1,45 1,45

4 x 4 0,000430 8 x 8 0,000404 Cerniera 0,000393 16 x 16 0,000396

8,60 2,72 0,76

32 x 32 0,000069

1,45

32 x 32 0,000394

0,25

4 x 4 0,000025 8 x 8 0,000022 Incastro 0,000021 16 x 16 0,000022

16,00 4,55 4,55

4 x 4 0,000217 8 x 8 0,000200 Incastro 0,000189 16 x 16 0,000193

12,90 5,50 2,07

32 x 32 0,000021

0,00

32 x 32 0,000191

1,05

121

I risultati del Sap mostrano come aumentando il numero di mesh migliori il confronto con i risultati calcolati a mano. Il modello con 4 x 4 mesh e l’opzione thickplate, nel caso di carico concentrato, mostra una differenza rispetto ai risultati teorici del 23% per vincolo cerniera e del 33% per vincolo incastro. Lo stesso modello, con l’opzione thinplate, presenta differenze sensibilmente inferiori: 8% e 13%. Il modello con 32 x 32 mesh e l’opzione thickplate, nel caso di carico concentrato, mostra una differenza rispetto ai risultati teorici del 26% per vincolo cerniera e del 39% per vincolo incastro. Lo stesso modello, con l’opzione thinplate, presenta differenze sensibilmente inferiori: 0,25% e 1%.

4.1.5 Risultati in termini di tensioni Di seguito vengono riportati i diagrammi della tensione S11 nei vari casi estratti dal SAP2000, sulla destra vengono riportati le varie tensioni mediate, a sinistra le stesse sollecitazioni non mediate. Se il grado di raffittimento della mesh è adeguato non si dovrebbe assistere a brusche variazione di colore e dunque di sollecitazione tra un elemento e quelli adiacenti.

CARICO CONCENTRATO

4x4 Thin incastro non mediato

4x4 Thin incastro mediato

122





123



CARICO UNIFORMEMENTE DISTRIBUITO



124





125



CARICO CONCENTRATO

4x4 Thick incastro non mediato

4x4 Thick incastro mediato

126





127



CARICO UNIFORMEMENTE DISTRIBUITO



128





129



4.1.6 Conclusioni Partendo dall’analisi effettuata sugli spostamenti notiamo subito che nonostante il raffittimento della meshatura i risultati continuano a discostarsi in maniera sensibile da quelli teorici, ma questo è un risultato che ovviamente ci attendiamo in quanto le formule teoriche a cui si è fatto riferimento per il calcolo manuale degli spostamenti derivano da una trattazione che trascura gli effetti deformativi a taglio, ed è quindi logico che tali valori non siano pari a quelli calcolati con una formulazione thickplate in cui si tiene conto di tali effetti deformativi. Quindi per poter effettuare un controllo si deve far riferimento agli spostamenti calcolati con una formulazione thinplate; come si può vedere dalla tabella comparativa degli spostamenti tali valori così calcolati all’aumentare della meshatura tendono ad essere pressoché pari agli spostamenti calcolati con le relazioni teoriche con errori molto modesti, nell’ordine dello 0,25 – 1,5% a seconda della condizione di vincolo e di carico, nel caso di meshatura 32x32. Però si può anche notare che già per la meshatura 8x8 gli errori massimi riscontrati sono dell’ordine del 5% e quindi per una formulazione thinplate già questa meshatura porta ad una stima degli spostamenti accettabile dal punto di vista ingegneristico. Tale considerazione sulla sufficienza della meshatura non può essere estesa alla parte del calcolo delle tensioni, difatti se analizziamo le varie meshature thinplate vediamo che, benché le tensioni non subiscano bruschi salti fra un elemento e quello adiacente (anche con l’opzione di mediatura delle tensioni disattivata) già con la meschatura 4x4, l’andamento globale delle tensioni è molto differenti sia per forma che per valori fra le meshature più rade (4x4, 8x8) e quelle più fitte (16x16, 32x32) in entrambi i casi di carico considerati (cioè carico concentrato in mezzeria e carico distribuito sulla superficie della piastra). L’effetto è ancora più evidente nel caso thickplate in cui non solo la forma e i valori delle tensioni sono molto diversi fra le maschiature più rade (4x4, 8x8) e quelle più fitte (16x16, 32x32) in entrambi i casi di carico considerati, ma nelle meschature più rade si notano anche bruschi salti nella continuità delle tensioni fra un elemento e quello adiacente (ovviamente con l’opzione di mediatura delle tensioni disattivata). Quindi in conclusione si può dire che una meshatura 32x32 è adatta a modellare il problema oggetto della nostra analisi sia in termini di accuratezza nella valutazione degli spostamenti che in termini di calcolo tensionale sia per il modello thinplate che per quello thickplate. Detto tutto ciò è importante sottolineare che le valutazioni sugli spostamenti si fatte sono state riferite al modello thinplate in quanto si dispone di modelli teorici per tale formulazione, ed è quindi

130 l’unica con la quale possiamo fare un confronto; però nel calcolo effettivo degli spostamenti è meglio far riferimento alla formulazione thickplate in quanto (come riportato nello stesso manuale del SAP2000) è più accurata nel caso in cui lo spessore superi una dimensione compresa fra un decimo e un quinto della lunghezza (che è proprio il caso in esame visto che s = 0,25 m > l/10 = 2/10 = 0,20 m), tenendo però conto che tale formulazione è molto soggetta ad errori nel caso si abbiano elementi molto distorti, che però non è il caso in esame visto che si sono utilizzati elementi bidimensionali regolari. Quindi per modellare la piastra di fondazione in esame il modo migliore (sia in termini di calcolo degli spostamenti che delle tensioni) è scegliere una meshatura 32x32 utilizzando elementi shell con formulazione thickplate.

5 Allegati Si riporta di seguito la compilazione del programma Matlab in tutte le sue parti.

ESERCITAZIONE Con Confronto Matlab-SAP2000

Recommend Documents