Relazioni Esercitazioni Simulazione Numerica

Relazioni di simulazione numerica per l'ingegneria meccanica AA 2019-2020

Esercitazione 1

È richiesto di calcolare, mediante il metodo degli spostamenti, le reazioni vincolari e gli spostamenti nodali di un sistema reticolare, composto da dieci elementi, vincolato tramite cerniere in tre dei suoi nodi, conoscendo a priori i carichi applicati al sistema ed i loro punti di applicazione. Tale metodo prevede la creazione di una matrice di rigidezza relativa all'intero sistema, determinata a sua volta a partire dalle matrici di rigidezza dei singoli elementi, ottenute tramite la nota relazione in cui il parametro E rappresenta il modulo di Young del materiale costituente le travi, il parametro A la sezione trasversale (comune a tutte le travi), L la lunghezza dell'elemento considerato, e θ l'angolo di inclinazione di ogni elemento rispetto all'asse di riferimento orizzontale (asse x).

Nel nostro caso, porremo L pari alla lunghezza dell'elemento 10, dalla quale si potranno calcolare, tramite semplici relazioni trigonometriche, le lunghezze relative agli altri elementi. In particolare, avremo:

• E=207000 Mpa
• A=188.5
• L=2000 mm

Procediamo dunque con la numerazione globale degli elementi, dei nodi, e dei gradi di libertà relativi a questi ultimi: in cui sono stati progressivamente numerati inizialmente i nodi liberi e successivamente quelli vincolati. Tale numerazione risulta rilevante, come vedremo, nella determinazione dello sforzo computazionale (al variare della larghezza di banda della matrice di rigidezza del sistema) necessario alla risoluzione del sistema di equazioni risultante dal sistema stesso.

Dalla figura a sinistra, si ricavano le lunghezze dei singoli elementi:

• (θ) = = ;1 2 2
• = = ;3 4 2
• = = = (θ) ;5 6 9
• = = = 7 8 10

Il passaggio successivo è quello di determinare l'angolo di inclinazione θ, e quindi i nodi sinistro e destro di ogni singolo elemento (nodi 1 e 2 della numerazione locale). Si noti come, per la determinazione di ciascun angolo di inclinazione, si sia scelto l'angolo positivo minore ottenibile.

Determinazione matrici di rigidezza k di ogni elemento

k₁:
 0 0 0 0 1 0 1 0 1EA   2k  1 0 0 0 0L 2 4 3  0 1 0 1 45 6 1 2 5 1 0 1 0  60 0 0 0EA  k  2 1 0 1 0 1L 2 4

k₂:
  51 2 1 2 1 2 1 2   61 2 1 2 1 2 1 2EA  k   3 1 2 1 2 1 2 1 2L 2 3   1 2 1 2 1 2 1 2 43 4 7 8

k₃:
  31 2 1 2 1 2 1 2   41 2 1 2 1 2 1 2EA  k   4 1 2 1 2 1 2 1 2L 2 7   1 2 1 2 1 2 1 2 85 6 7 8

k₄:
  50 0 0 0 0 1 0 1 6EA  k  5 0 0 0 0L 2 2 7  0 1 0 1 89 10 5 6

k₅:
  91 0 1 0 0 0 0 0EA 10 k  6 1 0 1 0L 2 2 5  0 0 0 0 65 6 11 12

k₆:
  51 2 1 2 1 2 1 2   61 2 1 2 1 2 1 2EA  k   7 1 2 1 2 1 2 1 2L 11   1 2 1 2 1 2 1 2 129 10 7 8

k₇:
  91 2 1 2 1 2 1 2   101 2 1 2 1 2 1 2EA  k   8 1 2 1 2 1 2 1 2L 7   1 2 1 2 1 2 1 2 811 12 7 8

k₈:
 1 0 1 0 11 0 0 0 0EA 12 k  9 1 0 1 0L 2 2 7  0 0 0 0 87 8 13 14

k₉:
 1 2 1 2 1 2 1 2   81 2 1 2 1 2 1 2EA  k   10 1 2 1 2 1 2 1 2L 13   1 2 1 2 1 2 1 2 14

Assemblaggio della matrice di rigidezza K del sistema

