multi7.txt

Quello che faremo oggi in particolare è molto importante per i MEMS perchè è praticamente il principio di funzionamento spiega le basi il principio di funzionamento di moltissime strutture quindi la dinamica dobbiamo includere gli effetti della dinamica. Ne faremo due cose, la dinamica analisi modale che è l'argomento che direttamente poi impatta i membri, e poi comunque vedremo anche la dinamica nel tempo, la dinamica nel tempo, come si può veramente giocando e facendo piccole modifiche, quello che abbiamo, affrontare dei problemi di complessità crescente. Quindi cominciamo con questo. Nelle dispense avete già questo capitolo, capitolo 4, che è l'inclusione degli effetti di inerzia e quindi procediamo un po', avremo un paio di esercizi da fare, un po' di cose da fare insieme e un po' di teoria. Quello che io metto nelle slide, come avete ormai capito, sono sempre pezzi delle disperse, quindi non c'è assolutamente nulla di particolare. Allora, come dobbiamo arricchire il principio delle potenze virtuali che abbiamo considerato le ultime due lezioni in elasticità statica? Allora, sostanzialmente la novità è la presenza dei termini di inerzia, quindi quel primo termine lì contiene la potenza delle forze di inerzia, ro zero è la densità del solido, la densità iniziale del solido, u è lo spostamento, la derivata seconda quindi la sua accelerazione, si intende la derivata della velocità materiale quindi la velocità della particella proprio la derivata totale rispetto al tempo della velocità quindi la vera accelerazione magari lo vedremo più avanti ma se invece di scrivere le equazioni per un solido scrivessimo quelle per un fluido scrivessimo le equazioni di Stokes la struttura di questo principio rimane molto simile cosa cambia? cambia la legge costitutiva cioè non è vero che sigma è uguale ad A per epsilon cioè lo sforzo non dipende solo dalla deformazione, ma dipende anche da due termini, lo sforzo dipende da due termini, uno è la pressione e uno è il gradiente della velocità, quindi non dallo spostamento ma dipende dalla velocità. Cambia però la legge costituttiva. La seconda cosa è che in meccanica dei fluidi l'accelerazione non si esprime in questa maniera, in questa maniera lagrangiana, cioè dicendo che quella è la derivata totale della velocità, ma si usa l'approccio meridiano, quindi l'accelerazione ha un termine di derivata parziale, cioè posizione fissa, più termine convettivo, fa entrare la velocità, scavare il gradiente della velocità, sono i termini di convettione, ma è un modo diverso di scrivere l'equazione. Allora, quindi questo dipende dal tempo e questa equazione deve essere risolta per ogni valore del tempo in un intervallo che si vuole analizzare, ad esempio da 0 fino a tf, che è il tempo finale dell'analisi. E poi ci sono delle condizioni iniziali, adesso è diventato un problema iperbolico, abbiamo delle condizioni iniziali sullo spostamento, quindi per ogni x è il tempo 0, abbiamo un campo assegnato di spostamento e la velocità anche ecco qui ho sbagliato avrei dovuto mettere la derivata parziale ancora qui come mai è rimasto questo ah è rimasto però sulle slide non sulle... ah no, anche sulla dispensa lo devo segnare che devo cambiare per essere più coerente con la notazione. Il senso è chiarissimo, però la notazione sarebbe più corretta, una derivata, visto che ho usato la derivata parziale là, dovrebbe essere usata anche qui, numero 6. Ok, quindi vi spiego semplicemente qua, abbiamo la dipendenza da x e da t, quindi invece di usare il puntino è più corretto usare la derivata parziale rispetto al tempo quando poi invece avremo la separazione delle variabili delle quantità che dipendono solo da x oppure dal tempo allora metterò il puntino per indicare la derivata rispetto al tempo quindi era un po' un mix usavo questa annotazione non mi piaceva quindi derivata parziale di U rispetto al tempo, quindi velocità, uguale a un campo assegnato di velocità, quindi U0 e U0 punto sono due campi completamente indipendenti uno dall'altro che rappresentano lo spostamento iniziale e la velocità iniziale. Vedo delle altre piccole cose da cambiare. Ok, allora, questo quindi è il principio delle potenze virtuali che impone l'equilino per ogni valore del tempo nell'intervallo di analisi. Ok, decidiamo di affrontarlo con il metodo degli elementi finiti e quindi avremo un'interpolazione per gli spostamenti, interpolazione che non riporto più qui ma che abbiamo visto l'altra volta in cui c'è UH, cioè lo spostamento interpolato, è uguale alla sommatoria su tutti i nodi di questo U.J che è il generico coefficiente J della collezione delle incognite per le funzioni di forma globali. Abbiamo passato un po' di tempo settimana scorsa per chiarire come facciamo questa interpolazione globale. Le U e le P sono esattamente le stesse. Quando passiamo alla velocità, le funzioni di forma globali che utilizziamo sono esattamente le stesse di prima, cioè decidiamo di interpolare la velocità usando le stesse funzioni di forma globali e semplicemente chiamiamo in modo diverso i parametri delle interpolazioni, mettiamo il punto sopra per dire che sono le derivate temporari degli spostamenti nodali quindi quando io trovo uj in punto intendo la derivata temporale della jesima componente del campo di spostamento che è quindi funzione del tempo il campo è qui questo nella dispensa è giusto qui nelle slide invece è rimasto così l'interpolazione allora questo è un campo controlerò questo dipende da x da t si attribuisce la dipendenza dal tempo ai valori nodali e si attribuisce alla dipendenza spaziale solo alle funzioni di forma, quindi che composizione, l'espressione moltiplicativa del campo, come spesso si fa. Stessa cosa per l'accelerazione, visto che abbiamo bisogno di metterla all'accelerazione, la generazione sarà data dalle stesse funzioni di forma che moltiplicano le derivate seconde degli spostamenti nodali. Questo termine a noi non interessa all'interno di questa formulazione perché non abbiamo termini che dipendono dalla velocità, se ci fosse, lo vedremo più avanti, una dissipazione avremmo bisogno anche di questa interpolazione, in realtà quello di cui abbiamo bisogno là dentro è l'interpolazione dell'accelerazione. Per il resto il principio delle potenze virtuali è identico a quello che abbiamo visto l'altra volta, quindi la sola novità è la presenza dell'accelerazione. Quindi che cosa otterremo? Otterremo che i termini a noi noti questi daranno un contributo alla matrice di rigidezza ah sì perché questa ipotesi in alto sarà molto importante per semplificare alcune procedure tutti i problemi di dinamica che noi affrontiamo oggi hanno spostamenti imposti uguali a zero. Quindi vedremo ad esempio delle travi incastrate che vibrano e gli spostamenti sono imposti ma sono nulli. E così via. Quella è un'ipotesi che ci permette di andare un pochino più veloci nelle procedure di assemblaggio, non cambierebbe assolutamente lo schema generale della procedura. Per cui questa è l'interpolazione totale dello spostamento. Manca quello che avevamo invece in elasticità statica, la parte assegnata UD che noi ipotizziamo essere uguale a zero, solo per semplicità. Ok, quello che cambia è quello che si aggiunge a questo primo termine, questo primo termine che deve essere anche un fiscatizzato e ovviamente non c'è sto ghisce che generi un nuovo contributo di questa forma. Da dove viene? Rifacciamo i passaggi che avevamo già visto l'altra volta. Vse W che è il campo di test. W dipende solo dallo spazio. Cioè dato a un certo istante imponiamo la forma debole del equilibrio e scegliamo un capo test che dipende solamente dalle coordinate spaziali W quindi sarà la sommatoria su i che va da 1 fino ad n grande di questa W collezione dei valori natali di WI che moltimica la funzione di forma phi funzione di X questa è l'espressione che avevamo ottenuto alla fine della lezione scorsa la interpolazione globale della W che usa le stesse identiche funzioni di forma dello spostamento dei valori nodali diversi. Qui poi abbiamo due punti, x tempo, che è uguale invece alla sommatoria, uso una lettera diversa, quella che è anche lì indicata j, che dunque contiene le derivate nodali dello spostamento, i moltiplicata per le funzioni di forma phi che dipendono dallo spazio si intende che queste quantità o due punti sono dipendenti dal tempo se noi dobbiamo valutare l'integrale sul domino omega h di 0 derivata seconda di u udo anche seconda della U rispetto al tempo, derivata seconda della U rispetto al tempo, quindi la potenza sviluppata dalle forze di inerzia sul campo test di velocità di omega, beh, se immaginiamo di sostituire queste interpolazioni qua dentro, allora avremo una sommatoria doppia, sommatoria su i e j che vanno tutte e due da i ad n, di che cosa? Allora la sostituzione di w farà apparire questa w i che tengo a sinistra perché classicamente questa è un'abitudine, non vedete perché tenere i condizioni che riguardano la w a sinistra. Poi avremo il vero integrale che riguarda che cosa? Il prodotto scalare è chiaramente puntativo, quindi io adesso scrivo prima w e poi scriverò le forze di inerzia, avremo quindi la FI che dipende dallo spazio mettiamo la densità davanti, 0 che moltiplica scalarmente FIJ che dipende anche l'esolo dallo spazio integrata in globica poi chiudo la parentesi perché l'altra quantità che devo scrivere ovvero questo non dipende dallo spazio ma dipende solo dal tempo quindi abbiamo un due punti ok? questo oggetto qua questo oggetto qui dipende chiaramente dagli indici EJ è una funzione solo dello spazio e dei parametri del materiale, in questo caso il parametri del materiale da cui dipende è solo la densità, quindi per calcolare queste integrali avremo bisogno di che cosa? Avremo bisogno dell'informazione sulla densità e avremo bisogno dell'informazione sulla geometria. Quindi ci sarà una nuova routine elementare, così come c'era in elasticità per la rigidezza elementare, per le forze e così via, ci sarà una nuova routine che calcolerà la matrice di massa. Questa la chiameremo matrice di massa coefficiente ij. Ed è questa roba qua che c'è scritta qua. Assomiglia moltissimo a quello che abbiamo visto in elasticità. L'unica differenza è effettivamente la presenza di queste derivate seconde sui valori nodali. Questo è l'elemento nuovo. Cioè abbiamo per ora, si dice, effettuato una semidiscretizzazione in spazio, mentre per quanto riguarda il tempo la dipendenza è ancora continua non abbiamo fatto nessun tipo di discretizzazione questo è il concetto molto semplice come si calcola? come si costruisce la matrice M? Solita strategia. Si parte da un elemento e poi si fa l'assemblaggio. L'assemblaggio è identico, lo vedremo proprio nel codice, è inserito esattamente all'interno delle stesse righe in cui si fa l'assemblaggio della matrice di rigidezza. L'assemblaggio della matrice di massa è identico. Però vediamo un po' una piccola differenza a livello del calcolo della matrice elementare. Il primo passaggio è, dato un elemento, calcoliamo il contributo dell'elemento alla matrice di massa. E poi facciamo l'assemblaggio. L'assemblaggio non lo commenteremo più perché è veramente inutile, è sempre la stessa cosa. Però quando io voglio calcolare sull'elemento generico il numero E il contributo alla potenza delle forze di inerzia, come faccio? Ripassiamo perché magari questo è utile, poi arriviamo a capire qual è la piccola differenza che lo contraddistingue dalla matrice di resistenza. Sul singolo elemento noi abbiamo sempre rappresentato U e W in una maniera leggermente diversa abbiamo sempre rappresentato U in questa maniera attraverso una lista abbiamo trasformato il vettore in una lista che contiene le due componenti U1 e U2 magari ci mettiamo l'H a pedici per dire che è il anche qui dovrebbe essere più corretto mettere l'h a pedici per dire che è il, anche qui dovrebbe essere più corretto mettere l'h a pedici per dire che è il campo interpolato agli elementi finiti. E l'abbiamo espresso come prodotto di una matrice n, funzione delle coordinate parametriche, che moltifica il campo w, il tentore w, la lista w, anzi. Vi ricordate? Allora,? questa era una matrice questa è una lista 2 per 1 in una colonna che contiene due coefficienti questa è una matrice in nozioni 2 per 12 le due linee e 12 perché deve moltiplicare questa lista che contiene i valori nodali che è di 12 elementi cioè le component componenti sul primo nodo, le due componenti sul primo nodo seguiti dalle due componenti sul secondo e così via. Questo abbiamo detto non lo giustificheremo mai più e lo useremo sempre. Ora, evidentemente, quando abbiamo visto per calcolare il nuovo contributo alla matrice di inerza, questo verrà utilizzato anche per la derivata seconda dello spostamento. Quindi quando dovrò fare la derivata seconda di OH rispetto a T, creerò una lista con le due componenti dell'accelerazione e questa si esprime esattamente così. Cioè le funzioni di forma sono le stesse e la derivata rispetto al tempo la attribuiamo tutta alla lista dei valori nodali. Vi stupisce qualcosa? Non fatevi problemi a chiedere se dico qualcosa che mi sembra strano o da giustificare meglio. Quindi, l'interpolazione a livello dell'elemento si fa esattamente come si faceva lo spostamento. E a maggior ragione ancora per la W. Allora, creiamo una lista W che ha le due componenti e questa sarà ancora data da n e alla W. Questo non ha bisogno di nessuna modifica è lo stesso che avevamo bene, se inseriamo questo all'interno del nostro integrale non facciamo tutti i passaggi ma questo risulterà essere uguale di E tras WI trasporto per l'integrale su omega, quindi avremo il trasposto della N, N trasporto per N, allora N trasporto è relativo a questo W che abbiamo trasporto per portare fuori WI. Quest'altro N invece è quello relativo all'accelerazione, qui ci mettiamo davanti la densità, poi dovremmo scrivere l'integrale in D omega, questa è la matrice di massa elementare per UE. Questa è la matrice che chiamiamo M, matrice di massa elementare. Quindi dobbiamo essere capaci di effettuare sul generico elemento questo integrale che è decisamente più semplice degli integrali che avevamo imparato a fare nella matrice di rigidezza perché non ci sono i gradienti, c'è solo lo spostamento come si fa? stessa tecnica che avevamo adottato per gli altri elementi uno, primo passaggio lo si mappa sull'elemento di riferimento. E, non lo scrivo qui perché abbiamo già l'espressione qua, cioè lo si mappa dall'elemento fisico sull'elemento di riferimento delta E, cioè il triangolo simplesso standard bidimensionale. E quindi avremo l'apparizione dello Iacobiano nell'elemento di superficie. Quindi quello che noi dobbiamo calcolare è questo oggetto. Come si fa? Tecnica numerica. Vi ricordate come avevamo detto che si faceva, si utilizzavano le regole di integrazione numerica. Ecco, qui arriva la piccola differenza, quindi facciamo attenzione. Quando noi dobbiamo valutare l'integrale solo di un elemento triangolare, il siplesso standard in 2D, cioè questo elemento nel piano parametrico a 1 e a 2, questo è 1 e questo è 1, questo si chiama siesto, è il nostro delta per un movimento triangolare, sempre lui, non cambia mai, allora ci sono delle regole, ci sono delle regole numeriche, le regole di Gauss-Hammer, e facciamo un piccolo ricollo, avevamo tante regole che vedono un numero di punti di Gauss crescente, e ogni regola è associ, non dico niente di nuovo è la stessa identica cosa che ho letto l'altra volta è solo un ripastino questo P è la cosa che più mi interessa dice che qual è l'ordine del polinomio in A1 e A2 che può essere integrato in maniera esatta. Quindi la prima regola che prevede un solo punto di Gauss, qui penso che ci sia, per questo qui è certamente un errore, questo 3 qui lo devo correggere 1, perché questo indica la moltiplicità, ma non ci interessa, questo qui, questa è la colonna, possiamo anche non guardarlo. La prima regola che prevede un solo punto di Gauss che sta nel baricentro dell'elemento con peso un mezzo può integrare un polinomio lineare in maniera esatta. Se invece passiamo alla regola tre punti che abbiamo poi inserito nel calcolo della matrice di rigidezza elementare, riusciamo a calcolare un polinomio di ordine 2, cioè vuol dire un polinomio che contiene a1 quadro, a1 per A2, A2 quadro e poi tutti i termini lineari. Quindi un polinomio di ordine 2 nei due parametri di superficie A1 a 2. Molto spesso non è sufficiente, allora dobbiamo salire, dobbiamo salire, utilizzare le regole ad esempio con 6 punti che garantiscono un'accuratezza del quarto ordine, cioè qualunque polinomio di ordine 4 può essere integrato accuratamente. E se vi ricordate, vi avevamo anche suggerito la regola che tipicamente negli elementi finiti viene utilizzata per scegliere tra queste integrazioni, perché molto spesso qui da integrare non abbiamo polinomi, ok? Non abbiamo polinomi, però per scegliere la regola, perché ci sono dei denominatori che quindi complicano, abbiamo rapporti fra polinomi, per scegliere la regola di integrazione si fa l'ipotesi che lo jacobiano sia una costante e ci si chiede, nel caso, questa è la procedura standard, nel caso in cui lo jacobiano sia una costante, cioè che l'elemento si possa deformare ma non distorcere, possa ingrandirsi ma non distorcere, qual è la regola che ci permette di integrare in maniera corretta il nostro integrando. Bene, proviamo a vedere. Supponiamo che ρ0 sia una costante, cioè che la densità sia costante all'interno del nere. Se j è una costante, allora dipende tutto da questo prodotto n. Ed n che cosa sono? n sono le funzioni di forma, vi ricordate come erano le funzioni di forma quadratiche è una funzione di forma che ognuna delle n quindi è una funzione quadratica in 1 e 2 se facciamo dunque il prodotto abbiamo un polinomio all'ordine 4 in generale dunque se vogliamo integrare in maniera corretta questa vanticea di massa dovremmo utilizzare una regola con 6 punti rigaio ed è effettivamente la scelta che viene fatta prima di andare a vedere tutto per non lasciare così vi ho aggiornato la distribuzione dei file e adesso avete anche il file no, questo non volevo aprire il dynamics perché questo è il secondo ma modal... va bene, comunque adesso non cambia grazie ci metto un istante come sempre ecco qua dopo vediamo la funzione model ma quello che ci interessa guardare adesso è la routine che calcola la matrice di massa elementare che siuss quindi se andate a confrontare questi sono gli stessi coefficienti che abbiamo visto nella slide e definiamo quindi sei punti di gauss faremo un loop dunque su sei punti di gauss e non più su tre come era invece il caso della matrice di rigidezza se guardiamo la matrice di rigidezza invece vi ricordate che aveva solamente 3 punti di Gauss perché in quel caso era sufficiente perché avevamo un polinomio solamente di ordine 2 per il resto che cosa cambia? non cambia quasi nulla però guardiamola si comincia a inizializzare una matrice di 0 e 12 per 12, perché M ha le dimensioni di 12 per 12, poi si fa un loop su tutti i punti di Gauss, prende la coordinata parametrica corrente A1 e A2, continua esattamente come per la matrice di residenza, calcolando la derivata delle funzioni di forma, la gradiente f. A chi servono, visto che non dobbiamo calcolare i gradienti delle funzioni? Servono per calcolare lo jacobiano. Il calcolo di d e di f, in questo caso, per la matrice di massa, serve solo e esclusivamente per arrivare al calcolo dello jaciano per effettuare l'integrale numero. Poi però vi ricordo che dobbiamo integrare questo oggetto qua, n trasposto n, dove n è quella matrice, la matrice che contiene due righe e dodici colonne che moltiplica i valori nodali. Vi ricordate com'era la formula della matricenne? Noi possiamo vedere andando a vedere le slide dell'altra volta. Ecco, questa è la matricenne, quella che contiene le funzioni di forma, le sei funzioni di forma alternate. Nella prima riga ci sono dei zero in posizione pari, nella seconda riga ci sono dei zero in posizione dispari. Perché? Perché? Semplicemente per la causa dell'ordine che abbiamo introdotto in UE, che contiene la prima componente e la seconda componente per il primo nodo e poi così via. Quindi se vogliamo la prima componente dello spostamento dipenderà solo dalle prime componenti su tutti i nodi, quindi tutte le posizioni dispari nella lista, per cui ci sono zero nelle posizioni dispari. E viceversa, nella seconda ci sono 0 in tutte le posizioni dispari quindi la n è quella funzione lì la n però non l'abbiamo mai implementata non l'abbiamo trovata mai in un codice e adesso effettivamente invece nella matrice di massa troviamo proprio questa questa è la matrice n che abbiamo appena visto prima definisco la lista delle funzioni di forma nl n lista n lista, che contiene le sei funzioni di forma, quelle che abbiamo calcolato al primo giorno, e poi compila la matrice n, 2 per 12, che poi ci permette di fare l'integrale numerico, ρ per n trasposto, per n per lo jacobiano, per il peso di Gauss, che è esattamente quello che, questo lo giudano, non mi serve più, che è esattamente quello che devo fare per effettuare questo integrale numericamente. Niente di nuovo, quindi. L'unica differenza rispetto alla matrice di rigidezza è veramente il fatto che abbiamo dovuto utilizzare una regola di Gauss più ricca, perché il polinomio che dobbiamo integrare è di ordine più elevato. Ok? Bene. Alcune osservazioni sulla matrice di massa, perché sarà fondamentale, perché se la matrice di massa non avesse questa proprietà, tutto quello che diremo dopo, ma anche il funzionamento di tutte le nostre strutture non sarebbe lo stesso. Allora, prima ci ricordiamo da dove viene la matrice di massa. Noi l'abbiamo visto, l'abbiamo discusso, viene da questa discretizzazione. Bene, la matrice di massa appare però anche se calcoliamo un'altra quantità, ovvero l'energia cinetica. Perché è importante l'energia cinetica? Perché innanzitutto è sempre positiva. L'unico caso in cui l'energia cinetica può essere zero è quando la velocità è identicamente nulla in tutti i punti. Quindi, escludiamo questo caso, perché supponiamo che il corpo sia in movimento, e allora ci aspettiamo che questa, che io ho chiamato KH, cioè l'energia cinetica discretizzata agli elementi finti, sia una quantità positiva, sempre, qualunque, questa è la base. Qual è la definizione dell'energia cinetica? Scuso che ho mischiato un po' S, un po' omega, non ho usato S nel codice ma omega nella teoria generale perché era fatta per il tridimensionale mentre S nel codice era usato perché il codice è bidimensionale però meglio fare un po' meno confusione a questo livello. Comunque, non è possibile. S è un integrale di volume o di superficie di ron zero per la velocità al quadrato ora, se vogliamo integrare questa cosa è meglio vedere questo modulo al quadrato come accelerazione, prodotto scalare per l'accelerazione scusate, velocità, prodotto scalare per la velocità cioè il modulo di un vettore al quadrato è il prodotto scalare del vettore con se stesso quindi facciamo velocità scalare velocità quindi noi avremo da valutare questo passaggio non l'ho messo nella slide però lo possiamo commentare molto facilmente qui dobbiamo integrare sul volume questa quantità, rho zero per la derivata di u h rispetto al tempo, con toscalar la derivata di u h rispetto al tempo. Quell'oggetto di S è d'accordo? Di stessa cosa, no? Di omega z. Questo è quello che noi dobbiamo calcolare. Quello che noi abbiamo calcolato prima invece è il termine che c'è là sopra, quindi l'integrale di rho omega su rho omega a di ρ omega su ρA della derivata seconda di UH rispetto al tempo scala W in di omega e avevamo fatto vedere che questo era effettivamente W trasposto per M per U due punti, no, l'abbiamo commentato pochi minuti fa come abbiamo fatto? Abbiamo inserito all'interno di questo integrale la discretizzazione agli elementi finiti con le funzioni di forma. Se facciamo la stessa cosa per la velocità, che ovviamente sarà, torniamo a, dove sarà il prima, la velocità è espressa nella stessa identica maniera, qua, no? Come sommatoria delle stesse punte di forma globali per i valori nodali derivati, se inseriamo la stessa identica discretizzazione qua dentro, bene, mi permetto di saltare questo passaggio, è ovvio che otterrò U.MU. Cioè, faccio la stessa procedura, quello che cambia è solo l'interpretazione di questi due campi, è ovvio che la distedizzazione porterà all'apparizione della stessa identica matrice di massa. Siamo d'accordo? L'altra differenza, la piccola differenza, è che per definire l'energia cinetica abbiamo l'un mezzo davanti. Ma questo non cambia la sostanza. L'energia cinetica è un mezzo di quella quantità. Ora, se io voglio che l'energia cinetica sia una quantità sempre positiva, è ovvio che sto imponendo delle condizioni sulla quantità di massa, m, perché se qualunque punto, io dico che dire che l'energia cinetica è sempre positiva vuol dire che qualunque punto questa quantità è estremamente positiva ed è zero solamente se tutti i valori nodali di questa roba sono zero. Beh, questa è la definizione stessa di matrice definita positiva. m è quindi una matrice che non solo è simmetrica, perché la definizione di M e J fa vedere che scambiando I con J otteniamo lo stesso coefficiente, quindi la matrice M è simmetrica. Ma cosa ancora più importante è una matrice definita positiva. E non ha gli stessi problemi della K, ricordate la K era una matrice che era definita positiva se aveva bloccato tutti i moti rigidi. Infatti vi ricordavate che aggiungevamo i vincoli puntuali per impedire che ci fossero le traslazioni? Ora, se siamo in dinamica, un corpo in generale si può anche muovere, lo vedremo. Voi, nell'esempio, l'avete visto nella prima slide, simuleremo la caduta di un corpo in generale si può anche muovere lo vedremo, voi nell'esempio l'abbiamo visto nella prima slide simuleremo la caduta di un corpo se un corpo cade non ha vincoli ci sono i moti rigidi possibili quindi la K in dinamica in generale a meno che non ci siano delle condizioni particolari è solo semi-definita positiva cioè se ci sono dei mot rigidi possibili, questi motri rigidi danno energia elastica zero, per cui la K è solo semi-definita positiva, cioè non è fatta così bene come le matrici di massa. Però va bene, noi circondiamo quindi a partire da questo momento che quando terminiamo la nostra discretizzazione arriviamo ad un'espressione di questo genere. Da dove arrivano questi termini? Questo dalla discretizzazione delle potenze virtuali delle forze di inerzia, questo dalla potenza virtuale degli sforzi interni sull'arte di sigma e questo colleziona tutto il contributo delle forze esterne. Questo valeva anche in elasticità, cioè forze di volume, forze di superficie, vanno a finire tutte quelle, con le condizioni iniziali sui valori dei spostamenti nodali, quindi il valore iniziale di U è questo campo assegnato, il valore iniziale della velocità è quest'altro campo assegnato. Proprietà fondamentali, la matrice di massa è simmetrica e definita positiva, la matrice K invece è simmetrica ma solo semidefinita positiva. Ok? Quindi il primo step è terminato, il primo step consiste nella semidiscretizzazione in spazio. Lo ripeto perché qui la dipendenza di U dal tempo è ancora continua. Questo è un enorme sistema di equazioni differenziali ordinarie del secondo ordine, perché appare la derivata seconda. E dobbiamo cercare di risolverlo. Questo è il secondo step, la seconda fase di ogni problema di dinamica. Una volta fatta la discretizzazione nel senso degli elementi finiti, a partire da questo momento non c'è nessuna differenza rispetto alle tecniche che probabilmente avete già visto in qualche corso per integrare qualche equazione di differenza ordinaria. Magari ne avete viste una, due, un sistemino piccolino, quindi tipicamente invece abbiamo migliori di equazioni accopiate tra di noi. Però il principio è lo stesso. Allora, a partire da questo è dunque il nostro problema di riferimento. Ci dimentichiamo da quest'istante in poi da dove arriva questo. Siamo benissimo, arriva dagli elementi finiti, però per quello che faremo dopo non ci interessa più, tranne in qualche esercizio, e adesso vediamo quali sono le tecniche di soluzione di queste equazioni, e vediamo le due cose di prima, l'analisi modale, che cos'è, e l'integrazione del tempo, che cos'è e come funziona, con degli esempi pratici per tutte e due le cose. Abbiamo anche due codici separati per la soluzione di questi problemi, il primo si chiama modal e l'altro si chiama dynamics. Analisi modale. Scusate, una cosa che non commenterò, ma ci ritorneremo magari le settimane prossime. Qui manca un termine di dissipazione. In genere, la dissipazione è difficilissima da modellare, si aggiunge in maniera artificiosa un termine che sia proporzionale alle velocità, si dice dissipazione alla Rayline, che contiene un termine proporzionale alla massa, un termine proporzionale alla rigidezza moltiplicata per la velocità. Però, non è quello anticipo, ma lo vedremo più. Lo vedremo anche nel caso scalare. Analisi modale. Questo è il punto fondamentale, le strutture. La cosa che differenzia enormemente le strutture e i MEMS da altre cose come i fluidi. Nei fluidi non si fa l'analisi modale perché le cose che troveremo in cosiddetti modi non hanno la stessa importanza drammatica che hanno nelle strutture. È un'esperienza abbastanza comune, che se si prendono delle strutture magari un po' flessibili, se le si sollecita ad una certa frequenza queste rispondono staticamente, ma se le si eccita con una frequenza maggiore a un certo punto queste cominciano a vibrare tantissimo. Ci sono delle frequenravano tutte, poi se aumentavate a 40 la vibrazione passava, ma gli aerei stessi quando decollano c'è sempre spesso un momento in cui sembra che tutto vibri, poi quando i giri motori aumentano, la velocità di riduruzione del motore aumenta e tutto passaioè, ci sono delle frequenze di sollecitazione che eccitano tremendamente le strutture. Le mandano, si dice, in risonanza. Le mandano in risonanza. È un'esperienza continua, magari non ci avete mai fatto caso, ma adesso tra un po' ci farete caso. I MEMS, nella stragrande maggioranza dei casi, funzionano a risonanza. C'è solo la categoria degli accelerometri capacitivi che non funziona a risonanza, ma funziona quasi staticamente, cioè sente l'accelerazione, si spostano in maniera quasi statica, si misura la variazione di capacità e si dice di quanto si sono spostati. Ma i giroscopi vibrano sempre. I clock, i risonatori che servono per il clock 5G sono sempre in risonanza. I microspecchi sono sempre in risonanza. Tutte queste strutture funzionano secondo un modo proprio di vibrare. Che cosa diavolo è? Allora, è facilissimo. Cioè, in realtà no, ci vuole un po' però alla fine vedrete che calcolarli sapendo fare le due cose che abbiamo imparato finora è drammaticamente facile poi il difficile sono gli algoritmi che a partire da MK calcolano gli autovalori e gli autovettori che perché sono raffinatissimi, non ce ne occupiamo, tutti i grandi codici li hanno implementati, sono delle librerie che lo fanno. Proviamo questa equazione, questa qua, che è quella di prima, e mettiamo f uguale a 0. Mettiamo f uguale a 0, cioè noi vogliamo vedere se esistono delle soluzioni di quell'equazione senza eccitazione esterna e come standard in tantissimi campi della fisica matematica cerchiamo soluzioni quindi avendo portato questo tale a zero cerchiamo soluzioni di questo genere cioè cerchiamo la U come prodotto di una quantità che chiameremo modo di vibrare questa è la collezione dei valori nodali di quello che chiameremo modo di vibrare che moltiplica questo esponenziale compresso l'esponenziale di omega t dove l'unità immaginaria omega la chiameremo frequenza propria ci sono mille nomi sono di alto valore e alto vettore ma in meccanica li si chiamano modi di vibrare, frequenze proprie ma sono alto vettore e alto valore d'accordo? allora quindi immaginiamo di inserire questo x è una quantità che non dipende dal tempo sono dei valori non alifissi l'immedizzare il tempo è tutto sparato all'interno dell'esponenziale l'esponenziale complesso lo buttate all'interno della vostra equazione e cosa ottenete? è che m u due punti immagino che non sia la prima volta che trovate una cosa di tutto genere, ma diventa meno omega quadro per e alla i omega t, no? Che moltiplica x. Invece se lo buttate all'interno di k, questo rimane semplicemente e alla i omega t per x. Allora, e alla i omega t appare da tutte le due parti, la si butta via e quindi otteniamo meno omega quadro m più k x uguale a 0. Questo è il sistema che dobbiamo risolvere. Quindi, le x, le coppie di x e omega che risolvono questa equazione sono i modi di vibrare x e omega che risolvono quella sono i modi di vibrare e le frequenze proprie e le frequenze proprie della struttura giroscopio un giroscopio mens a che frequenza? 20.000 Hz Hz. Viene costruito, viene disegnato per avere un modo intorno ai 20.000 Hz. La sua frequenza propria è quindi 20.000 Hz. In generale hanno più modi alla stessa frequenza, vedremo perché. Avendo già avuto voi tutti più o meno un corso MEMS con Corigliano, queste cose dovreste averle in testa. Adesso si vede come calcolare effettivamente. Allora, quindi otteniamo questo sistema qua. Questo problema, anche all'interno di MATlam, si chiama problema generalizzato agli autovalori. Perché noi siamo abituati a risolvere fin dal primo anno i problemi agli autovalori, ma i problemi agli autovalori in generale sono messi così. A per x, almeno a me mi insegnavano così, no? L a lambda x. Questo è un problema standard dagli autovalori. Invece noi abbiamo un problema leggermente diverso perché appaiono matrici, non solo abbiamo l'equivalente della A che è diventata la K, ma invece di avere lambda per x abbiamo, questo potrebbe benissimo giocare il ruolo di lambda, però abbiamo anche la matrice di massa che interviene. Quindi non è un problema standard agli autovalori, però lo si può ricondurre molto facilmente a un problema standard, grazie alle fantastiche proprietà della matrice di massa allora come si fa? allora m siccome una matrice è definita positiva la possiamo esprimere come m alla 1,5 per m alla 1,5 come si fa alla radice di una matrice? l'avete mai trovata? si diagonalizza la matrice si esprime nella base cioè come si calcano tutti gli autovalori della m che sono positivi e sarà la matrice con gli autovalori radice di quella iniziale comunque per definire le potenze e le radici delle matrici si passa sempre nella versione diagonalizzata delle matrici. Si diagonalizza, poi si prende la radice della diagonale e poi si ritorna indietro. Questo si può fare perché M è una matrice definita positiva. Tutte le matrici simmetriche definite positive si possono esprimere così. Allora, come si fa? Allora, inseriamo, rappresentiamo m come m alla un mezzo per m alla un mezzo, e poi sostanzialmente introduciamo e definiamo m alla un mezzo m alla un mezzo che è x introduciamo un nuovo vettore che teniamo y ok se introduciamo questo nuovo vettore che è definito così allora possiamo ved, avremo, facciamolo ma, avremo, contemporaneamente premoltiplichiamo a sinistra per m alla meno un mezzo tutte le equazioni, quindi avremo m alla meno un mezzo, che moltiplica k, poi invece di scrivere x, scrivo m alla meno mezzo per y. Ok, questo è come ho modificato questa qua. Poi ho meno omega quadro, che moltiplica m alla mezzo per e alla mezzo, quindi se premoltiplico per m alla meno mezzo, il primo m allo mezzo sparisce, rimane semplicemente m allo mezzo per x che è x. Quindi con questo trucchettino, quindi di definire il nuovo vettore y in questa maniera e con il fatto che ho premoltiplicato tutto per mello mezzo, ho trasformato il mio problema agli autovettori generalizzato in un problema agli autovallori standard. Perché questa è diventata la matice A del problema che ho scritto prima. Se chiamo questa A, allora è proprio un problema agli autovalori standard ogni codice usa le sue librerie per calcolarsi questo problema agli autovalori generalizzato vedrete che ci sarà una linea all'interno del nostro codice che non si occupa di fare questa decomposizione semplicemente dà MATLAB le matrici K ed M e chiede a MATLAB di calcolare gli autovallori di autovallorio. Per noi, indico come utenti, è una cosa semplicissima. Però ci sono delle probabilità fondamentali. Da quello che io ho fatto discendono queste probabilità che sono di importanza drammatica per quello che dobbiamo fare noi. Tutti gli autovettori sono m ortonormali. Cosa vuol dire? Che se premoltiplico e postmoltiplico per lo stesso autovettore ottengo 1. Se invece prendo due autovettori diversi ottengo 0. Cioè non sono ortogonali, ma si dice m ortogonali perché sono 0 solamente se in mezzo al prodotto scalare c'è la m. Quindi l'autovettore xi per m per xj da 0 invece se sono lo stesso autovettore si ottiene 1 questa probabilità è importantissima e discendere da quella relazione là anche che una probabilità simile vale per k se faccio xi k xi ottengo omega quadro cioè l'autovalore quadrato se Se invece faccio x, i, k, x, j con i diverso da j io tengo 0. Quindi le proprietà sono estremamente simili per le due matrici m e k, l'unica differenza è che è normalizzata in m, invece per k vengono fuori le frequenze proprie. Ecco, provate a tenere a mente queste proprietà perché sono fondamentali nel piccolo passaggio che adesso dobbiamo fare, che spiega quello che dicevamo prima, che spiega tante cose e spiega anche il funzionamento sostanzialmente dei mezzi, di molti mezzi, di tutti i mezzi vivanti. Ok, allora quanti autovettori abbiamo in generale? E se la matrice è di ordine n, si possono calcolare n autovettori. Nessuno mai si metterà a calcolarli tutti quanti, perché è difficilissimo ed è computazionalmente estremamente oneroso. È facile calcolare, ad esempio, alcuni autovettori, quelli associati alle frequenze più basse, che sono in generale quelli che ci interessano. Questo è abbastanza facile da fare, anche abbastanza rapido, ed è proprio quello che chiederemo a Matela di fare, calcolami i primi 5, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i primi, i prim questo è un esercizio teorico che facciamo, di conoscere tutta la base costituita dagli autovettori, perché gli autovettori costituiscono una base per lo spazio di dimensione n. E esprimiamo, quindi, decidiamo di esprimere il vettore u come una combinazione di questi vettori della base, cioè degli autovettori. Quindi facciamo una sommatoria su j che va da 1 ad n dove è grande la dimensione della matrice di questi coefficienti che chiamiamo uj sono degli scalari, vedete, sono degli scalari uj che dipendono dal tempo che si chiamano i coefficienti modali proprio perché sono i coefficienti che moltiplicano i modi allora uj che moltiplj, quindi un coefficiente che moltiplica un modo, un altro coefficiente che moltiplica un altro modo e così via. Ok? Immaginiamo di farlo. E se lo facciamo con lo spostamento è evidente che per l'accelerazione varrà questa decomposizione, in cui usiamo gli stessi identici modi e scriviamo due punti sopra il coefficiente j, esattamente come abbiamo fatto con gli elementi finiti di prima. Ok? Allora, teniamo a mente questo, ecco, no, adesso questo dopo lo ritorneremo, no, come mai non... Cos'è? Eccolo qua, scusate, ho saltato questo. Allora, facciamo un riassunto, stiamo risolvendo quella roba là, con l'equazione, abbiamo decomposto lo spostamento e accelerazione nella base modale, abbiamo scelto di fare questo, e utilizziamo adesso queste proprietà, le proprietà di ortonormalità. Ok, cosa facciamo? Due cose, prendiamo queste cose qua e le sbattiamo qua dentro, e poi premoltiplichiamo per il generico modo x, i. Facciamolo, vediamo cosa viene. Premoltiplico per x, i. Ok? Per x, i. E poi scrivo quello che ho là. Quindi avremo m. Poi, per u2 punti, ma u2 punti lo scrivo su questa base, quindi avremo sommatorio su j, attenzione avrò la sommatoria su j di xj che poi moltiplica u2 punti j, non scrivo la dipendenza dal tempo perché vale sempre poi più lo faccio anche per k chiaramente quindi avrò xi che moltiplica k che moltiplica la sommatoria su j di che cosa? di xj che moltiplica uj non perché l'ho messo sia il numeratore e poi ci sarà anche il termine noto stessa identica cosa quindi avremo scusate la x i che moltiplica che moltiplica ok? io ho dimenticato i trasposti e li metto subito bene quindi ho usato la rappresentazione nella base modale e ho semplicemente riscritto la prima equazione là sopra dopo averla premoltiplicata per xi ok? chiaro quello che ho fatto? adesso utilizziamo le proprietà di ortonormal, xd per m per xj dico che è uguale a delta ij. È uguale sempre a zero tranne quando i è uguale a j. Ma siccome poi devo fare la sommatoria su j di delta ij per 1,2 j, scriviamo. Poi sommatoria su j. Qui invece di avere delta i j avrò omega quadro j per delta i j per 1 j. E questo coefficiente invece lo chiamo banalmente fj. La proiezione che cos'è? Cosa rappresenta l'ultimo termine? Fj trasposto F rappresenta la proiezione delle forze nodali sul modo. Quindi questo Fj viene... Scusa, perché l'ho chiamato J? Perché c'è solo i invece di... perché qui ho scritto malissimo, è fi, perché questo è l'indice i. Ora, che cosa otteniamo qua? Sommo su j, ma questo coefficiente è sempre zero, tranne quando io uguale a j, quindi ottengo ui. Ui due punti, più omega, stesso discorso, sempre uguale a zero, tranne quando io giungo all'I, quindi avrò omega I al quadrato per ui uguale ad fI. E questo è vero qualunque I, perché la scelta di x è stata arbitraria, quindi lo posso fare per qualunque altro vittorio. Cosa ho ottenuto? Quindi delle equazioni differenziali e ordinarie completamente disaccopiate fra di loro. Prima aveva un sistema, questo sistema di partenza, però là in alto, era un sistema che accopiava tutte le incognite Accoppiava tutte le incognite, quindi era difficilissimo da risolvere. Ho rappresentato nella base modale il vantaggio che queste equazioni sono tutte equazioni scalari, completamente disaccoppiate una dall'altra. Cioè la prima è indipendente dalla seconda, dalla terza, dalla decima, dalla miglionesima. E quindi è infinitamente più semplice risolvere, perché questa lo sappiamo risolvere tutti. Questa è la classica, è la prima equazione che ci hanno insegnato, non so se si chiama ancora da voi, meccanica razionale, l'oscillatore armonico, la cui risposta va, se abbiamo una, se analizziamo questa, con una forzante armonica, la risposta di questa roba parte da uno, poi va all'infinito quando la frequenza della forzante è uguale ad omega i e poi dopo riscende. Questa è la risonanza. Quindi, cosa vuol dire? Questo è un esercizio accademico, quello di conoscere tutti quanti i modi, ma questo risultato finale è valido anche se ne sappiamo pochi, perché ammettete di sapere solamente il primo modo, allora conoscete questo, e se questa operazione la fate solamente per il primo modo, beh, ottenete la prima equazione. Le altre non le conosciamo, perché non siamo bravi a calcolare gli altri modi, però sappiamo che la prima equazione è questa. Una volta che so la frequenza propria del primo modo e poi conosco il primo modo per poter fare la proiezione della forza sul primo modo, posso costruire questa equazione scalare che riguarda il primo modo. Se io so che forzerò la struttura con una forzante a questa frequenza, so per certo che risponderà a questa frequenza e basta. Che la risposta di questo è caratterizzata da VUM, dal picco, no? In corrispondenza di omega uguale omega I. Quando omega nella frequenza della forzante è uguale alla frequenza propria, la risposta spara ad infinito. Tutti i mensa vengono costruiti proprio per lavorare attorno ad un modo, ad un modo particolare. Bisogna capire volta per volta qual è il modo che ci interessa. Bene, qui abbiamo detto una storiella, in realtà nelle strutture reali c'è un po' di dissipazione e questo sarà importantissimo la settimana prossima, quando parleremo di dissipazione e faremo il primo esempio di dissipazione, dissipazione termoelastica, cos'è? In generale a questa equazione si aggiunge, in base a qualche informazione sper sperimentale ormai siamo abbastanza bravi anche a calcolarlo però molto spesso è meglio avere un'indicazione sperimentale un termine che è proporzionale alla velocità modale e che viene espresso sempre in questa maniera come omega la stessa omega modale diviso per un coefficiente che chiamiamo Q e che è noto come fattore di qualità, quality factor. Quality factor quando vi dico a un membro che frequenza lavora, che Q ha? Finito. Le tue informazioni che tutti vogliono sapere di un oggetto vibrando che frequenza ha e che fattore di qualità. Perché la frequenza vi dice chiaramente l'informazione di qual è il regime a cui state pensando e di quanta energia dovete mettere dentro il MEMS per tenerlo in vibrazione. Perché se non ci fosse dissipazione qui non ci sarebbe necessità di energia. Vi date una condizione iniziale, quello va avanti all'infinito a vibrare. C'è un po' di dissipazione, nei MEMS pochissima, ma c'è, e allora bisogna quantificare i cubi per sapere quanta energia si mette sulla dissipazione e torneremo la settimana prossima. Oggi ci dimentichiamo della dissipazione. Ok, allora, quindi, avendo visto questa cosa, andiamo a vedere effettivamente che cosa succede all'interno del codice. Aprite il codice modal. Ecco. Nel display se troverete alla pagina 70, 4, 2, 1, questo esercizio che qui ho riportato, che è un pezzo della mia presa. Vogliamo calcolare i modi e le frequenze proprie in due casi. Uno è il più semplice, cioè non è vero che che è più semplice però è mi ho gelato le informazioni nel testo allora lanciamo modal e scegliamo come fai v input model bim cc cosa vuol dire clamped clamped cioè è una trave che è incastrata agli estremi. Questi saranno effettivamente i modi di vibrare che troverete. La trave è incastrata a sinistra e incastrata a destra. E vibra. E questi sono i primi tre modi di vibrare che dovremmo trovare facendo andare il codice. Allora, apriamo e vediamo il primo ris, ecco questo dovrebbe essere il risultato che trovate in Gmesh. Ha funzionato? Adesso sto ancora facendo, sto estraendo. Ah ok, allora aspettiamo un secondo. Intanto io apro il file diput per vedere cosa c'è dentro. Ok. Ok. Cominciamo a vedere... avete provato a farlo andare e è arrivato più in fondo. Dovrete aprire automaticamente Gmash con l'immagine di una travinzessione. Siete riusciti? Problemi? Stavo ora spessando di sapere una flip. Ah ok! Anche voi. Faccio run con il codice modal. Modal, sì, modal, esattamente come la st stessa di lo lanci, vi richiedo un file input e prova a lanciare questo no? adesso modal che è a livello superiore nella rete superiore non ho capito dove scusi
In realtà tu non vedrai magari la capacità, leggerai una corrente, insomma, ma le due cose sono legate l'una all'altra. Vedrai a schermo una cosa di questo genere, qui hai la omega della forzante, cioè l'omega della via C che stai facendo variare e qui metti addirittura di capacità, qui parte da da da da da da da da da darostatica, distribuita su tutto l'elelettro faremo in una delle ultime lezioni questi conti dove sperimentalmente come l'unica cosa che vuoi calcolare sono le frequenze proprie perché corrispondono ai picchi tu avrai uno spettrogramma che sarà fatto così picco, picco picco, ecco questi sono le varie frequenze proprie della tua strada. Io avrei una domanda molto specifica. Se a livello teorico noi abbiamo un solo picco di risonanza, ma in realtà a livello sperimentare ne abbiamo più di uno molto vicini. Questo è il problema tipico. Molto spesso ci sono dei picchi spuri che vengono dall'elettronica. In alcuni casi, purtroppo, per nostra fortuna, sono i problemi che abbiamo più frequentemente, il comportamento non lineare dell'oggetto fa sì che più modi si combinino fra di loro. E allora ci sono quelle che si chiamano le risonanze interne e più modi che a priori dovrebbero essere indipendenti si mischiano. Vi faremo vedere, mi fermo magari in seminario su questo. Quando c'è un rapporto intero fra le frequenze, avete un modo che ha frequenza 1 e un altro modo, ad esempio, che ha frequenza 2, poi pensate di citare solo il modo 1, ma parte anche il modo 2, per le non linearità che ci sono nel sistema, e questa si chiama risonanza 1 a 2, l'abbiamo vista tantissimo nei giroscopi, le risonanze 1 a 2. Non so se avete mai sentito parlare dei DRG, che sono i disk resonator gyroscopes, che sono fatti per avere due modi praticamente identiciici rotati di 45 gradi, lì ci sono delle risonanze uno a uno cioè voi citate un modo e si mette a partire il modo a 45 quindi a volte sono delle porcherie di questo genere a volte è semplicemente che non risponde solo alla meccanica, vedete dei modi elettronici anche e i nostri amici in S in ST in generale sono bravissimi a capire subito quali sono i modi meccanici e quali sono i modi elettronici. Io non lo so, ma loro capisco lì immediatamente, in base al modo 40 Hz, in modo qui, sanno la loro elettronica e che modi genera. La loro elettronica e i modi genera. Ma queste sono cose che impari sul campo ah perché lo scorso semestre ho fatto il laboratorio con il professor Ghisi ah, claro ci aveva dato quel MEMS con diciamo per il testa fatica di film metallici sottili e appunto il problema era questo quindi noi abbiamo fatto un report cercando di indagare questi quattro picchi e poi lui cosa va dentro? il problema è che anche lui questo meme si era un po' incasinato ce ne sono tanti, vi faremo vedere degli esempi loro spesso hanno dei fenomeni meccanici che non capiscono allora vengono da noi a cercare aiuto per capirli magari ci sono tre modi nell'ultimo esempio era un giroscopio in cui c'era una combinazione di modi, uno a 20 kHz, uno a 50 e uno a 90, o uno a 20, proprio che devono essere l'insomma, quindi la frequenza del terzo era la somma dei primi due, quindi se 20, 50, 70, qualcosa del genere, non mi ricordo più, ogni volta che ci sono queste combinazioni di interi che generano una frequenza come somma delle altre due, allora possono essere questi fenomeni di risonanza interna. Siccome ci sono tantissimi modi in una struttura meccanica, purtroppo più la struttura è complessa, più queste porcherie sono probabili. Però il vostro problema secondo me era completamente diverso, più da elettronica, e allora si va a chiedere a Giacomo Lanfelder che cosa succede, oppure a quegli stessi IST. Comunque, prima di tutto, per ogni device loro vanno prima a fare lo spettro lineare, cioè a calcolarsi tutte quante le frequenze lineali con queste dengue. Siete riusciti a far andare? Allora, vediamo un po' quali sono le informazioni e quali sono le piccolissime differenze rispetto al vecchio codice, anzi direi che è una versione insentificata. Questo è il file input relativo alla beam doppiamente incastrata che ha bisogno, come sempre, anche lui di una mesh. La mesh è generata nella stessa maniera. Avete un rettangolo, vi ricordate l'avete fatto anche l'altra volta più o meno per fare il calcolo della rigidezza della fold quindi allora materiale qui ha tante cose in realtà perché poi l'ho utilizzato lo utilizzerò in futuro anche per fare la termoelasticità quello che vi interessa a voi sono i primi tre yang, quassone e densità fittizi. Poi però una volta ci mettiamo lì veramente a vedere cosa viene fuori con l'unità MEMS, questi numeri veri, cosa vengono, se no rimangono sempre numeri fittizi. Poi qui avremo anche i parametri termici, capacità termica K, alfa, ma li vedremo la prossima volta. Poi stessa cosa dell'altra volta, dobbiamo associare un solido a un materiale, quindi avremo un elset fisico 3 che rappresenta la superficie associato al primo materiale. Poi abbiamo delle condizioni a contorno che bloccano questa linea e bloccano quest'altra linea nella prima e nella seconda direzione. Evidentemente questa qui è l'elset fisico 1, questo è il set fisico 2, la superficie è l'elet fisico 3, quindi sul primo esset fisico in direzione 2 mettiamo spostamento 0, sul secondo 1, 2, spostamento 0. Quindi l'abbiamo incastrato a destra e a sinistra. Quindi questo è l'input. Il model va, non so se avete preso un po' di tempo per andare a ripassare la stessi ma la prima parte di model è la fotocopia di quello che fa la stessi cioè legge la mesh fa l'enumerazione delle incognite alloca la matrice sparsa l'unica differenza è che adesso non alloca solo a k ma all alloca anche una seconda matrice, che è la matrice di massa, la m. E poi lancia la procedura di assemblaggio. La procedura di assemblaggio è copia-incolla quella di elasticity, con la piccola differenza che adesso si calcola anche la matrice di massa elementare, e poi assembla la matrice di massa. Vedete che però le istruzioni sono identiche a quelle della matrice di rigidenza. Quindi l'unica cosa che è cambiata è proprio quella routine elementare che abbiamo calcolato prima. Va bene, alla fine dell'assemblaggio ha tutto quello che ci serve, perché se vogliamo calcolare gli autovalori, la matrice K abbiamo la matrice M utilizziamo questo comando EX se volete sapere che cosa fa andate sulla riga di comando di MATLAB e scrivete help EX e lui vi dirà, vi inonderà di informazioni su che cosa fa tutte le opzioni eccetera eccetera quindi imparate come si può fare in particolare se gli si danno due matrice vedete qua se AX su A e B risolve il problema agli autovalori generalizzato che è proprio quello che serve a noi quindi gli diamo la matrice K che è semidefinita la matrice M che è definita il numero di autovalori che vogliamo, i più bassi, perché gli diamo come opzione SM, gli smallest. Quindi lui si calcola i 5 autorettori associati agli autovalori di valore più basso. Quindi i primi modi di vibrare. Poi fa la radice della D per calcolare le ω e lo scrive lo scrive a schermo guardiamo insieme come come va cosa mi fa vedere se avete chiesto 5 auto vettore avete 5 output la mesh verdina è quella originaria quella colori e la deformata modale ora deformata modale a volte viene fuori scala allora bisogna andare a risettare la scala con il solito trucchetto questo però nei modi è molto importante perché bisogna fare options si va sulla vista che noi vogliamo ad esempio la vista 0 e si va in aspetto e in aspetto c'è il displacement factor che può essere modificato ad esempio se faccio 10 vedete che è enorme quindi amplifico lo spostamento devo calibrare un po' a mano purtroppo ci sono delle tecniche anche automatiche, ma non ho avuto voglia di farlo, calibrate con questo fattore l'ampiezza. Allora, questo è il primo modo di vibrare a flessione, a flessione perché la flessione domina. Un risonatore di questo genere lavorerà senza il minimo dubbio a questa sua frequenza. I modi, man mano che saliamo, diventano sempre più nervosi, hanno sempre più onde. Questo è il secondo modo flessionale. E poi così via, terzo modo flessionale. Attenzione, non ci sono solo i modi flessionale ma ci sono forse l'abbiamo trovato vedete che adesso è completamente coperto io cerco di snidare sarà un po difficile snidare il modo amplifico molto non basta ecco sì che è bastato è andato'altra parte un po un po orrendo allora vediamo non so vediamo come è che potrei fare dovrei togliere la mesh però non mi ricordo come si fa dovrei togliere la mesh sovrastante è un modo assiale cioè c'è una deformazione assiale non è una formazione frizionale che lo fa uscire. I modi assiali che hanno delle deformazioni lungolaste sono generalmente a frequenze superiori agli altri. Allora, che cosa vi propongo di fare? Ok, avete la stima delle Omega, sputata fuori da Matlab, e quello che vi chiedo di fare è di andare a verificare la formula 412 che se voi andate su internet potete benissimo provare e dire Bending Frequencies of Clam Clam Beam vi vengono date queste formule qua. Provate a vedere almeno la prima. Le sono formule analitiche notissime. Che cosa abbiamo? La frequenza di librare, usate matelo stesso per farlo, fate uno scrittino. L'anda sono queste cose, L è la lunghezza della trave, dove la trovate? All'interno dello geo. L'ho della trave, rho è la densità, quella cosa che abbiamo messo 1, e il modulo di Young, h è lo spessore della trave. Provate a vedere se funziona, quanto siamo lontani dalla verità, spero che siano giuste le poveri, io le ho copiate ma non le ho verificate. Vediamo. Magari la cosa più semplice che potete fare e scrivere all'interno di bim della stessa modalità l'imm a un po più la formazione però voi perché guardatevi che la validazione è sempre la parte più importante di un codice così magari siete costretti ad andare a vedere come è definita la geometria, tutte queste cose qua. Oh, hai qua da... a Sì Wow, è andata bene. A me è venuto 1.02.60C, quindi direi che vorremmo dire lì. Provate a... La formula si trovava in Grazie. qual'era la formula? la formula è... trovo una che è una su due pi greco radice quadrata di e per i su ok grazie poi cambiando lambda potete vedere tutti i modi questi sono i primi modi flessionali valgono queste ponuli solo per i modi flessionali Ombra Sì, è 10 Rho è 1 Quindi se sostituite mi dovrebbe essere 4.73 diviso 10 moltiplicato square root di tutto questo e 12 e rho diviso 12 e Rho diviso 12 e Rho era 1 ah giusto le parentesi quadra non le prende se qualcuno è riuscito provate a fare il secondo modo ah è il quadrazzo andiamo a vedere se funziona anche per il secondo modo rimane tutti identi basta cambiare la 7,853 2,81 va bene ma che ci questo ci serve per fare qualche commento. Allora, vediamo un secondo. Queste sono le frequenze calcolate, no? Queste sono, segniamoceli da qualche parte, 1,0265 e 2, sono i primi due modi, 2,7728. Se andiamo a vedere le approssimazioni analitiche che io gli ho fatto calcolare, vengono 1,0212 e 2,814. Anche se il secondo modo fa un errore del 3x1000 però il primo modo era molto accurato questa è una regola, i primi modi sono sempre più accurati e man mano che salite con i modi diventa meno accurato magari è un effetto della mesh? Boh, non so, vediamo a vedere come era fatta la mesh lc1 e v diviso 4, credo che fosse già abbastanza fine, però esageriamo, facciamo V diviso 8, vediamo se cambia qualcosa. Poi la sorgente del terrore potrebbe essere che queste sono delle stime fatte con la teoria delle travi, noi stiamo usando la meccanica dei solidi che è diversa. Allora questa è la mesh, quindi molto più fine e vediamo se lancio modal, se cambia un po' la frequenza. Io non mi aspetto di vedere grandi cambiamenti nelle prime frequenze proprie però l'ha fatto, non ha avuto grandi problemi, andiamo a vedere 2,77, niente, 2,81 vedete, l'errore quindi è rimasto invariato, i due modi primi erano già accuratissimi l'errore è rimasto invariato è dovuto al fatto che stiamo implementandoorie diverse, teorie delle trami opposte alla teoria della meccanica dei solidi, senza poi trascurare che adesso qui non stiamo ipotizzando che il problema sia plane strain e non è detto che effettivamente sia veramente plane strain, insomma è un'approssimazione del tridimensione. Allora, questo vi lascio un paio di esercizi, perché come vedete finché le racconto io le cose hanno un senso, poi mettersi a farle ha sempre un senso diverso, questo esercizio è ripetere la stessa cosa per una cantilever, un esercizio molto semplice, cioè che cos'è la cantilever? È un linguaggio tecnico che tiene una trance che è incastrata solo da un lato, cioè è incastrata solo qua. In teoria potete fare una micromodifica all'input per la cl la clan clan per ottenere un cantilever. Non per ottenere un cantilever. Basta togliere l'incastro qua. Basta togliere l'incastro qui e diventa un cantilever. Ecco, poi fate andare il codice, vedete i risultati e fate un esercizio di ricerca su internet delle sequenze proprie di una cantilever. Di una trave cantilever. E parilivere e paragonate quello che mi viene. Quindi sappiate che su internet ormai si trova qualunque cosa. Allora, l'altro esercizio più interessante è l'uso delle simmetrie, perché anche questo è sempre drammaticamente importante, soprattutto quando i merici sono molto complicati e sono tridimensionali, si cerca di risparmiare un po' introducendo delle simmetrie. Allora magari a qualcuno viene l'idea di dire ok, ritorniamo la nota che è il tab e invece di discretizzarla tutta quanta dico che è una struttura simmetrica, la taglio a metà, tolgo questo e metto dei carrelli, che dovrebbero rappresentare le condizioni di simmetria. Allora la domanda che vi faccio è, innanzitutto provate a farlo, facilissimo, basta fare garanzia. Proviamo a farlo assieme, questo è quasi divertente. Vediamo quello che si deve fare' diversa cambio solo le condizioni al contorno invece cos'è che devo modificare sulla prima linea devo bloccare gli spostamenti in tutte e due le direzioni, nella seconda devo bloccare gli spostamenti in direzione 1, per cui chiudo la cosa qua, facciamo così, mi commento se no se ne perdo dopo, faccio 2, 1, 0 e 2, 1, 0 e chiudo ho cambiato quindi le condizioni di contorno senza modificare la lunghezza della trave quindi paragonerò i risultati analitici con una trave di lunghezza 20 perché se se 10 rappresenta la metà, la lunghezza della trave dovrà essere 20. Vediamo che cosa succede se faccio questo giochino, che cosa mi aspetto. Innanzitutto se funziona, facendo le cose così, mi rischia che non funzioni. Ecco. Allora, il primo modo, mi dice che il primo modo è questo. Il secondo modo è quest'altro. E ovviamente tutti i modi che trova rispettano questa condizione, che nell'estremità di destra questo rimane piatto, sostanzialmente verticale, si può spostare solamente in verticale. Cioè, se andiamo a paragonare con quello che trovavamo prima, che è rappresentato qua all'interno delle nostre slide, che cosa succede? che se cerchiamo di risolvere una trave che ha come condizione che questa sezione rimanga diritta e si possa spostare solo in verticale certamente risolveremo in maniera corretta alcuni autovettori questo perché rispetto a questa condizione ah no, scusate, questo no però questo qui no neanche questo, dire la verità questo prevede una sezione che si inclina qua quindi se forziamo la simmetria questo modo lo perdiamo completamente calcoliamo solamente i modi simmetrici riusciamo a calcolare i cosiddetti modi simmetrici se questi sono i primi tre modi non me lo ricordo il secondo lo perdiamo, ilzo lo perdiamo. E ricominceremo a trovare i modi a partire da quelli più alti, quando c'è un numero pari di onde e la sezione e mezzo resta verticale. Quindi attenzione a utilizzare le condizioni di simmetria perché bisogna essere sicuri che non solo la geometria sia sim vi potete anche divertire a guardare abbiamo preparato se poteteciare modal, dovrebbe essere abbastanza potente il vostro computer. Lanciate modal e andate a prendere modal gyro. Avete già il cantilever, come vedete, modal cantileilever è già lì pronto, dovete solo andare a verificare. Il mod al gyro è abbastanza pesante come analisi, se mi ricordo. Contiene già 13.000 di nemico di tipo. Ecco, il mod al gyro ha, se non mi ricordo male, una definizione di proprietà sensata apriamo, andiamo a vedere com'è l'input, è il punto M perché non lo apre? perché non sta aprendo il mod al gyro? ecco sì allora vedete per la prima volta troviamo delle unità mems corrette cioè yang 160 gigapascal per il poli silicio. E poi lì c'è scritto 1,6 per 10 alla quinta. Vabbè, il Poisson è un numero puro, quindi 0,22. La densità è 2330 kg al metro cubo, eppure c'è 2,32 per 10 alla meno 15. Pa zero questa è un'altra storia è uno sforzo iniziale che va bene non serve nel codice non è utilizzato come mai ci sono questi numeri? le dimensioni del gyro saranno certamente assegnate la creazione della mesh gyro se entrate nel punto G non è così facile non è così ci vuole un po' di tempo modal gyro ecco se andate a vedere è un file abbastanza lungo che definisce tutte le proprietà non so se era Valentina Anziga che aveva fatto questo input non so, un sacco di cose e la domanda a cui dovete cercare di rispondere a questo esercizio è se le dimensioni sono date in micron quale può essere una scelta ragionevole delle altre unità? In questo caso dobbiamo segnare la lunghezza e l'abbiamo fessata in linea. Poi abbiamo sempre bisogno di dare il modulo di Young, che è una pressione, una forza divisa unità di superficie. Ma se la lunghezza l'abbiamo fessata, allora dobbiamo fessare unità di forza. Quali unità di forza abbiamo usato? Newton, mega newton, ok. Poi abbiamo bisogno di dare la densità, e la densità in unità SI si dà chilogrammo al metro cubo, però ormai abbiamo abbandonato il metro, abbiamo deciso di usare il micron se fosse chili al micron cubo se la densità di partenza è 2300 allora bisognerebbe ovviamente a ogni micron cubo ci sono ben pochi chili quindi dobbiamo prendere 2300 e dividere per 18 allora può essere? eh essere. Allora, in realtà in generale non si usa il chilo. Se ho usato il chilo è solo per mettere un po' in difficoltà, si usano delle unità di massa più piccole. Ora adesso se abbiamo fissato la lunghezza, abbiamo fissato la massa e l'unità della forza dovrebbe essere conseguente. Allora, provate quindi a chiedervi, non adesso, con tranquillità, qual è l'unità della forza dovrebbe essere conseguente. Provate quindi a chiedervi, non adesso, con tranquillità, qual è l'unità di forza coerente? No, voi direte non basta, devo fissare anche il tempo. Se fesso il tempo, allora anche la forza è fissata, oppure viceversa, fesso la forza è fissato il tempo. Una volta che avete fissato queste quattro cose, devono essere coerenti fra di loro, non sono tutte libere. La forza è una massa per lunghezza nel tempo quadro, per cui sono legato fra di loro. Forza uguale a massa per accelerazione. Per cui se fissate la massa, fissate la L, fissate il tempo, laza viene fissata e così via chiedetevi quindi quali sono le unità con cui abbiamo dato questo compito comunque andiamo a vedere il risultato in modo giro adesso arriverà un po' ecco si è calcolato i miei modi beh non si vede niente questa è la geometria è la geometria del gyro queste sono le molle folded è molto semplice però basta per spiegare il funzionamento di un giroscopio di un semplice giroscopio adesso vi dico come funziona ci sono questi quadratini che sono gli ancoraggi cioè l'oggetto è ancorato sul strato qua. Poi, le molle folded sono quelle che abbiamo già visto, sono gli elementi, non ci sono cerniere, non ci sono svincoli, non c'è niente, le rigidenze e le cedevolezze sono date dagli elementi elastici, quindi facciamo delle molle lunghe ma ripiegate per non occupare troppo spazio, come abbiamo visto l'altra volta, ci sono delle formule che vi danno le rigidezze di queste molle. Ci sono quindi due set di molle, che sono questi due e questi due, che danno rigidezza in direzione orizzontale, ci sono due set di molle identiche che danno rigidezza in verticale. Poi c'è il cosiddetto shuttle, il rotore, che è la massa in mezzo che è collegata alle molle, che quindi è libera sostanzialmente di muoversi nel piano come vuole, si può traslare, può ruotare, per effetto della celevolezza di queste molle. Che cosa rappresentano invece questi piatti rettangolari? Sono degli elettrodi, sono lì fissi, sono messi per bellezza, sostanzialmente a questo livello non entrano assolutamente nell'analisi, perché è un'analisi sicuramente meccanica. Questo dovrebbe essere un semplici, lo scopo è capacito di guardare come funziona. C'è un cosiddetto modo di drive. Innanzitutto andiamo a vedere come sono fatti i primi modi, i primi modi di vibrare. Vediamo innanzitutto le frequenze. Vedete che le prime due sono praticamente identiche, sono molto poco diverse tra di loro. Andiamo a vedere che forma hanno i modi di vibrare. La geometria è talmente simile in orizzontale e in verticale che ci aspettiamo dei modi uguali in uno spostamento orizzontale e in uno spostamento verticale però dobbiamo modificare l'entità del fattore moltiplicativo quindi andiamo su aspetto e moltiplichiamo vediamo se 100 basta no, vedete è così piccolo che devo andare... ecco, 100.000. Ho messo... scusate, 10.000. Allora, qui è difficile un po' vedere perché si è sovrapposto, però vedete l'apparizione di questa massa rossa. La massa rossa è quella sottostante e si è tutta spostata verso sinistra. Una specie di traslazione rigida. Vete come si sono deformate le molle, però manca totalmente lo spostamento in verticale, questo è ottimo, quindi la massa si sta spostando verso sinistra. Se fate la stessa cosa per il modo numero 2, dovete modificare alla stessa maniera il fattore moltiplicativo, ancora 10.000, e vedete la stessa cosa in verticale. Quindi il secondo modo è un modo di traslazione verticale. Proviamo a vedere cos'è il modo 3, chissà se è una rotazione, andiamo innanzitutto qui, dobbiamo cliccare sul modo 3. Poi andiamo a modificare l'aspetto e anche qui lo moltiplichiamo per 1000, 1000 è troppo poco, 10ila, diecimila è tanto. Vedete che il terzo modo è un modo di rotazione, è un modo di rotazione del piano. E così via, potreste calcolare tutti quanti i modi e andare a vedere cosa sono. Ma ritorniamo quindi sul primo modo di vibrare. Questo è il cosiddetto modo di drive, almeno io immagino che possa funzionare così cioè questi elettrodi danno una via C questo e questo qui e mettono in movimento e ci sarà un primo modo di vibrare danno una via C esattamente alla frequenza del primo modo e questo comincia a vibrare così questo viene messo in movimento di vibrazione ovvero il giroscopio avrebbe una frequenza intorno ai 20 kHz in genere, come vi ho detto. Adesso immaginate che ci sia una rotazione con un componente ωZ, cioè questo oggetto è all'interno della Ui, la Ui sta ruotando, sta ruotando così, quindi l'ω è un vettore che esce da piano, ha solamente la componente ωZ. Allora qui c'è la famosa forza di Coriolini, no, c'è la famosa forza di Coriolini, che è una forza apparente che mi dice che all'interno del sistema non inerziale, questo device è attaccato al package che si sta muovendo, quindi noi siamo all'interno del package, siamo in un sistema non inerziale, c'è una forza apparente che è data da due volte, e il due ce lo dimentichiamo, anzi o meno due, omega, velocità relativa, quindi omega, velocità relativa, forza, viene una forza verticale, la forza di Corioli genera una forza verticale. Quindi cosa succede? Che questo qui, con il petto dell'omega, comincia a oscillare anche in direzione verticale. Questi elettrodi misurano la variazione di capacità per effetto di questa vibrazione verticale e vi dice, a ritroso, ottenete l'identificazione dell'omiga. Questo è un giroscopio. Poi, i veri giroscopi, li vedrete, sono molt sono molto più complicati perché servono per misurare tutte e tre le componenti di velocità angolare questo è fatto solo per la per la omega z non funziona però se avete omega y omega y questa semplice struttura non funziona deve essere più complesso però vi dice chiaramente come... ecco, è un giroscopio quindi l'omega z si basa sulla presenza di due modi di genere cioè la stessa frequenza, uno orizzontale e uno verticale, perché è ovvio che io voglio dare poca tensione via c, permettono di vibrazione risontale, ma voglio che una piccola omega eccidi tantissimo il modo verticale e questo avviene se siamo alla frequenza giusta cioè se il modo orizzontale continua la velocità relativa varia con la omega del primo modo e quindi anche la forza in verticale varia con laω del primo modo quindi se il secondo modo verticale ha la stessa frequenza si eccita tanto per questo servono i due modi si dicono de genere con frequenze molto semili fra di noi ok? questo è un po' il principio di funzionamento e di un banalissimo di un giocattolo scusi volevo solo informare che se si tiene premuto il pulsante sinistro su displacement factor e si muove il mouse da sinistra verso destra si mettono i valori che consiglia Gmesh perché a me il primo viene 2061,49 allora andiamo tools, options andiamo su aspetto displacement factor, tengo schiacciato e si muove il mouse verso destra, da sinistra a destra da sinistra verso destra e lo aumenta così, no veramente sì sì lo aumenta però non è che mi dia quello consigliato, mi dice quello che vuole lui no perché a me il primo che esce è carino cioè 2061,49 esce carino anche il secondo che è tipo ah ma no ok sono multipli di quello del primo che arriva cioè sono multipli di 2061 e 49, peccato va bene quindi ogni struttura che voi avrete per prima cosa noterete che tutti fanno l'analisi modale, in ST usano ANSYS, un codice molto costoso per farlo, durante le tesi qui si usa più COMSOL che altro, che è quello che imparerete a usare con Matteo alla fine di... Questa è la prima cosa, l'analisi modale, che per i MEMS è la più importante. Facciamo una pausa, ci sono minuti, ma mi rimane un po', e poi facciamo la seconda parte, che invece è la parte più, se volete, divertente, ma meno utile per i MEMS, e cerchiamo di finire abbastanza presto.
Siamo al sistema semi-discretizzato. Ok, che è questo qui. Ha questo sistema. Ovviamente la cosa che si fa più frequentemente, non con i MEMS, ma con le strutture normali è prendere questo sistema di scrittati e cercare di integrarlo nel tempo. Cioè a partire dalle condizioni iniziali che sono assegnate sul spostamento e sulle velocità dobbiamo seguire l'evoluzione di questo campo degli spostamenti. Questo l'abbiamo già visto in realtà, qualcosa di molto simile, la primissima lezione quando abbiamo fatto la diffusione 1D sulla sfera. Abbiamo ottenuto un'unicazione discretizzata in questa maniera, solo che non avevamo le derivate seconde dell'oppostamento, in quell'esempio c'era la derivata prima della temperatura, perché era un problema parabolico, non era un problema parabolico. E quindi vedremo qualcosa di abbastanza simile. L'approccio più semplice è prendiamo l'intervallo di tempo che ci interessa, l'intervallo di tempo che va da 0 a TF, e lo dividiamo in time step, in passi temporali. Nel caso più semplice che facciamo noi è il caso più banale, passi di tempo tutti di lunghezza uguali, quindi prendiamo un incremento che chiameremo delta t, l'incremento di tempo, t0 sarà 0, t1 sarà delta t, tm sarà m per delta t. E vogliamo calcolare la soluzione a ciascuno di questi istanti temporali. Cioè, abbandoniamo l'idea di trovare una soluzione continua, ci accontentiamo di trovare una buona stima, la migliore stima possibile della soluzione a questi timestamp. Allora, è evidente che... Allora, supponiamo, sc evidente che supponiamo di conoscere tutta la soluzione al tempo n, dove questo p di cn vuol dire il valore di quella variabile al tempo tn. Cerchiamo di formulare il problema sul passo. Cos'è il problema sul passo? Data la soluzione al tempo tn, quindi se conosciamo un punto n e conosciamo un due punti n, quindi la stima dello spostamento, velocità, accelerazione al tempo tn, vogliamo fare la transizione al tempo tn più 1 e quindi calcolare la miglior stima possibile di un n più 1, u punto n più 1 e u due punti n più 1. Chiaramente dobbiamo conoscere l'evoluzione dei carichi, cioè quella F che appare a termine destra deve essere assegnata, i carichi applicati alla struttura devono essere noti. Allora, vedremo due algoritmi, il primo è il più facile, l'esempio su cui giocheremo dopo è basato proprio su questo metodo, ma vedremo che è un elemento di una famiglia molto più generale. Il metodo è le differenze centrali, che è un metodo esplicito e dopo commento su cosa vuol dire il metodo esplicito. Allora, cosa si fa? Si prende questa equazione e la si impone, io quindi sto facendo quel problema lì, la impongo al tempo tn e quindi avrò m per un due punti più k per un uguale fn fn è noto, un è noto, un due punti è noto ok, allora dice ma che mi serve? mi serve per calcolare la stima al basso successivo utilizzando questa formula le differenze finite centrali. Cioè io approssimo l'accelerazione un n due punti con questa formula che probabilmente avete già visto e che è accurata a seconda ordine, richiede di conoscere, è espressa in funzione di un n più 1, quindi lo spostamento al passo successivo, un lo spostamento al tempo attuale che è noto e un n meno 1 che a maggior ragione è noto perché è già stato calcolato nei passi precedenti. C'è un primo problema, è che se vogliamo partire al tempo 0, cioè vogliamo applicare questa equazione al tempo 0 per calcolare la prima stima di un n più 1, non conosciamo, non sappiamo neanche cosa vuol dire un n meno 1. Beh, allora un n meno 1 si attiene da, se io conosco la velocità al tempo zero, allora la posso immaginare come un zero meno un n meno 1 diviso delta t da questa relazione. Questa lo conosco sulle condizioni iniziali sullollo spostamento queste sono le condizioni iniziali sulla velocità quindi riesco a calcolare la stima di questo un meno 1 da lì parte la procedura e poi continua in maniera autonoma allora quindi se io prendo questa espressione la sostituisco qua dentro va bene, con qualche piccolo passaggio riesco a ottenere un più 1 in funzione di questo right inside in cui c'è fn meno un contributo della rigidenza al tempo n più 1 e poi questi contributi che vengono da questo oggetto basta sostituire, cambiamo un pochettino che cosa devo saper fare? devo sapere risolvere un sistema con la matrice M-1 quindi risolvere un sistema in cui la matrice è M perché questo è il termine noto quindi M-1 per questo termine noto è la risoluzione di un sistema e poi aggiungo questo contributo ora, così come abbiamo fatto nell'esempio sulla diffusione, molto spesso si è visto che l'accuratezza di questo processo rimane ottima se invece di usare la matrice M calcolata come abbiamo visto prima, prendiamo un'approssimazione diagonale. Se M fosse diagonale, immaginate che veramente M sia diagonale, la nostra soluzione è immediata, è banalissima, e non richiede nemmeno la soluzione di un sistema, se M è diagonale, per questo si chiama esplicito. Abbiamo la soluzione in forma esplicita senza neanche dover risolvere un sistema. Questo processo di ridurre la matrice di M a diagonale si chiama processo di lumping della matrice. Lumping. E la matrice a diagonale si dice la matrice lumped. Come si fa? Ma ci sono tantissime tecniche. La matrice di massa in generale avrà dei coefficienti diversi da zero e però è sparsa. La cosa più semplice, si fa la somma di tutti gli elementi di una riga e si attribuisce il valore all'elemento diagonale e così via si procede e si trasforma in matrice diagonale. È possibile mostrare che così facendo si preserva la massa della struttura. Se applichiamo questa semplicissima procedura a elementi quadratici, troviamo dei coefficienti negativi sulla diagonale e non va assolutamente perché la matrice di massa è definita positiva e vediamo che rimanga tale. Una matrice di massa diagonale con elementi negativi non è più definita positiva. Allora ci sono tecniche alternative, non è il mio obiettivo spiegarvi quali sono, ma sappiate che tecnica corrente quando si usa questo metodo di differenze centrali che è esplicito, usare le tecniche di lumping della massa. Se si usa questo, ogni time step è velocissimo. Il vantaggio di questa tecnica è che è strarapida. E l'altro vantaggio è che per risolvere ogni time step, Fn, questo sarà importante per il nostro problema della ruota che imbalza Fn dipende solo dalla soluzione nota al tempo tn questo è importante perché è valutato al tempo tn il P di cn vuol dire che questa quantità è valutata a cn quindi se abbiamo una forza come la forza di contatto che dipende dalloostamento, Fn dipende solo dallo spostamento al tempo tn, quindi è noto, quindi può essere calcolato. Mentre se dipendesse dallo spostamento finale, allora sarebbe un problema fortemente non lineare, perché la forza di contatto dipenderebbe dalla soluzione ancora incognita, e quindi sarebbe un bel pastic pasticcio quindi quando ci sono problemi fortemente non lineari, algoritmi strasemplici come questo hanno il vantaggio di richiedere la valutazione delle forze solo al tempo n quando lo spostamento è noto ed è per questo che utilizzeremo il codice che ho implementato questa facilissima versione allora sappiate lo dico però solo a livello di curiosità, che però altri algoritmi sono preferiti. Il tipo di algoritmo utilizzato in meccanica dei solidi è quello della famiglia di Newmark, che si basa su queste due ipotesi. Guardiamo un po' cosa c'è scritto. Partiamo da queste sotto. Noi sappiamo che in generale si può fare un'espansione in serie di tela, ma se io voglio ottenere un più uno, posso partire da un, aggiungere delta t per la velocità, poi un mezzo delta t quadro dell'accelerazione e commetto un errore dell'ordine di delta t quadro e potrei andare avanti e lo stesso lo posso fare con la velocità. Quindi queste sono delle espansioni in serie di Taylor. Il metodo di Neumark applica queste espansioni, infatti vedete che i primi termini sono uguali, questi qui sono gli stessi, questo è lo stesso, solo che modifica questo termine qua. Invece di usare delta T4 mezzi per l'accelerazione al tempo n, la modifica con un mix dell'accelerazione nota più l'accelerazione incognita al tempo Tn più 1. E il mix è calibrato da questo parametro, il parametro beta. Notate che facendo beta è compreso fra z e un mezzo, se beta vale 0 otteniamo esattamente questo. Se beta vale un mezzo, allora questo scompare e tutto dipende solo dall'accelerazione alla fine del passo. Quindi sostanzialmente Newmark fa una sorta di espansione in serie di Taylor sostituendo all'accelerazione al tempo tn una media pesata dell'accelerazione al tempo tn al tempo tn+. E' la stessa cosa per la velocità. Invece di avere l'accelerazione nell'espansione di Taylor della velocità, invece di avere l'accelerazione al tempo tn, fa una media pesata fra le due accelerazioni ottenendo in funzione di un parametro diverso, gamma, gamma in questo caso invece è compreso fra 0 e 1, se gamma è 0 ottengo esattamente questa, se gamma è 1 invece ottengo l'accelerazione solo al tempo successivo. Questa è un'ipotesi, se prendete questa equazione e la imponete al tempo tn più 1, cioè prendete all'interno di questa espressione, mettiamo lo spostamento, la velocità in questa forma non c'è, se prendiamo al posto di questo, mettiamo questa cosa qui, otteniamo dopo un po' di lavoro qualche passaggio al genico e la soluzione di questo sistema. Otteniamo una equazione per stimare l'accelerazione alla fine del tempo, alla fine del passo temporale, che è funzione, però, e questo è il problema che accennavo prima, è funzione di fn più 1, cioè delle forze alla fine del passo. Se per caso queste dipendono dallo spostamento, allora sono dei pasticci che devo iterare. E poi queste sono quantità invece note, perché sono tutte espresse al tempo tn. Per cui il passo, la transizione dal tempo tn al tempo tn più 1, si fa in questa maniera, calcolando l'accelerazione al tempo finale, utilizzando questa relazione, e poi ritornando in queste relazioni per calcolare lo spostamento e la velocità. Una volta che si è calcolata la nuova accelerazione, tutto il resto è noto, si calcolano lo spostamento e la velocità allaa del tempo. E così l'algoritmo va avanti all'infinito. La grande differenza rispetto all'algoritmo di prima è che vedete che dobbiamo risolvere un sistema in cui la matrice è una combinazione di m e di k. Se m si comporta bene anche se diagonalizzata, con k è assolutamente impossibile e quindi in generale questo è un algoritmo implicito si dice algoritmo implicito perché necessita sempre della soluzione di un sistema lineare a ogni passo temporale e va bene, questo poi c'è un riguardo ma non vogliolio entrare all'interno di questo, perché come si possa utilizzare quella relazione per procedere passo dopo passo è assolutamente ovvio, siccome poi non analizzeremo un codice fatto con questa tecnica, è del tutto inutile analizzarlo. Volevo solo dire due parole su un paio di argomenti che però sono fondamentali, perché ne abbiamo già visto l'impatto drammatico sulla diffusione, la stabilità e la instabilità, esplosione della soluzione e non esplosione della soluzione. C'è un famoso teorema che vale nei problemi lineari, però vedremo che poi la sua validità rimane abbastanza vera anche nei problemi non lineari, è un teorema che si chiama teorema di Lax che garantisce che la soluzione che si ottiene utilizzando l'alberino di Neumark converge alla soluzione esatta delle equazioni differenziali, cioè per me il problema è questo, il problema di riferimento è questo qui. Mi dimentico il fatto che ho già discretizzato uno spazio, che sto già facendo degli errori a causa della discretizzazione dello spazio. Io voglio risolvere questo problema che per me è il problema esatto, del resto non mi interessa più. Quindi questo teorema dice che la soluzione dell'algoritmo di Neumark tende alla soluzione esatta di quel sistema se l'algoritmo è stabile e consistente. Di queste due cose trovate dimostrazione all'interno delle dispense, sia la stabilità sia la consistenza. La consistenza è un concetto che ci interessa poco, garantito automaticamente. E' molto più interessante la stabilità, perché impatta drammaticamente sul risultato dell'analisi. Cosa vuol dire stabilità? Si può lavorare un pochettino sul sistema e trasformare il nostro sistema in questa forma. Cioè, che cosa succede qua? Guardiamo questa relazione qua sotto, poi i passaggi, se volete, ve li guardate. È possibile esprimere un n più 1 e un punto n più 1, cioè lo spostamento della velocità alla fine del passo, come il prodotto di una matrice che moltiplica la stessa cosa, lo stesso vettore di stato, spostamento di velocità, al tempo n. Quindi questa Q più un termine che dipende dalla quantità note al tempo tn. Questa matrice Q viene chiamata matrice di transizione, cioè noto il vettore di stato al tempo n, fa la transizione al vettore di stato al tempo tn più 1. Uno schema time stepping è stabile o instabile a seconda che questa Q amplifichi o riduca un eventuale errore presente in questo vettore di stato. Cioè, supponete che questo vettore di stato contenga un errore. Questo errore viene amplificato o viene ridotto dalla Q? Magari nella slide dopo c'è il... Ecco, qui, questo. Se io ho un errore sulle condizioni iniziali, sullo spostamento e sulla velocità, questa Q amplificherà eventualmente questo errore o tende a snussarlo? E vabbè, allora c'è tutta un'analisi particolare che passa attraverso gli autovettori e gli autovalori di quel sistema lungo, sono una decina di pagine e sarebbe impossibile da far qua, l'analisi quindi di questa matrice di transizione si fa sulle componenti modali di spostamento e velocità per disaccoppiare tutte le equazioni un'un dall'altra e alla fine quello che mi interessa è solo questo messaggio. A seconda dei valori di gamma e di beta, l'alcorismo può essere incondizionatamente stabile, condizionatamente stabile o instabile. Che vuol dire? Se gamma è minore di un mezzo, il schema di Newmark è instabile, cioè se c'è una piccola perturbazioneazione questa viene amplificata e la soluzione numerica letteralmente esplode. Vi ricordate il caso della sfera, lo vedremo dopo nella nostra soluzione. Se gamma, quindi gamma minore di un mezzo, ovviamente è inutilizzabile, ma se gamma è maggiore di un mezzo, a seconda del valore di beta, l'algoritmo può essere condizionatamente stabile. Che cosa vuol dire? Che il delta t che noi utilizziamo per avanzare nell'analisi è limitato superiormente da un valore che dipende da questa quantità. Cioè noi possiamo avanzare con un delta t a condizione che questo sia sufficientemente piccolo. Questi cosa sono? Gli omega j sono gli autovalori del sistema. Quindi se io devo... e gamma e beta sono i coefficienti utilizzati. Quindi se io devo prendere il minimo di questa quantità è ovvio che devo andare a prendere l'omega massimo. Quindi l'omega massimo, il massimo autovalore della mia analisi, limita il massimo delta t che devo usare e dovete sapere che in genere l'omega è legata alla taglia dell'elemento più l'elemento è piccolo più il massimo omega è grande più raffino la mesh più il massimo omega diventa grande quindi questo vincolo potrebbe essere veramente stringente se stiamo lavorando con una mesh fine in cui il massimo omega è veramente una quantità altissima. Comunque vale questo ragionamento. Se invece utilizziamo un beta tale che questa quantità sia maggiore di zero, lo schema di numeri è che è incondizionatamente stable, cioè possiamo procedere con il delta t più gamma. Sarà più o meno accurato, ma non esplode. E quindi tipicamente si utilizza gamma uguale a un mezzo e beta uguale a un quarto. Ci si mette lì vicino all'interistabilità, gamma uguale a un mezzo e beta uguale a un quarto è una scelta molto praticata. Vale la pena di dire che se prendete beta uguale a zero e gamma uguale a zero ottenete esattamente lo schema delle differenze centrali di prima. Quindi se entriamo beta uguale a gamma uguale a un mezzo è esattamente il limite inferiore di validità per gamma. E che succede? Che se beta uguale a zero, allora siamo in questa situazione. Quindi differenze centrali è condizionatamente stabile. E possiamo avere una stima analitica di stabilità, chiediamo a MATLAB di calcolarsi l'autovalore più grande, facciamo due diviso e otteniamo il delta T massimo con cui possiamo procedere. L'idea è che abbiamo capito un pochettino come succede, cosa succede, andare nei dettagli richiederebbe un sacco di lezioni dedicate, se siete curiosi vedete la procedura, a me interessava semplicemente darvi un po' l'idea di cosa sia la stabilità, la consistenza neanche velocità vediamo praticamente come può impattare sul nostro problema allora abbiamo questo pneumatico qua che tutto è codificato all'interno del codice Dynamics provate a lanciare Dynamics si dovrebbe arrabbiare a un certo punto dovrebbe darvi un errore di questo genere cominciate comunque a lanciare comunque l'abbiamo trovato andiamo a vedere, andiamo a analizzare. Perché questo è un piccolo esercizio, adesso vediamo se abbiamo tempo di farlo, se non abbiamo tempo vi dico come risolvere subito il problema. Che cosa fa questo codice Dynamics? Beh, se lo scrollate vedete esattamente tutta la parte iniziale che è identica, tal quale a modal. L'unica differenza la trovate per la prima volta alla linea 99. La linea 99 è quella routine che noi dovremmo completare e che tiene conto delle forze di volume, perché la gomma cade rispetto al peso e noi il peso non l'abbiamo mai messo negli analisi, quindi dobbiamo essere capaci di introdurre il peso. Se riusciamo a assemblare la K, la M e la F, che sono i tre ingredienti di quel sistema che abbiamo visto prima, poi l'algoritmo procede passo a passo con l'implementazione delle differenze centrali. Allora, vediamo prima di tutto questa cosa. Allora, un commento su questo esercizio, veramente l'esercizio dovrebbe richiedere 5 minuti massimo, vediamo, se in 5 minuti ce la facciamo bene, se no vi do la soluzione. Richiede di calcolare questo termine qua, il peso, la gomma è inizialmente separata dalla superficie rigida e comincia a cadere, comincia a cadere il perfetto del peso. Dobbiamo introdurre il peso proprio, cioè dobbiamo calcolare, abbiamo sempre detto che anche quei termini avrebbero generato un contributo alle forze nodali che abbiamo chiamato FE, però non l'abbiamo mai calcolato effettivamente. Diamo un occhio solo a come si può, esempio prendere in considerazione il peso cioè nel caso del peso andiamo a vedere nel file di input apriamo il file tire al dynamics al primitutto dynamics-tire.m. C'è la mesh che è generata nell'opportuno file, il materiale contiene tre informazioni, yang, poisson e densità, solide informazioni sul solido e il materiale, però devo cominciare a dare anche le condizioni iniziali. Impongo che la velocità iniziale sia verso il basso V0 uguale a meno 0,2. Non ho segnato lo spostamento iniziale, perché lo spostamento è dove ho generato la mesh. E poi introduco queste body forces BF che sono 0 e meno 3, cioè una forza di volume che non ha componenti in direzione X ma ha una componente in direzione verticale, perché è meno 0,3 totalmente fittizio. Quindi una forza di volume negativa. Queste altre quantità sono quantità che vedremo dopo, sono associate al contatto, adesso non ci interessa. Come possiamo introdurre queste body force? Quando MATLAB si legge il file input, trova quanto vale questo vettore, il vettore BF che è associato alle forze di volume, che sarà una lista di due componenti in cui la prima è 0 e la seconda è meno 0,3. Per cui dobbiamo essere capaci di integrare quella quantità lì. Per cui, come sempre, si fa l'integrale elemento per elemento, quindi dobbiamo essere capaci di effettuare questo integrale su di un elemento, e poi si assembla, e la procedura di assemblaggio è identica a quella della matrice di rigidezza e di massa ecco F è un spieghiamo solo questa espressione poi andiamo a vedere cosa dobbiamo scrivere all'interno della routine perché se capiamo questo come modificare la routine anzi come aggiungere come mettere le due linee che mancano nella routine è banale. Allora, noi vogliamo integrare questa cosa qua. Adesso probabilmente dovreste riuscire a intuire anche voi, perché facciamo un scambio di l'etat, ormai è una cosa abbastanza ripetitiva, comincia a diventare abbastanza ripetitiva, si fa sempre alla stessa maniera. Devo integrare sui elementi questo oggetto. Allora, innanzitutto vi ricordo che invece di usare il vettore preferisco usare la lista n, la lista W, che sarà W per n per W, la famosa matrice n delle funzioni di forma. Poi invece questa qui la metto all'interno di una lista F che contiene le componenti di questo vettore, che quindi erano in mio caso 0 e meno 0,3. È quello che il matravaletto invierte, cioè mettiamo queste due componenti. All quindi questo è noto quindi se inserisco lì dentro avrò WE trasposto per l'integrale sull'elemento di che cosa? DN trasposto per F ma F è noto, per cui devo semplicemente effettuare questo integrale questo qui lo chiamo FE questo è il vettore di forza nodale quindi devo essere capace di effettuare questo integrale come si fa? come indicato qui sempre la stessa maniera mappa dall'elemento fisico all'elemento paramedio quindi apparizione dello jacobiano e l'unica cosa che dovrò inserire saranno queste funzioni di forma andate a vedere andiamo a vedere perché ci ha dato errore il simpaticone di codice andiamo aprire c'è dato errore perché non riusciva a capire bene che cosa combinare in T6FE. Ha scritto questo, ha detto T6FE line 26 colonna 3, c'è qualcosa che non va. Andiamo a aprirlo. T6FE, eccolo qua. Sì, certo che dirà perre, perché io ho fatto uno scherzo e ho tolto due linee. Quindi lui trova i puntini e dice no, non va. Allora, provate a spaccarvi la testa tre minuti e a dire, anche per somiglianza con le altre cose che abbiamo visto sinora, come dovete riempire, io vi metto le due slide vicino alla file matrim, il file matrim ce l'avete voi, è inutile che lo lasci, me lo metto qui di fiango, provate a spaccarvi la testa tre minuti per stabilire che cosa dovete scrivere, manca una riga qua, una riga e dovete completare quest'altra riga. Provate a pensarci 5 minuti, veramente, non più, non più di 5 minuti perché sono già alle 5.05, quindi diciamo alle 5 e alle 12 di do la soluzione. Provate a pensarci se volete fare domande io sono me La cosa divertente è che se riempite bene con le due linee il codice va e una a sei, è una lista, ispiratevi alla matrice di massa che abbiamo visto prima. Un'altra cosa del BF è proprio questo vettore qua, questo rozzero per questa lista si chiama BF, viene dato in ingresso, viene dato in ingresso alla routine, questo BF. Quindi quello con cui avete bisogno sono J e N trasposto, sostanzialmente. Grazie. A di calo come si... è di Gauss però di G a Gauss Bene. per l'inverso Ok. Vediamo se... E questo... E' quello che stavamo facendo girare? Dynamics! Ok, grazie! NL trasposto per... cosa fa di per l'inversa di F? andiamo assieme allora prima di tutto abbiamo bisogno di creare questa N che non è NL perché NL è una lista N invece è quella matrice che oggi abbiamo visto tante volte la matrice 2 per 12 dove l'abbiamo vista? la matrice 2 per 12 all'interno della matrice di massa quindi apriamo la funzione ah ecco perché mi dava errore pensavo fosse già era da ricreare quindi andiamo a vedere la matrice elementare di massa che fa questa cosa qui e copiamo ricordatevi la stragrande maggioranza delle cose viene copiata poi dopo apriamo ME vedete che qui segue la stessa cosa prima definisce NL la lista e poi si crea questo N questo N riempiendo le posizioni dispari le posizioni pari la matrice N è sempre lei adesso mi copio anche quello che fa l'integrale. Vado qui, vediamo se funziona. Copio questo. Quindi definisco n. Poi ho portato anche un'altra cosa che mi serve da copiare, però per il momento non la commuto. Quindi per prima cosa ho copiato la definizione della matrice N. Ho copiato la matrice N, però dopo devo fare l'integrale. Vedete che l'integrale FE richiede di fare N trasposto, cominciamo a compilarlo. N trasposto è N con l'apice. Poi, che moltiplica BF, la forza di volume, che moltiplica lo jacobiano J ma lo jacobiano J ce l'ho già perché l'ho calcolato qua sopra che moltiplica il peso di Gauss lo copio dalla cosa sopra J per Gauss cancello questa riga che mi ero cambiato vediamo se funziona e a volte bisogna vedere se bf è una riga o una colonna non mi ricordo mai io provo semplicemente a farmi andare e no vedete allora non va bene perché le dimensioni non sono giuste, allora faccio il debugging, decido di entrare in FE e mi metto al posto giusto qua, lo bloccherò a questa linea in modo da vedere che cosa diavolo succede, faccio partire Dynamics e lei si ferma proprio a quella cosa lì, allora andiamo a vedere, Bf è 1,2, vedete è una linea, è una riga, invece ho bisogno di una colonna, beh allora semplicemente bisogna fare Bf trasposto. Ho sbagliato, avrei dovuto ricordarmelo, abbandono del buggy e faccio Bf, semplicemente Bf trasposto con BF appice vediamo se basta vediamo parte, fa l'integrale vediamo cosa succede ecco vedete qui fa apparire il pneumatico poi per farlo andare questa è una time history quindi bisogna schiacciare questo pulsantino per far partire la time history questo comincia ad andare giù e scompare diavolo ma non doveva fermarsi sulla superficie allora adesso ci manca da capire il contatto comunque funzunque funziona, cioè accelera verso il basso. A me rimbalza già. Adesso vediamo. A te rimbalza già? Sì. Forse c'è già il file. Ecco perché, perché io oggi ho modificato il parametro. Allora, fate così. Andate su Dynamics Tire. Vedete questo lambda. Voi evidentemente avete lambda uguale 10, io l'ho messo a 0. Se mettete lambda a 0, fate così, mettete lambda a 0 e fate a 0. Dopo vi spiego cos'è questo lambda. Volevo prima partire da questo. Se mettete lambda a 0 dovrebbe succedere la stessa identica cosa che succede con me d'accordo? ok, allora il primo commento la prima cosa che mi interessa è il delta t stabile il codice, vedete qui si è calcolato con la formuletta che abbiamo visto prima questo è un algoritmo di differenze città, si è calcolato il massimo delta t stabile, 0,0033. Allora qui voi potete, qui dentro noi definiamo i parametri dell'analisi dinamica, il tempo finale, cioè arriva avanti fino al tempo 10, integrare, e poi decide qual è il delta t da utilizzare. Vedete che io ho usato un delta t che è leggermente più piccolo del delta t stabile? Allora us che è già molto più vicino anzi facciamo così non andiamo per i grandi se no non finiamo più utilizziamo un delta t che è leggermente superiore a quello del 0,34 0,34 e rifate andare l'analisi quindi è 0,0034 quindi veramente un epsilon più grande della dovrebbe succedere l'irreparabile già lo vediamo perché è da dei nan ecco questa è l'immagine di un'esplosione quindi è successo certamente qualcosa che non va bene Niente, si rifiuta persino di farmela vedere. Allora, facciamo così, lo inganniamo un po', usiamo un delta T, questo era il motivo per cui prima era così, ecco, usiamo questo delta T che è vicinissimo, un po' più piccolo, lui dovrebbe rimanere stabile perché è leggermente più basso. Però vedete che già si stabilizza. Poi esplode. E poi a un certo punto comincia ad andare giù e poi esplode. Questa è l'esplosione. Cioè se ci sono delle piccole perturbazioni nell'analisi, ci sono sempre, l'algoritmo se è condizionatamente stabile le fa esplodereiamo a un delta T anche marginalmente superiore a quello di stabilità questo è il primo concetto e poi andiamo a vedere mettiamo un delta T più basso quindi rimettiamo il delta T 0,02 che stavo utilizzando prima cioè ci teniamo a distanza di sicurezza dal limite di stabilità e andiamo a vedere cosa succede col contatto. Perché manca il contatto, chiaramente. Allora, come facciamo a simulare il contatto? Simulare il problema è una maniera semplificata di simulare il contatto, che ha una maniera esplicita, anche questa, che funziona benissimo insieme alle differenze centrali. Queste sono le relazioni che governano il contatto, sono le relazioni seguenti, un problema fortemente non lineare che dice la seguente cosa calcoliamoci la y2 che è la coordinata corrente di un nodo che è data da x2 più 2 che queste sono x2 è la coordinata iniziale di un nodo e 2 è lo spostamento, cioè calcoliamo la coordinata y, cioè la coordinata attuale del nodo. Allora y2 è governato dalla condizione che y2 deve essere un maggior uguale di 0, ma se la superficie rigida è in 0 c'è questo vincolo, c'è questo vincolo. Poi c'è un'altra condizione sulla forza, la forza assegnata in direzione 2, la forza che si è citata dalla superficie rigida sulla gomma, può essere solo positiva. Non ci può essere una forza verso il basso. Siamo d'accordo? Al limite può essere o zero o positiva. Però le due cose non sono indipendenti tra di loro, perché la forza si può attivare solamente se y2 è uguale a zero. se invece y2 è maggiore di 0 questa è 0 questa alternativa si esprime in questa maniera y2 per f2 è uguale a 0 sempre il prodotto tra queste due cose è sempre 0 questo è un problema stranoso che governa mille problemi nella fisica e matematica è un di complementarietà, è un problema che si chiama problema di complementarietà, un problema estremamente, fortemente non lineare, difficile da risolvere. Allora in maniera molto pragmatica noi sostituiamo la superficie rigida, perché poi le superficie rigide in realtà non esistono, con una superficie simulata con una specie di letto di molle, per cui abbiamo che la forza, scusate la lettera, qui T, lì ho usato F, T è data da lambda per delta, dove lambda è un coefficiente, il coefficiente penalty, quello che trovate definito nell'input, e rappresenta la rigidezza delle molle. Delta è la penetrazione della superficie della gomma all'interno della superficie, perché x2 più 2, cioè x2 più 2 è la coordinata corrente. Se questa diventa negativa, allora applico questa forza. Quando c'è compenetrazione, quindi controllo il valore di questo delta. Se questo delta diventa negativo, allora applico, anzi, se meno delta diventa negativo, applico una forza verso l'alto, che è proporzionale alla compenetrazione. Quindi sto simulando questa non più come superficie rigida, ma come un letto di molle. Come faccio a... questo algoritmo è tutto... perché funziona bene con le differenze centrali? Per il motivo che vi ho detto prima, che l'ho già commentato ma lo ripeto, a ogni passo temporale io ho bisogno di calcolare le forze Fn che sono note al tempo Tn. Fn è funzione dello spostamento al tempo n, ma se lo conosco posso farlo. Quindi posso calcolarlo al tempo n. Calcolare la stima al tempo finale, la stima dello spostamento al tempo finale in funzione delle forze esercitate al tempo Tn. Quindi a ogni istante temporale n vado a vedere se c'è compenetrazione, se c'è compenetrazione applico una forza verso l'alto e quindi questo simula la reazione del solido. E questo, vabbè, non vi teglio con l'ennesimo integrale, perché questo comunque sarebbe la prima volta che troviamo un integrale di superficie di questo genere, e quindi spieghiamo un pochettino che questo si calcola esattamente con le stesse tecniche, cioè lunica differenza è che questa è un integrale di superficie cioè l'integrale sul bordo di un elemento un integrale sul bordo non è un integrale di superficie però si fa esattamente come il resto vediamo che così dopo vi lascio e qui trovate se siete curiosi entrate in questa routine di 3f contact che fa l'integrale ma mi premeva di più farvi vedere come funziona andiamo a mettere l'anda diverso da 0 fatelo anche voi così vedete la differenza mettiamo l'anda 10 quindi attiviamo il contatto per che a me non rimbalzava, perché gli avevo messo rigidezza delle molle 0, quindi non applicava forze. E quindi gli facciamo andare la soluzione, voi l'avete già vista. Facciamo andare la time history. E' evidente che lui alla fine rimbalza, però penetra tanto la superficie rigida. E come si fa per farlo penetrare di meno? Noi abbiamo un solo strumento per farlo. Dobbiamo aumentare lambda. Dobbiamo aumentare la rigidezza. Quindiiamo lambda a 100 e vediamo cosa succede. Questo si chiama coefficiente di penalty del contatto. Aumentando a 100, poi sono tutti parametri che devono essere calibrati, noi li stiamo fissando un po' così a caso. Vedete che la condizione di contatto è rappresentata in maniera molto migliore, insomma il rimbalzo della gomma è quasi realistico. Chiaramente più è rigida la superficie, più il contatto è impulsivo e voi sapete che un impulso eccita tutte le frequenze della struttura. Vedete che infatti la risposta della gomma diventa sempre più nervosa. Diventa molto nervosa perché eccita tutte le frequenze, anche quelle più elevate. Allora, se tu devi andare, per carità vai. Io ne continuo solo 5 minuti per fare vedere una cosa curiosa. La domanda che vi lascio alla vostra riflessione è Matteo ma perché funziona? non stiamo facendo tutto per degli spostamenti infinitesimi non sono le ipotesi del capitolo 2 non è una trasformazione infinitesime è una cosa che lascio alla vostra riflessione perché funziona? non dovrebbe funzionare dovrebbe succedere qualcosa di drammatico allora per portarvi sulla strada giusta cambio il file di input invece di prendere la ruota con superficie diritta prendo una ruota che rimbalza su una superficie inclinata e ci metto anche l'attrito il fatto sta che quando la ruota cade riceve una forza anche tangenziale in direzione opposta alla velocità che la mette in mutazione come un vero pneumatico vediamo se le cose sono calibrate bene vediamo cosa succede poi potete andare a vedere negli input come sono definito l'inclinazione, in questo caso mi sembra che sia pigreco sesti. Qui vedete la gomma, no, qui mi sembra più pigreco quarti, comunque bisogna andare a vedere. E questa è la superficie rigida su cui rimbalza. Andiamo a vedere che cosa succede. Lui cade. No! Come è possibile? Scusate, no. Allora è solo il disegno. Devo andare a cambiare. Scusate, devo cambiare. Evidentemente ho messo il parametro. Oppure devo cambiare qua sotto perché l'ho fatto un po' di giorni fa. Ah sì, devo cambiare la routine, che tiene conto, scusate, che tiene conto del fatto che la sommetica è inchinata, dovrebbe essere, credo solo questo, vedo solo se le parole sono giusti, sì, vedete che c'è angolo meno P6, quindi evidentemente la visualizzazione non è perfetta, ma c'è un attrito con efficiente 0,3 quindi non avevo solo cambiato la routine elementare di calcolo cioè non metteva l'attrito e metteva solo una forza verticale per cui lui ribalzava esattamente come prima ok adesso dovrebbe funzionare tac, ved vedete, lui schiaccia così, cosa sui? Beh, che succede qui? Una cosa assolutamente non fisica, questo cerca di ruotare, però, io non ho mai visto una gomma a rimbalzare diventare grande così, c'è un problema, che è legato anche alla domanda che mi facevo prima. Perché prima funziona e perché quando funziona? La risposta è la stessa. Provate a darvi voi la risposta, però vi do un suggerimento. Tutto è legato al gradiente di U. Il gradiente di U deve essere piccolo. Deve essere piccolo perché tutto valga. Tutte le nostre cose le stiamo facendo sui gradienti di U. Piccolo. Perché in un caso va bene e nell'altro non va bene? Perché non è verticale, quindi subendo spostamenti... Il gradiente di U è sull'asse che cos'è che viene filtrato quando calcolo un gradiente? qui c'è pure l'angolo e fa no ma cominciamo a dire perché funziona nel primo caso al limite se fosse rigido U scenderebbe tutto sarebbe costante rispetto allo spazio avremmo uno spostamento costante una traslazione quindi quello non è violato qui viene pure accelerato qui invece viene ruotato e la condizione per che quello sia vero è che la rotazione sia infinitesima ma qui è alto che infinitesima qui c'è che il graio è più gira infatti e più gira più sbaglia quindi quello che stiamo facendo è valido solamente se la rotazione è infinitesima ma qui la gomma cerca di fare un giro completo e quindi a un certo punto ci dice cose assurde perché stiamo usando una formulazione che non non va bene. Per simularlo bene, bisognerebbe utilizzare il famoso tensore di Green Lagrange per le deformazioni, che ho la Kirchhoff di stress, ma diventa tutto molto più complicato e soprattutto diventa un problema fortemente non lineare. Magari l'ultimissimo giorno vi dirò qualcosa, perché questo è straimportante per i MEMS, soprattutto i MEMS moderni, per riuscire a tenere in conto i grandi spostamenti. Però è veramente complicato in un corso di introduzione esagerato. Va bene, dai, finiamo qua.