Una volta calcolate le matrici di rigidezza dei vari elementi, ed aver assegnato ad ogni riga e colonna di queste ultime degli indici riferiti ai gradi di libertà corrispondenti ai nodi sx e dx dell'elemento stesso, per l'assemblaggio della matrice globale K non resterà che sommare, per ogni termine, tutti i coefficienti di rigidezza (sparsi nelle varie matrici k delle singole travi) riferite a quello stesso grado di libertà. La matrice K sarà quindi:

 2 2 0 0 0 2 2 0 0 0 0 0 0 0 0 0  0 2 2 0 2 2 0 0 0 0 0 0 0 0 0 0   0 0 2 0 1 1 1 1 0 0 0 0 0 0      0 2 2 0 2 2 2 1 1 1 1 0 0 0 0 0 0      2 2 0 1 1 3 2 3 / 2 1 / 2 0 0 2 0 1 / 2 1 / 2 0 0      0 0 1 1 1 / 2 2 3 / 2 0 2 0 0 1 / 2 1 / 2 0 0       EA 0 0 1 1 0 0 2 2 1 1 / 2 1 / 2 2 0 1 / 2 1 / 2 K        L 0 0 1 1 0 2 1 2 2 1 / 2 1 / 2 0 0 1 / 2 1 / 2    0 0 0 0 2 0 1 / 2 1 / 2 2 1 / 2 1 / 2 0 0 0 0   0 0 0 0 0 0 1 / 2 1 / 2 1 / 2 1 / 2 0 0 0 0     0 0 0 0 1 / 2 1 / 2 2 0 0 0 2 1 / 2 1 / 2 0 0  0 0 0 0 1 / 2 1 / 2 0 0 0 0 1 / 2 1 / 2 0 0   0 0 0 0 0 0 1 / 2 1 / 2 0 0 0 0 1 / 2 1 / 2   0 0 0 0 0 0 1 / 2 1 / 2 0 0 0 0 1 / 2 1 / 2

Per l'ottenimento della quale ci si è affidati al programma di calcolo MATLAB®. La relazione che lega il vettore delle forze esterne, agenti nei nodi, agli spostamenti dei nodi stessi risulta:

{F} = [K] {d}

Da cui, dividendo la matrice [K] ed i vettori {F} e {d} in modo da differenziare i carichi applicati ai nodi liberi “F_L” (noti) dalle reazioni sui nodi vincolati “F_V” (incognite), e gli spostamenti dei nodi liberi “d_L” (incogniti) da quelli vincolati “d_V” (noti e nulli), si può ottenere:

d_L = K_LL^-1 · (F_L − K_LV · d_V)

Le sottomatrici K_LL, K_LV (=K_VL) e K_VV possono essere estratte dalla K una volta specificate le condizioni al contorno, e saranno:

K_LL =  2 2 0 0 0 2 2 0 0 0 0 2 2 0 2 2 0 0 0 0    0 0 2 0 1 1 1 1      EA 0 2 2 0 2 2 2 1 1 1 1 k     LL L 2 2 0 1 1 3 2 3 2 1 2 0 0     0 0 1 1 1 2 2 3 2 0 2   0 0 1 1 0 0 2 2 1      0 0 1 1 0 2 1 2 2

K_LV =  0 0 0 0 0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0 0 0 0 0  k   LV 2 0 1 2 1 2 0 0L  0 0 1 2 1 2 0 0     1 2 1 2 2 0 1 2 1 2    1 2 1 2 0 0 1 2 1 2

K_VV =  2 1 2 1 2 0 0 0 0 1 2 1 2 0 0 0 0   EA 0 0 2 1 2 1 2 0 0  k VV  L 0 0 1 2 1 2 0 0 0 0 0 0 1 2 1 2   0 0 0 0 1 2 1 2

In cui ogni sottomatrice avrà dimensione (nxn), con n pari al numero di gradi di libertà (Gdl) dei nodi a cui si riferisce, ovvero 8 Gdl per i nodi liberi (sono 4) e 6 Gdl per quelli vincolati (sono 3), avendo ognuno di essi 2 Gdl sul piano.

Risoluzione del sistema

Considerando l'equazione vista in precedenza:

d_L = K_LL^-1 · (F_L − K_LV · d_V)

Si possono ricavare le seguenti relazioni: −1{d}_L = [K_LL] · {F}_L − [K_LV] · {d}_V, ovvero le soluzioni dei due vettori incogniti del problema, rappresentanti rispettivamente gli spostamenti dei nodi liberi ed i carichi applicati ai nodi vincolati (reazioni vincolari).

Essendo i vettori {d}_V e {F}_V noti, in base alle condizioni al contorno del problema, e pari a:

• {F}_L = [-30000 -15000 10000 0 0 0 0 0]
• {d}_V = [0 0 0 0 0 0]

Sostituendo, nella prima delle relazioni sopra citate, il vettore {F}_L si ottiene:

d_L = [K_LL] · {F}_L

Andando similmente a sostituire il vettore appena calcolato nella seconda relazione, si può ottenere il vettore delle reazioni ai nodi vincolati:

{F}_V = [K_VL] · {d}_L

Ottenendo così entrambi i vettori inizialmente incogniti. Si noti come non sia possibile procedere alla risoluzione delle reazioni vincolari senza aver prima calcolato gli spostamenti dei nodi liberi, essendo il vettore di questi ultimi calcolabile a prescindere dalla conoscenza di quello delle reazioni vincolari, e non altrettanto il contrario.

Dagli spostamenti di ogni nodo è possibile infine calcolare le deformazioni di ogni asta, tramite la relazione ε = (ΔL / L), e quindi lo stato tensionale agente in ognuna di queste, essendo σ = E * ε. Di seguito sono mostrate le soluzioni del problema, ottenute mediante uso di MATLAB®.

Analisi sulla larghezza di banda della matrice K

Considerando i vari termini che compongono la matrice di rigidezza K, ovvero K_ij, si avrà che, quando per ciascun coefficiente della matrice forza e spostamento saranno relativi a nodi appartenenti ad elementi diversi, esso risulterà nullo: ciò dipende dalla definizione stessa di questi coefficienti, che rappresentano i valori dello sforzo che si desta al nodo i a seguito di uno spostamento unitario del nodo j dello stesso elemento, essendo nulli tutti gli altri spostamenti nodali. La matrice K, applicata ad un sistema di più elementi, in cui determinati nodi saranno quindi non direttamente collegati fra loro, è quindi caratterizzata dalla presenza di molti termini nulli: essa risulta una matrice sparsa.

Una raffigurazione schematica di tale matrice può essere ottenuta tramite il programma di calcolo, in cui sono evidenziati tutti i termini diversi da zero (senza specificarne il valore): MATLAB®. Una stima della sparsità, in termini della sua ampiezza di semibanda, si può ottenere come:

β = 2(maximum difference between node numbers of the same element + 1)

, dove rappresenta la massima differenza tra la numerazione dei nodi appartenenti ad uno stesso elemento, mentre i gradi di libertà disponibili per ogni nodo. Come brevemente citato in precedenza, la larghezza di banda può variare con la scelta effettuata per la numerazione dei nodi: essa è una grandezza di interesse in quanto incide sul numero di operazioni che il calcolatore deve eseguire per risolvere il relativo sistema d'equazioni. Maggiore sarà questa, maggiore sarà lo sforzo computazionale del calcolatore. La condizione ottimale si avrà dunque nel caso di elementi il più possibile raccolti intorno alla diagonale, condizione raggiungibile minimizzando la differenza tra i numeri dei nodi connessi allo stesso elemento.

Di seguito sono mostrati alcuni esempi dei valori di larghezza di semibanda al variare della numerazione dei nodi, in cui sono raffigurati soltanto i termini diversi da zero delle matrici K, calcolate di volta in volta.

• Ampiezza di semibanda pari a 14, con massima di 6 tra i nodi 1 e 7, avendo sostituito il nodo 1 con il nodo 4
• Ampiezza di semibanda pari a 10, con massima di 4 tra i nodi 2 e 6, avendo sostituito il nodo 2 con il nodo 3
• Ampiezza di semibanda pari a 12, con massima di 5 tra i nodi 2 e 7, avendo sostituito il nodo 4 con il nodo 7

Esercitazione 2

È richiesto di analizzare la struttura reticolare trattata nella precedente esercitazione mediante l'utilizzo del software Ansys Apdl, calcolando nuovamente le deformazioni, gli spostamenti nodali, le reazioni vincolari, i carichi e le tensioni delle aste. Successivamente, i risultati dell'analisi andranno confrontati con una stessa struttura avente gli elementi trave al posto degli elementi asta. La procedura da seguire, all'interno di Ansys, consiste nel definire inizialmente il tipo di elemento, le caratteristiche del materiale, le costanti del problema (sezione degli elementi), la geometria del sistema, i carichi ed i vincoli, per ottenere infine una soluzione numerica formata da risultati in forma grafica (in termini di deformata della struttura) e da file di testo contenenti i carichi e le tensioni agenti sulla struttura. Di seguito verranno mostrati i passaggi necessari ad ottenere tali soluzioni di output attraverso l'interfaccia grafica del programma.

Analisi con “Link180”

Preprocessor

La struttura verrà inizialmente analizzata tramite un elemento monodimensionale a due nodi, ciascuno con tre gradi di libertà, capace di trasmettere solo sforzi assiali in trazione e compressione: l'elemento asta, indicato dal programma con la sigla “Link180”. Per poter scegliere tale elemento, è necessario selezionarlo dal corrispondente menu di Ansys, seguendo il percorso:

Preprocessor > Element Type > Add/Edit/Delete > Add

Che porta alla seguente finestra di dialogo:

Alla scelta dell'elemento segue l'impostazione riguardante il materiale da utilizzare:

Preprocessor > Material Props > Material Models

Nelle seguenti finestre di dialogo devono essere selezionate le specifiche riguardanti il comportamento del materiale (a sinistra) ed inserite le conseguenti caratteristiche elastiche (a destra): in cui EX rappresenta il modulo di Young, mentre PRXY il coefficiente di Poisson. Accedendo inoltre al menu:

Preprocessor > Section > Link > Add, 2188.5

Si può specificare la sezione delle aste, nel nostro caso pari a:

Una volta specificati quindi il tipo di elemento, il materiale e la sezione trasversale delle aste, si procede importando la geometria del sistema da una sequenza di comandi contenuta in un opportuno file di testo da copiare sul prompt del programma, contenente le coordinate dei keypoints, ovvero i nodi del sistema: fra i quali verranno inserite delle linee per rappresentare i vari elementi del sistema, tramite il comando:

Preprocessor > Modeling > Create > Lines > Lines > Straight Line

Utilizzando l'opzione “Pick” per selezionare i nodi da collegare. Costituita quindi la geometria, occorre creare l'opportuna mesh, che può essere inserita utilizzando il "wizard" per le mesh del programma, ricordando di specificare il numero di parti in cui debbano essere suddivise le aste (nel caso in esame pari ad 1, dunque asta intera), di selezionare i relativi elementi ed infine di avviare la procedura di meshing. Il comando è raggiungibile da:

Preprocessor > Meshing > MeshTool

Ed una volta selezionati gli elementi per cui eseguire il processo (nel nostro caso tutti) il programma porta alla seguente finestra:

Prima di passare alla fase “solution” per l'ottenimento dei risultati della prova, la correttezza del processo fin ora adoperato può essere verificata facendo ricorso alla funzione “plot-element”, che mostra tutti gli elementi del sistema (come linee congiungenti i nodi), oppure alla “element-list”, ovvero un resoconto sugli elementi ed i nodi inseriti nell'attuale analisi:

Prima di poter avviare la soluzione del sistema, per completare la simulazione è necessario definire vincoli e carichi (entrambi contenuti nel sotto-menu “Loads” di Ansys) sugli opportuni nodi. Entrambe le operazioni possono essere eseguite sia nella zona Preprocessor che in quella Solution, anche se, nel nostro caso, saranno avviate dal relativo sotto-menu nel Preprocessor:

Preprocessor > Loads > Define Loads > Apply > Structural > Displacement > On Keypoints

In cui, dopo aver selezionato i keypoints interessati, si accede alla finestra di dialogo nella quale andranno specificati i Gdl da vincolare (ad esempio uno spostamento lungo una direzione o tutti gli spostamenti):

In questo caso verranno vincolati con un incastro i nodi connessi alla parete (nodi 5, 6 e 7).