att4.txt

 però vediamo in questa settimana prossima non c'è la lezione perché il primo novembre approfitto magari di questa pausa per pensare di più un po' di più all'esame finale allora l'argomento di oggi è duplice alla fine arriveremo a trattare una cosa che per i MEMS è molto importante e che è la cosiddetta dissipazione termoelastica che cosa significa e in che casi è effettivamente molto importante. Vedremo che per parlare di dissipazione termoelastica è necessario innanzitutto ricordarci il concetto di modi meccanici che abbiamo introdotto l'altra volta. Il modo meccanico è un elemento base per il funzionamento di tutti i microsistemi, quindi sapere da dove viene, come si calcola è molto importante. Poi voi non implementerete mai un codice che calcola questi autovettori, autovalori, però dentro ai codici che userete, che siano ANSI, se sarete in ST, se è CONS, se starete a lavorare un po' in università o comunque altri codici, però insomma è fondamentale sapere che cosa sono, che proprietà hanno eccetera eccetera poi abbiamo bisogno della diffusione del calore che è la prima parte della lezione però la diffusione del calore la faremo molto è il capitolo 5 attualmente delle dispense, però lo faremo molto rapidamente perché in fondo era il proprio primo argomento che avevamo visto la prima lezione per cui è vero che in quel caso eravamo limitati a un problema 1D perché avevamo assunto simmetria sferica quindi avevamo lavorato in coordinate sferiche e tutto dipendeva solo dal raggio però la sola differenza rispetto a quello che faremo oggi quindi molto rapidamente con un drive e qualche simulazione da vedere sarà la prima parte, la forma debole che poi fra l'altro avevamo sviluppato come prima cosa nella lezione allora quindi prima ricordarci un pochettino della diffusione del calore poi una prima applicazione è la termoelasticità non accopiata nel senso che se c'è una variazione di temperatura all'interno di una struttura questo induce degli sforzi se la struttura non è completamente libera di espandersi e ci sono quindi gli sforzi di origine termica di cui non abbiamo ancora parlato ma che spesso sono importanti nelle microstrutture soprattutto per quanto riguarda il package le strutture sono attaccate sono ancorate al package il package molto spesso è di natura plastica o simile va a finire all'interno del cellulare o vicino al motore di una macchina si riscalda quindi questo si deforma e altera completamente il comportamento del vex grave grave problema motivo per cui la configurazione ideale del mens è quella fungo in cui c'è un solo ancoraggio c'è un singolo ancoraggio e tutta la struttura è sospesa questo ancoraggio per cui anche se il pacchetto impegno ci si deforma L'effetto sul comportamento del mens è minimo Soprattutto agli inizi è stato un problema enorme questo delle deformazioni termiche del ha un suo senso però non è l'obiettivo di oggi perché dopo aver visto queste cose passeremo a utilizzare tutti questi concetti nel capitolo 6 sulla dissipazione un tema enorme importante oggi vedremo la dissipazione termoelastica poi magari volte prossima avremo dei seminari sulla dissipazione fluida e ne vi racconterò anche io direttamente alcune cose sulla dissipazione fluida. Quindi questa è la parte più importante della lezione, la parte più nuova perché il resto è nuovo per modo di dire, c'è un po' ormai in edizione di cose che dovrebbero essere abbastanza vediamo cominciamo con dal cominciamo in ordine andiamo in ordine e quindi cominciamo con la diffusione ricordiamo alcuni concetti già visto il che avevamo già visto il la prima lezione ok che succede? va bene ho fatto una piccola modifica ho visto che mi serviva soprattutto per trattare i problemi accoppiati se vi ricordate il se vi ricordate il campo test in elasticità la velocità testa la indicavo con W che era una lettera diversa dalla velocità la velocità era lo spostamento in realtà era U e il campo testa era W siccome poi ad esempio quando parlo della temperatura ho bisogno del campo della temperatura, del campo testa di temperatura, alla fine si usano 100.000 lettere diverse per uniformare in questo contesto almeno il campo testa in unico con lo stesso simbolo che uso per la variabile principale con una tilde sopra ok? così poi dopo pian pianino trasformerò anche il resto con la stessa notazione è solo un cambiamento di presentazione non cade assolutamente niente nella nella nella procedura ok allora l'equazione forma forte del problema forma forte del problema è un problema di evoluzione in cui però a differenza della dinamica strutturale abbiamo la derivata prima solo della temperatura, ρ0 è la densità, c è la capacità termica, quindi questo ci dice che ρ0 per c per la derivata della temperatura rispetto al tempo è uguale a questo termine, leggiamolo, tra parentesi c'è k per il gradiente di t, ma K per il gradiente di T a parte il segno è ciò che origina il flusso termico infatti se noi prendiamo K per il gradiente di T e facciamo un protoscolo normale al contorno questo deve essere proprio uguale al flusso termico quindi qua sopra a parte il segno questa è la divergenza del flusso termico in elasticità lì c'era la divergenza di sigma qui invece c'è la divergenza del flusso Q e questo deve essere uguale a una sorgente per generalità ho lasciato questa sorgente poi vedremo che nella seconda parte della lezione questa sorgente assumerà un significato molto preciso e sarà legata alla deformazione del solido è un po' come il gas che quando si espande che succede? d'estate buttiamo acqua che si evapora e si raffredda perché assorbe calore e quindi l'espansione induce un assorbimento di calore e una contrazione quando invece il gas si condensa cede e quindi a seconda della contrazione o della dilatazione anche nel solido questo avviene, un solido che si dilata è una specie di pozzo di calore e un solido che si contrae invece è una sorgente di calore in genere questo effetto è molto piccolo e lo si trascura, invece nei nostri microsistemi diventa un effetto molto importante però per adesso questa è la diffusione classica in cui con G si intende una classica sorgente di calore, qualcuno che va lì e scalda allora questa equazione deve essere valida per ogni punto del dominio e per ogni tempo in un intervallo che si vuole analizzare tipicamente parte da zero e arriva a tf tempo finale dell'analisi ok questa quindi è la condizione diciamo l'equazione in forma forte che deve essere imposta per tutti i punti del dominio poi abbiamo le condizioni al contorno che sono un po' più semplici dell'elasticità perché è la temperatura muscolare quindi non ha una sola componente allora c'è una porzione del contorno che noi chiamiamo SQ che noi chiamiamo SQ intesa sulla configurazione iniziale in cui è il posto di flusso quindi abbiamo una condizione di Neumann quindi QD sta per dato non esattamente come nella siccità quindi su quella porzione abbiamo il flusso imposto in un'altra parte sulla superficie di rischia invece abbiamo la temperatura imposta e poi valgono le stesse identiche considerazioni che in elasticità solo che avendo una sola componente è più facile da esprimere l'intersezione di queste due superfici deve essere nulla l'unione di queste due superfici deve dare tutto il contorno del corpo cioè in ogni fisicamente cosa vuol dire che in ogni punto del nostro contorno dobbiamo imporre o la temperatura o il flusso in alternativa quindi queste sono le condizioni al contorno e poi abbiamo le condizioni iniziali cioè dobbiamo segnare il valore iniziale della temperatura per ogni x al tempo 0 deve essere uguale a una funzione assegnata t0 in tutto quello che faremo oggi t0 è sul posto uniforme quindi non dipende da x ma è una temperatura costante per tutto il corpo Va bene? Ecco, questa è la forma forte, il passaggio che abbiamo fatto nel primo giorno, abbiamo preso questa qui, l'abbiamo moltiplicata con una funzione test che adesso chiamo t tilde, abbiamo effettuato un'integrazione per parti, tra virgolette, un po' più complicata, perché c'è un'integrazione per parti con una teorema della divergenza, che trasforma questo termine, trasferisce sostanzialmente una di queste derivate la derivata della disorgenza la trasferisce al campo test e quindi abbiamo la prima parte che rimane tal quale perché semplicemente la capacità termica è moltiplicata per la temperatura test integrata sul volume quindi esattamente il problema qui invece c'è stata l'integrazione del bel parte in cui rimane K per i gradienti di T che però è moltiplicata scheramente per i gradienti di T tilde la divergenza è stata ribaltata sul campo test. Questo ha generato, poi c'è il termine di sorgente, che è semplicemente C per la temperatura test integrata, e poi c'è il termine che deriva dall'integrazione per parti, che è QD per T tilde, cioè il flusso imposto per la temperatura test. questo è come in elasticità è imposto per ogni scelta possibile della funzione test nello spazio che chiamo tau0 cioè lo spazio di quelle funzioni che sono sufficientemente continue e che si annullano dove? si annullano sulla superficie sth laddove vengono imposte le condizioni di rispetto e poi c'è questa condizione iniziale c'è una perfetta similarità somiglianza con l'elasticità infatti molto spesso anche noi abbiamo fatto prima questo problema della diffusione perché essendo un problema scalare tutto è più facile tutto è più semplice a maggior ragione adesso dopo aver visto la complessità dell'elasticità ritornare al problema scalare sembra effettivamente una semplificazione. Quindi tutto ciò che noi abbiamo fatto per la elasticità dai concetti di elementi isoparametrici, integrali elementali, assemblaggi, tal quali identici, solo che è un po' più facile. Solo che è un po' più facile. e ecco però qui ho fatto un po' il riassunto della procedura allora qui nel riguardo rosso c'è quindi la forma debole la forma debole in cui però abbiamo già anche messo l'equivalente di quello che abbiamo, che abbiamo chiamato legge costitutiva perché la legge costitutiva nel caso della diffusione dice semplicemente che il flusso Q è uguale a parte segno a K per il gradiente questa è la legge costitutiva ed è già inserita lì dentro per cui effettivamente tutto il problema è formulato in quella maniera cerco una temperatura nello spazio tau td cioè nello spazio delle funzioni che sono ammissibili con i dati tali che per ogni valore del tempo e per ogni scelta possibile di t tilde nello spazio tau 0 sia verificata questa questa identità date le condizioni iniziali poi questo è più che altro quasi un ripasso di quello che abbiamo visto in elasticità, si decide di discretizzare di scegliere T all'interno di uno spazio di elementi finiti per cui appare il Pdch e la si esprime come la somma di una funzione che rispetta le condizioni al contorno ed è una funzione che forse dovrei ricordarvi di mettere le parentesi perché forse in elasticità mettevamo le parentesi di sopra nella nella d verificare insieme a proprio sì per essere più omogenee quindi cercherò di metterle una parte una funzione che prevede che rispetta le condizioni a contorno e le stende a tutto il dominio più una parte che invece dipende dalle dalle incognite e questa è l'interpolazione globale per Per effettuare l'integrazione globale abbiamo bisogno di modificare, di arricchire leggermente la struttura dati. Perché abbiamo bisogno di dire in ciascun nodo, perché questa è una somma su tutti i nodi, una somma su tutti i nodi che moltiplica le funzioni di forma globali per le temperature nodali, abbiamo bisogno di sapere se in un dato nodo la temperatura è incognita oppure assegnata dalle condizioni al contorno. e in questo caso decido di chiamare questo campo, che è un campo scalare, Tdof della temperatura, il grado di libertà della temperatura. Se il Dof è maggiore di 0, allora è incognita e se invece il Dof è minore di 0, allora la temperatura è assegnata con lo stesso criterio dell'elasticità. Va bene? Quindi la parte assegnata sarà semplicemente la sommatoria di tutte le funzioni globali per il valore della temperatura laddove Tdof è minore di 0 e invece la parte incognita sarà la sommatoria su tutte le funzioni globali laddove Tdof è maggiore di 0. Poi, come anche in elasticità, ricordiamo che noi mettiamo, in questo caso qui, in generale metto un meno 1, non è un negativo generico, ma metto un meno 1 per convenzione. Qui invece metto in tdof, quando è un'incognita, un numero j, un intero j, che rappresenta l'ordinamento globale delle incognite, cioè io creo sostanzialmente un vettore che indico con la doppia barra T che colleziona tutti i valori nodali incogniti, tutti quanti i valori nodali incogniti. Per cui è possibile mostrare che questa rappresentazione con le funzioni globali è un approccio alla Galerkin, che può essere scritto in questa maniera, ecco qui che ho usato giustamente non so perché la sopra non l'ho usata quindi una funzione assegnata più una sommatoria di funzioni di forma globali che moltiplicano le incognite queste sono funzioni scalari a differenza dell'elasticità e dal paragone con queste due espressioni troviamo che la funzione di forma globale j non è altro che n tilde n cioè una specifica funzione di forma globale sul nodo ennesimo se j è uguale a nodi n tdof maggiore di 0. Cioè se nel nodo n la temperatura è incognita e l'ordinamento globale assegnato a questa incognita vale proprio j. Questa è la corrispondenza fra quello che abbiamo chiamato interpolazione globale agli elementi finiti e approccio alla galerchia perché insisto su questo? perché la prima lezione mi ha fatto vedere che una volta che conosciamo questa interpolazione una volta che conosciamo questa prima parte e le funzioni di forma automaticamente si butta dentro questa interpolazione nella forma debole e si ottiene la forma discreta e si ottiene la forma discreta che in questo caso ha questa espressione. Allora, il t tilde è il vettore che colleziona tutti i valori arbitrari del campo test. M, io ho dovuto mettere un pedice per non confonderlo con la massa della meccanica, M termico è la matrice di capacità. Beh, d'altronde avrei potuto chiamarla C, però C è un simbolo che viene usato più frequentemente per la matrice di conducibilità e quindi ho semplicemente detto che questa è la massa termica, mth. Da dove deriva questa mth? Deriva dalla discretizzazione del termine che contiene la derivata della temperatura. Le derivate sono imposte sulla vista dei valori nodali, t punto. C per t invece deriva dalla discretizzazione di questo operatore qua, gradiente scalare gradiente. E poi F invece deriva da questa parte. In realtà, se ci fossero, in realtà probabilmente vi ricordate che se ci fossero dei valori nodali di temperatura assegnati diverse da 0 sul contorno, questi contributi finirebbero qua dentro, all'interno del vettore F, perché sono noti e quindi finiscono lì. Quindi la discretizzazione di questo riquadro genera questa forma discreta che deve però essere imposta per tutti i possibili vettori titilde, bene, questo vuol dire che si riduce a risolvere questo sistema che è un sistema di equazioni differenziali ordinarie del primo ordine molto simile a quello della dinamica nella dinamica qui avevamo la matrice di massa qui avevamo la derivata seconda dello spostamento qui avevamo la matrice di rigidezza, qui lo spostamento e il termine di moto però è la stessa cosa quindi la procedura è questa penso che sia del tutto inutile e sia una perdita di tempo nel senso che quando voi studierete bene o ripasserete e capirete come si calcola la matrice di rigidezza e elasticità per riuscire a calcolare questa matrice di conducibilità nel problema scalare è una ipersemplificazione di prezzo quindi è inutile stare a rivedere chiaramente tutti quanti i passaggi abbiamo per scontato che siamo capaci e effettivamente è facile arrivare a questa forma discreta e allora cosa ce ne facciamo? qui a partire da questo momento è tutto identico a quello che abbiamo fatto nella prima lezione cioè prendiamo la 1.13 che è la relazione con cui approssimiamo le derivate e quel parametro θ è un parametro arbitrario variabile fra 0 ed 1 quando ad esempio ricordiamo, vediamo quando θ è uguale a 1 allora questo termine si annulla e rimane solo questo cioè vuol dire che noi approssimiamo la derivata della temperatura alla fine del passo come il valore della temperatura al quindi vabbè perché l'ho lasciato no, questo è sfuggito eh sì, sono un'astro questi vettori qua questi qui dovrebbero essere dei vettori invece sono questi scheletri strani va bene, vuol dire che usiamo questa differenza finita cioè la temperatura alla fine del passo meno la temperatura all'inizio diviso delta T per approssimare in questo caso la derivata alla fine del passo. Altrimenti, se θ è uguale a 0, allora questo sparisce e rimane solamente questo. Se θ è uguale a 0, utilizziamo questa differenza finita per rappresentare la derivata al tempo tm. Un approccio che utilizza θ uguale a 1 si chiama implicito, un approccio che usa θ uguale a 0 si chiama esplicito. Va bene, usando le equazioni di equilibrino della trama di prima si arriva a formulare questo che è l'espressione del passo temporale, cioè noto tutto ciò che succede al tempo tn, quindi nota la temperatura e nota la forzante anche alla fine del passo, si ottiene una stima della temperatura Tn più 1 alla fine del passo. L'abbiamo implementato e abbiamo visto che il caso θ uguale a 0 era condizionatamente stabile. Vi ricordate che a un certo punto se usavamo un ΔT troppo grande tutto esplodeva. invece il caso θ uguale a 1 è incondizionatamente stabile cioè qualunque ΔT noi utilizziamo la storia integrata rimane stabile magari non accurata ma stabile e ormai questo non ci sorprende più il ΔT nel caso di condizionatamente stabile ha un limite superiore che dipende dal massimo autovalore dal massimo valore del problema associato al problema termico tal quale all'elasticità poi ho incluso nelle dispense, anche se adesso non ne parleremo come si arriva a stabilire queste cose perché nel problema di diffusione è abbastanza facile è molto più facile che nell'elasticità e chiarisce un po' il senso di stabilità, consistenza e vale ancora il teorema di Lax che abbiamo citato la settimana scorsa che afferma che per questi problemi se un algoritmo è stabile e consistente allora la soluzione numerica converge alla soluzione vera tutti questi schemi sono consistenti se siete curiosi di vedere cos'è leggete la dispensa la stabilità l'abbiamo appena commentata quindi se sono sia consistenti sia stabili allora convergono alla soluzione esatta che è quello che noi effettivamente speriamo facciamo un piccolo esempio giusto per non lasciare tutti i concetti astratti adesso prendiamo una trave che è costituita da due materiali diversi e adesso andremo a vedere quali sono questi materiali A un certo punto, in questa trave si trova tutta la stessa temperatura P0, la temperatura iniziale, che possiamo immaginare 0, e a tempo 0, qui a destra, imponiamo una condizione di Dirichlet sulla temperatura, segniamo una temperatura qua. su tutto il resto del su tutto il resto del del corpo assegniamo flusso nullo qt uguale a 0 cioè in realtà non diciamo nulla ma siccome non dire nulla equivale a imporre che qt sia uguale a 0 stiamo imponendo flusso è uguale a zero. Ecco, secondo voi, fate immaginare, ripeto, siamo alla stessa temperatura iniziale, fate il nuovo materiale di Mersi, impongo la stessa temperatura a destra, anche se, questo poi lo voglio verificare perché non mi convince vista così questa figura cosa succede alla trave a regime verso che condizione evolve la nostra tra della storia della temperatura mentre provate a pensare a questo andiamo a vedere il file input che è questo il codice che ha diffuso adesso devo vedere tutto e verifichiamo effettivamente anche le condizioni al contorno parliamo cioè a vedere come sono strutturate anche questo caso il fai l'incontro Grazie. utilizzeremo per fare questa analisi prima andiamo a vedere il geo così controlliamo la definizione di tutte le superfici fisiche allora innanzitutto diffusion beam è questo qua che utilizza come file mesh lo stesso nome quindi diffusion beam geo, adesso vado ad aprire vediamo un po' come è strutturato allora la lunghezza è 10 l'altezza è 0,5 beh, deve definire due due domini che sono i due materiali per cui ha bisogno di 6 punti giustamente perché Allora, cerchiamo di capire come è definito qui, perché anche io questo giro non l'avevo guardato da un po' di tempo. Allora, definisce il primo punto che è in 0,0, poi il secondo punto è ancora in 0, meno h, quindi è partito da qui e sta costruendo questo è il punto 1, questo è il punto 2, probabilmente adesso farà così e definirà il punto 3 qua, vediamo se è giusto, il punto 3 è in l, meno h, il punto 4 è in L0, quindi è questo, poi il punto 5 è in LH, quindi questo qua, e il punto 6 è questo. Ok, quindi ha definito i sei punti che costituiscono il nostro oggetto. Poi fa le linee, che sono quindi la linea 1, la linea... no, scusate, la linea 1 è questa qua, la linea 1 è questa. poi la linea 2 poi la linea 3 la linea 3 da 3 a 4 la linea 4 da 4 a 1 con questo verso qui e poi avremo la linea 5 da 4 a 5 la linea 6 e la linea 6 poi definisce le due curve chiuse che sono costituiti da 1, 2, 3, 4 quindi questo è il loop e l'8 quindi questa è la curva 8 e poi dopo il meno 4 perché devi invertire il verso del 4 per fare il ciclo chiuso in senso antiorario quindi è un 4, 5, 6 e questa è il 9 e poi crea la prima superficie che è quella che ha come contorno chiuso l'8 quindi questa è la superficie poi crea le entità fisiche ovviamente non cambia niente che si crei una mesh per fare un'analisi di elasticità un'analisi di diffusione l'importante è definire quali sono i domini a cui associare i materiali e definire le linee l'entità su cui vogliamo imporre le condizioni del contorno quindi allora abbiamo physical line 1 e 1 più 7 che è questa qua quindi questa è la physical 1 poi abbiamo la physical 2 che è 3 più 5 quindi alla destra physical 2 e poi abbiamo la superficie 3 che è la superficie fisica scusate che è la la fisica 3 è la superficie numero 1 quindi questa è fisica 3 e fisica 4 ok questa è dunque la nostra definizione, andiamo a vedere poi all'interno dell'altro file come vengono imposte le condizioni a contorno quindi questo finisce, lo salva salva la mesh in dmar.mesh la mesh è piuttosto la mesh è piuttosto rata vediamo in questa effettivamente in questa configurazione com'è no, sbagliato open, ok si direi che è molto rada però mi fa vedere due materiali in due colori diversi beh, d'altronde questo è anche il primo esempio in cui abbiamo usato due materiali diversi ma tutti i codici che abbiamo visto finora possono essere utilizzati con il numero di materiali che si vuole basta dividere il dominio in parti, in entità fisiche diverse, entità superficie diverse, associare ogni entità di superficie a un materiale diverso, quindi è un po' più generale di quello che possiamo immaginare e oltretutto potete anche fargli dire da Gmesh quali sono le entità fisiche basta giocare un pochettino con gli strumenti che si hanno andiamo a vedere adesso il file input m, il punto m allora, innanzitutto dice di andare a cercare la mesh nel file giusto, poi definisce il materiale due righe due righe perché sono due materiali diversi il primo è K il secondo è il prodotto Rho per C è vero che sono due parametri diversi ma appaiono sempre come prodotto per cui assegniamo sono valori del tutto fittizi uno ha il materiale a K uguale 10 Rho C uguale 1 l'altro invece a K uguale 1 Rho C uguale 5 una riga per ogni elemento poi solid solid dà la corrispondenza tra gli elementi fisici e materiali quindi il primo elemento fisico è associato al materiale nella prima riga e l'entità fisica 4 è associata al materiale numero 2. Poi abbiamo delle temperature, TBC qui vuol dire boundary conditions in T grande, speriamo che, perché non l'ho mai usato effettivamente, mischiando un problema in cui ho la termica e la elasticità con forze applicate che chiamavano TBC. Devo cercare di verificare che non succedano pasticci, però non c'è pericolo, diciamo, nelle nostre piccole applicazioni questo non succede. TBC dunque sull'entità fisica numero 2 applichiamo la temperatura 10, quindi effettivamente è quello che avevo affermato prima, quindi devo capire poi per l'immagine. Qua sopra a un certo punto applicheremo una temperatura uguale a 10, tutto il resto è a 0 e vediamo come evolve la temperatura. Poi dopo dovremo dare alcuni parametri che definiscono la nostra integrazione temporale, il tempo finale è 20, cioè se noi vogliamo vedere la situazione a regime dobbiamo mettere un tempo finale lungo, altrimenti se vogliamo vedere quello che succede in un valore di tempo intermedio lo abbasteremo vogliamo procedere con delta t uguale ad 1 qui è implementato solamente il metodo implicito che è quindi sempre stabile con delta t uguale ad 1 e poi siccome creerà un piccolo filmino deve salvare l'output ogni tanto e qui gli dico semplicemente ogni quanto tempo salvare l'atto quindi genera 20 snapshot se poi sarà un filmino infilando una polata tutti vediamo cosa succede facendo girare quindi di fusion prendendo come input è proprio diffusion diffusion diffusion bim bimat. Allora vedete è rapidissimo, ecco, ah sì era solamente un artefatto. Allora vedete qui innanzitutto che è rosso uniforme cioè effettivamente sta applicando una temperatura uguale a 10 qua, la soluzione è chiaramente continua in temperatura perché non può che essere continua ma io non ho fruttato i gradienti se vedessimo i gradenti sarebbero discontinui d'altronde non c'è nessun motivo che impedisca i gradienti della temperatura di essere discontinui e poi cliccando qua sotto cliccando qua sotto su questo simbolino qua potete vedere l'animazione ecco potete vedere l'animazione e l'evoluzione poi non si ferma ma continua e va avanti. Secondo voi, vi lascio andare avanti, che cosa vi aspettate? Anche se non in maniera esatta, perché il flusso è imposto in forma debole, ricordatelo, tutte le condizioni di equilibrio sono imposte in forma debole, non puntuali, mentre le condizioni in Dirichlet sono imposte in maniera forte, cioè io impongo che su tutti questi nodi la temperatura abbia il valore che io voglio e quindi se non è 10 è perché ho sbagliato da qualche parte e le condizioni di equilibrio, cioè l'equazione e il flusso non sono imposte in forma forte ma in forma debole, ricordatevi, andate a vedere come si passa dalla forma forte alla forma debole, anche le condizioni al contorno sono imposte in forma debole, ma più o meno che cosa vi aspettate qua? queste linee che sono le linee di isotemperatura lì come devono comportarsi se tutto è fatto bene bene che acquisti la stessa temperatura è un po' come diffusione di moli, di concentrazione cioè se c'è un grado di concentrazione costante da un lato mi aspetto che a t infinito sia raggiunga lo stesso valore di concentrazione che sto imponendo da un lato Allora questa è la risposta della mia domanda, io mi aspetto che a regime tutto prenda la stessa temperatura e infatti se andremo avanti vedremo che tutto arriva a 5 o a 10, l'ho fermato troppo presto. Ma la mia domanda è un'altra adesso e si può rispondere solamente guardando il transitorio. Quando non siamo ancora arrivati a regime e quindi qui c'è la temperatura e qui c'è la temperatura manifestamente diversa, Che condizione deve essere rispettata qua? Ho detto flusso nullo, no? Questo è quello che semplicemente si impone in forma debole. Quindi, a parte il segno, il gradiente di T deve rispettare questa cosa, almeno in forma debole, almeno abbastanza bene. Quindi questo cosa vuol dire? Che la temperatura non può avere gradiente in direzione normale. Ma normale cosa? Normale al costo. Quindi la normale alcol, in questo caso, è meno Ex, quella superficie. Quindi vuol dire che la temperatura non può avere gradiente in direzione x, può averla solamente in direzione y. Questo è quello che si dice la condizione in forma forte. Non vuol dire che la temperatura deve essere uniforme qua, ma vuol dire che queste linee devono arrivare tangenti a questa superficie. ora nella maggior parte dei casi sembra sembra nero alcuni in alcune situazioni sembra che non sia così però effettivamente procediamo passo passo poi ci rendiamo conto che effettivamente è così allora vediamo un passantino vedete come noi così vado avanti così anche qua così continua ad essere anche se qui si vede si vede meno ogni volta che una linea arriva sul contorno ci arriva praticamente tangente e che è un rispetto della condizione contorno in forma debole e comunque se non sono sempre verifiche che è ben opportuno fare così come anche sulla superficie superiore inferiore non è esatta nella mesce molto molto rara per cui ad esempio questa questa curva di livello non sembra arrivare verticale però se raffinassimo la mesh ci arriverebbe con ottima approssimazione verticale comunque è abbastanza evidente che questa condizione è rispettata quindi su tutti i bordi le linee di contorno devono arrivare ortogonali. Siamo sicuri che arriva a 10 allora facciamo una prova visto che va così così veloce aumentiamo il tempo finale esageriamo, lo facciamo andare fino a 200 e salviamo e poi rilanciamo la diffusione comunque lo fa ancora super rapido e questo è l'inizio vedete che a un certo punto ormai diventa tutto rosso allora blocchiamolo e procediamo qui dovremmo voglio andare a vedere step 19 mi fermo al step 19 e poi per vedere effettivamente qual è il range ormai lo sapete bisogna andare in options e modificare il range diciamo lo mettiamo da 9 fino a 10 e vediamo se ci siamo abbastanza prima Vedete che tutte le temperature quantomeno sono comprese fra 9 e 10, allora possiamo essere un pochino più restrittivi, ancora fra 9 e 9 e 10 va tutto bene, quindi veramente ormai siamo già arrivati a regime. Solamente se gli metto una 9,999 come minimo, allora la maggior parte non è ancora arrivata veramente a 10, ma è ovvio che stiamo convergendo a 10. Dunque, quindi, a risposta a questa prima semplice domanda, tutto pian pianino si porta a... quello che mi aveva preoccupato prima è che sembra che la temperatura rossa, questo rosso intenso non sia esteso anche qui, quindi ho detto forse ho fatto un pasticcio, magari ho clipato la figura troppo breve, non lo so. In realtà, andando poi a vedere la soluzione, si possono sempre fare delle considerazioni interessanti come ad esempio la verifica delle condizioni al contorno in forza, se sono rispettate, in forma debole, in forma forte, sono tutte delle considerazioni che ogni volta che voi vedete un'analisi, fatte, realmente finite, dovete fare perché magari vi siete sbagliati a imporre qualche condizione, se non fate un'analisi critica del risultato vi sbagliate. quindi, certo se è il vostro codice è impossibile sbagliarsi ma se è il codice di un altro può succedere quindi è sempre meglio controllare ok, allora, quindi questa è la diffusione ovvero lezione 1 rivisitata in due dimensioni rimane un po' il dispiacere di non fare nessun problema in tre dimensioni ma veramente, credetemi, mentre passare da una dimensione a due dimensioni in po comporta un sostanziale cambiamento di prospettiva poi passare da due dimensioni tre dimensioni sono solo piccole modifiche e che rendono soprattutto la visualizzazione molto più complicata ma siccome non è la nostra priorità la visualizzazione poi in realtà è un po più facile è un po più facile però se un codice generale se volete se stare quando si fa una nuova formulazione si testa la formazione si fa in due di se funziona la transizione 3d è una cosa che prende mezza giornata magari invece l'implementazione di due di ha preso una settimana per verificare errori non errori insomma quindi una cosa veramente molto molto semplice allora in realtà quello che a noi interessa per i mezzi ok ma veramente importante la vista pratico che che succede quando un MEMS che è un oggetto solido viene scaldato. Quando un MEMS che è un oggetto solido viene scaldato. E allora questa è la classica formulazione della termoelasticità che non è altro che l'applicazione della legge costitutiva che noi abbiamo ripassato la lezione numero 2. cioè lo sforzo dipende anche dalla temperatura qui se andate a vedere è la stessa identica legge che trovate nel capitolo 2 ho mischiato i parametri lambda e mu con invece questa perché questa espressione è talmente standard, talmente universale che si scrive sempre con il modulo di Young e con il recente di Poisson comunque la prima osservazione vediamo di ripetere un pochettino come funziona è che innanzitutto T grande rappresenta la variazione di temperatura rispetto alla temperatura iniziale cioè l'avremmo dovuto chiamare la T-T0 ma poi scrivere sempre T-T0 è una menata quindi T rappresenta proprio questa variazione non è la temperatura assoluta non è la temperatura assoluta Poi, il secondo luogo ha influenza solamente sulle componenti di sigma che hanno gli stessi indici. Cioè, se impongo una variazione di temperatura, avrà influenza su sigma xx, sigma yy e sigma zz. Mentre tutte le altre componenti sono totalmente non influenzate. Cioè, la temperatura provoca una dilatazione omogenea, omotetica del solido. se il solito non ha vincoli si deforma in questa maniera non induce variazioni di angoli è solo una dilatazione in questa formulazione della termoelasticità classica l'equazione della diffusione è del tutto indipendente dalle sforzi e dalle deformazioni è l'equazione che abbiamo visto prima quindi prima si calcola l'evoluzione della temperatura poi si vuole vedere qual è l'effetto di questa evoluzione sul Sole in termini di sforzi e quindi si implementa questa legge costitutiva all'interno dei nostri codici è un accoppiamento non completo nel senso che la temperatura influisce sull'elasticità ma non viceversa ed è questa la grande differenza rispetto alla termoelasticità accoppiata che vedremo dopo che è quella che induce dissipazione nel solido. Quindi prima di tutto si risolve un problema di diffusione della temperatura come quello che abbiamo visto e poi si prende il campo di temperatura calcolato dal colice e lo si utilizza per calcolare gli sforzi. Ecco, io, anche se non ve l'ho detto, alla fine del nostro calcolo se andiamo a vedere il codice di diffusione faccio proprio questo guardate un po' qua qui proprio alla fine aspettate che lo faccio più grande se no non vediamo ok allora vedete qui proprio le ultime righe mi preparo a utilizzare il campo di temperatura calcolato adesso siamo arrivati al regime quindi mi aspetto che siano state salvate tutte le temperature nodali uguali a 10 perché è l'ultima analisi che abbiamo fatto e capita bene perché poi avremo un esercizio proprio su quello voglio salvarlo quindi creo un vettore che ha la lunghezza dei nodi e poi faccio un loop su tutti i nodi salvo all'interno di la temperatura e poi salvo all'interno di questo file che sarà infatti apparso nella vostra directory, save temperature punto mat questo vettore che ho appena creato e quindi lui l'ha fatto e voi avete questo file binario all'interno della vostra directory allo stesso livello del file di Fusion dovreste verificare di averlo che non l'abbia salvato ce l'avete? non potete aprire il cavinario ma verificate di averlo quindi io ho salvato questa è la fase 1 ho salvato la temperatura e adesso voglio implementarla per calcolare gli sforzi e l'ho fatto in un file che è praticamente la copia di elasticity che ho chiamato termoelasticity che voi avete tutto completo perché non è l'obiettivo di oggi però vedete che termoelasticity comincia a fare una cosa semplicissima carica dallo stesso identico file il vettore delle temperature e salva nel campo si di ogni modo che non avevamo mai trovato perché non aveva mai fatto questo problema interno però c'è anche questo campo ti non c'è solo il campo che sono gli spostamenti ma c'è anche il campo T che è uno scalare, carica il valore nodale della temperatura. Quindi l'analisi comincia, l'analisi elastica comincia conoscendo la temperatura su ogni nodo. Questo è il punto di partenza, quindi ritorniamo adesso al powerpoint. Ok, quindi punto 1 calcola l'evoluzione della temperatura risolvendo il problema di diffusione del calore se abbiamo fatto le cose giuste se la mia previsione è giusta se io faccio andare fino a questo punto il codice termoelasticity lui mi chiederà un file che sarà il volantissimo dopo lo guarderemo si chiama termoelbeamat ok lui arriva esegue come è ovvio che sia fino alla linea 43 mette la frecciata verde e si blocca lì se io vado a vedere nella linea di comando il valore delle temperature ad esempio sul primo nodo punto T lui mi dà un numero che è molto molto vicino a 10 perché avevamo visto che eravamo arrivati praticamente a convergenza non so, sul decimo nodo da un valore estremamente simile l'avevamo verificato prima quindi lui ha effettivamente preso l'ultima analisi che abbiamo fatto, quella arrivata in convergenza l'ha sbattuta in quel file l'ho riletta e ho salvato in ogni nodo il valore della temperatura a partire da questo momento non so se avete mai preso tempo per guardare gli altri codici trovate tutta la prima parte che è identica, assolutamente identica al codice elasticity La numerazione dell'incognite, l'allocazione della matrice, l'assemblaggio della matrice di rigidezza. Qui arriva l'unica differenza del codice. Il resto è esattamente elasticity. Una differenza ci deve essere per forza perché dobbiamo prendere in conto questa temperatura che prima non avevamo. quindi dov'è che entra? Entra nel fatto che, ecco qui è improprio perché A in genere tiene conto anche dell'effetto termico, ma se questa è la A che arriviamo in elasticità quindi senza l'effetto della temperatura, nella forma debole non avremo solo questo termine qua che è il genere K ma avremo anche questo termine qua che deriva, adesso cerchiamo di spiegarlo, che deriva da questo termine di temperatura. Da dove arriva? Perché? Beh, là c'è un'espressione, qua sotto c'è un'espressione che è diversa, cioè la temperatura, ok, va bene, la temperatura la ritroviamo l'Ith, il campo disperdizzato, alfa E diviso 1 meno 2 meno 1 che è il coefficiente, quello ce lo aspettiamo, però questa traccia di Y non capiamo bene un po' da dove venga, quindi lo giustifichiamo adesso ve lo spiego allora perché la la partenza il punto di partenza è che se noi esaminiamo un problema in cui com'è in questo caso qua non abbiamo forze di volume non abbiamo forze di superficie non abbiamo nulla di forze esterne allora e soprattutto l'evoluzione e cos'è che manca qua? manca l'inerzia quindi siamo ritornati al primissimo codice in cui trascuravamo l'inerzia cioè supponiamo che tutto avvenga molto lentamente allora il principio delle potenze virtuali è integrale su omega h di che cosa? di epsilon valutato sulla funzione testa che oggi chiamo U tilde per coerenza a quello che stiamo facendo che moltiplica sigma che moltiplica il contratto doppiamente con sigma valutato sulla soluzione UH che cerchiamo siccome non c'è niente altro questo deve essere uguale a 0 è l'unico termine che rimane in principio della potenza virtuale non c'è l'inerzia, non ci sono le forze quindi questo deve essere imposto per qualunque scelta di U tilde nello spazio CH0 ok? nello spazio CAC0 questo è l'unico termine che c'è è l'unico termine che c'è nel principio delle potenze virtuali allora prendiamo la nostra legge costitutiva e la sostituiamo qua per cui a sigma dovremo sostituire che cosa? l'anda della traccia di Y che dipende da U da UAC più 2 nu per epsilon che dipende da vh e questa è la prima parte quella che qui viene chiamata viene chiamata a contratto con epsilon ok, quindi c'è scritta in forma estesa e poi c'è una seconda parte che deriva dal contributo termico quindi meno e diviso 1 meno nu alfa T per il tensore metrico questo è dunque il termine che devo sostituire qua allora, la novità è dunque la presenza di questo termine devo fare il prodotto doppiamente contratto con epsilon quindi devo chiedermi avrò E per alfa diviso 1 meno 2 nu per la temperatura che moltiplica 1 contratto con y valutato sul campo test che cos'è? che cos'era la doppia contrazione? era la somma dei coefficienti poi abbiamo due tensori facciamo il prodotto ovviamente contratto questo era supponiamo che siano simmetrici perché ci sentivano la vita ijvj cioè la sommatoria di questi cioè bisogna fare la sommatoria sia su i sia su j però in questo caso la sommatoria di Sinclicata è la fatto che il primo tensore è uguale ad 1 e solamente se i uguali a j è 0 altrimenti per cui alla fine questo prodotto dovrebbe essere uguale a y1 più y22 più y33 ma questa è la traccia di y questa è la traccia di y ed è per questo che il termine legato alla temperatura viene portato a destra perché è noto perché la temperatura è nota e fa apparire la traccia di y valutata sul campo test ok? questo è il senso dunque, noi dobbiamo risolvere questa equazione non ci sono forze di volume non c'è niente al limite ci sono degli spostamenti bloccati il termine che sollecita il corpo è questo perché questi sono parametri noti questo è legato al campo test questa è la temperatura che abbiamo calcolato prima e caricato come abbiamo appena visto quindi quello genera una sollecitazione sulla struttura l'unica novità della termolatticità di un codice di termolatticità è questa passiamo un po' di tempo a guardare questo perché poi ci sarà molto utile per gli esercizi che vi proporrò di fare nella seconda parte della lezione quella sulla termolatticità completa qui andiamo quasi tutto non voglio insistere troppo cioè non voglio perdere troppo tempo non voglio passare tutta la lezione su queste cose qui però su questo facciamo un po' attenzione perché come sempre quando vogliamo integrare un termine in genere lì sappiamo che ci sono due fasi dobbiamo essere capaci di calcolare quel termine lì su di un elemento cioè fare quelli che si chiamano integrali elementari fase 2, assemblare il vettore globale l'assemblaggio, che sia questo termine che sia una forza di volume, che sia una forza di superficie è sempre lo stesso, non cambia nulla ci interessiamo solo di vedere come si fa l'integrare elementare perché ha alcune piccole novità a cui magari non siamo abituati Allora, innanzitutto vediamo che lì dentro, quindi vogliamo fare l'integrale su un elemento, no? Quindi quello che io chiamo l'integrale, che molto spesso ho chiamato l'integrale sull'elemento E, numero E grande, che sta per elemento numero E. Ok, quindi, la temperatura, la temperatura TH, la temperatura TH ricordiamoci è nota, però non è nota in ogni punto, è nota nei nodi, è nota nei nodi. e io devo fare l'integrale, quindi devo conoscere la temperatura in ogni punto. No? Allora, quindi, prima di tutto, beh, cosa potrò dire di TH? TH è localmente, su ogni elemento, TH è interpolato usando le stesse funzioni di forma. Ad esempio, se abbiamo un elemento T6, avremo le stesse funzioni di forma. quindi n e i allora sommatoria per i che va da 1 a 6 di n e i per le facciamo così uso la k k si perché usavamo la k e questo è n k vediamo di spiegare che cosa è allora nk sono le stesse funzioni di forma dell'elemento t6 che avevamo e questi sono i valori nodali vi ricordate perché uso questo intero nk strano? perché devo andare a prenderlo dalla connettività dell'elemento quindi sarà t6 e nodi di k e questo è il numero globale del nodo questa è però l'interpolazione classica. Che differenza c'è rispetto all'elasticità? Siccome T è uno scalare abbiamo che T, H può anche essere espresso come N per Te, dove Te contiene la lista delle temperature nodali, le sei temperature ed n è solamente una riga in questo caso, è una riga che contiene n1, n2, n3 fino ad n6 che moltiplica tn1, tn2 fino ad arrivare a tn6. Molto più semplice perché in elasticità se vi ricordate avevamo n che era fatta da due righe in cui la prima riga sono riempite di coefficienti dispari la seconda è più certi pari invece avendo solo una componente la n è una riga è una riga e quindi in ogni punto dove ho bisogno di calcolare la temperatura la interpolo in funzione delle coordinate parametriche poi questo è il primo vedete qui sì la notazione è un po' impropria perché n a qui non è un vettore ma è una riga quindi preferisco metterla in quella forma lì eccolo quindi la prima cosa è questo corrisponde alla prima riga ma poi vado avanti a parte i coefficienti alfa e 1-nu vedete che saranno proprio i coefficienti che passo alla routine che dovrò utilizzare quindi passo E nu alfa e passo T E che sono le temperature nodali dell'elemento oltre alla geometria dell'elemento X allora, devo calcolare un ultimo termine che è traccia di Y di U tilde U tilde è il campo test U tilde è il campo test qui è paradossalmente più simile a quello che facevamo per allora se vi ricordate quando io devo valutare y su di un campo generico ad esempio questo campo un tilde preferisco esprimerelo con le liste no, esprimelo con le liste l'ho sempre fatto con queste liste in cui ho 3 componenti che sono y1, y22 e due volte y1, 2 vi ricordate questa lista e come lo esprimavamo? lo esprimavamo così v di a che moltiplica gli spostamenti nodali beh, siccome io lo sto valutando per il campo tilde in questo caso ci metto la tilde sopra per dire che sono gli spostamenti nodali del campo test vi ricordate questa b? la stessa è la stessa identica però qui non ho y non ho non ho y ma ho la traccia la traccia che cos'è? siamo in due dimensioni la traccia non è la traccia è la somma delle componenti con indici uguali ma y33 è nullo quindi è la somma di y11 e y22 quindi se ho le due righe che danno y1 e y2 devo sommare le due righe di B sostanzialmente quindi BV sarà nella rappresentazione della maniera BV sarà una riga che sarà data da B scriverò 1,2 punti cosa vuol dire B1,2 punti tutta la prima riga più b2,2 punti tutta la seconda riga quindi se conosco b e ormai so come calcolarlo magari non lo sono dimenticato ma basta che vada a copiare basta che apra una qualunque routine di elasticità copio come si calcolava b l'unica cosa che devo raggiungere è bv uguale a b1 più b2 no? e sono pronto e sono pronto, semplicemente poi se voglio scriverlo in quella maniera un trasposto F farò il trasposto di questa relazione, quindi allora scriveremo così che la traccia di Y valutato su T è uguale a BV per UE ma è anche uguale al trasposto quindi è uguale a 2E trasposto per BV trasposto andiamo a vedere come effettivamente eccolo qua una parte è già segnata qui di lì trasposto per NA per le temperature che sono note lo jacobiano il peso e il resto si fa come sempre ma guardiamo un secondo il MATLAB l'artista che apriamo quindi quella routine che è la termo e per te questa è la stessa intestazione di prima entriamo nella routine con tutti questi elementi le coordinate del nodo in modo di angola nulla fa e le temperature nodali usiamo una regola di Gauss a tre punti e poi facciamo una sommatoria per ogni punto di Gauss e si procede come prendendo il valore di Gauss creando la lista N che contiene le sei funzioni di forma una, due, tre, quattro, cinque, sei è una riga la temperatura per la data coordinata di Gauss è semplicemente uguale ad n per t è la riga delle funzioni di forma per la colonna delle temperature nodali e poi il resto è simile, abbiamo bisogno dello jacobiano per integrare quindi si calcola prima D, poi la matrice jacobiana e lo jacobiano è il determinante della matrice queste sono identiche a qualunque routine di elasticità poi si calcola l'inverso perché abbiamo bisogno di questa parte? per arrivare a y per arrivare al calcolo della b quindi abbiamo bisogno della matrice jacobiana inversa del gradiente e poi invece di calcolare prima la b e poi fare la somma delle due righe visto che la devo definire vado a fare la somma direttamente delle due righe potete andare a verificare guardando ad esempio all'interno della matrice elementare di rigidezza che la somma delle prime due righe di B è uguale a questo oggetto qua. È una riga sola, i tre punti vuol dire che va a capo ma non sta riuscendo la riga, quindi è una riga di 12 elementi che sono uguali alla somma delle due righe di B. E poi applica l'integrale numerico, come abbiamo visto prima, a parte il coefficiente che dipende dai parametri di V trasposto per la temperatura per il determinante, per lo jacoviano per il peso di Gauss e si fa l'intervento niente di assolutamente eccezionale quindi fatta questa modifica siamo pronti a vedere qual è l'effetto di questa cosa sulla sulla nostra trave aspettiamo ma lui sceglie automaticamente ah no, l'avevo già scelto io lo faccio andare prendo quindi termoel BIMAT ecco questo è il campo di spostamento finale componente 1 componente 2 e poi gli sforzi gli sforzi chiaramente sono discontinui perché c'è l'interfaccia se volessi un miglioramento della discontinuità dovrei raffinare la mesh ok, quindi questa è la soluzione ma associato a questo andiamo a vedere come è fatto il file di input termoelbinbmat insomma abbastanza abbastanza ah ma non me lo lascia se non lo chiuso qua ok i parametri elastici sono gli stessi per i due materiali 10 e 0 poi paragoneremo con una soluzione analitica che è nota per lungo e a 0 quindi metto lungo e a 0, 10 e 0, 10 e 0 e questi sono i due alfa, uno e il doppio dell'altro associo superfici fisiche a materiali impongo delle condizioni al contorno cerchiamo di capire cosa sto imponendo sull'elset fisico numero 1 che se vi ricordate era un bordo a sinistra in direzione 1 impongo spostamento 0 cioè impedisco che ci sia spostamento in direzione X sulla superficie di sinistra però in teoria questo potrebbe traslare verso l'alto se mi limitassi a imporre questa condizione e allora siccome devo filtrare sempre i moti rigidi se no la matrice K è singolare impongo una condizione concentrata sul nodo 1 che è il primo in basso anzi non è il primo in basso, è quella metà quella metà del nodo 1 in direzione 2 uguale a 0 cioè blocco lo sfruttamento verticale di quel nodo e così blocco il moto rigido di traslazione verticale perché ha fatto così? per riprodurre le condizioni dell'esempio, della soluzione analitica di cui adesso parliamo. Allora, c'è, quindi, ritorniamo sul... Ecco, questo è il esercizio che però non facciamo adesso, ma vi lascio da fare con calma per conto vostro. Ci sembra la stessa trama, i due materiali in condizioni di plane strain, la lower, la parte inferiore di avere un'alfa maggiore della parte superiore per ipotesi se no questa soluzione non va la soluzione analitica è nota nel caso in cui la temperatura con cui la trave viene riscaldata sia uniforme cioè questa TD sia imposta uniforme ma è un po' la situazione in cui siamo adesso perché abbiamo lasciato andare avanti la diffusione fino a che tutti i nodi erano riscaldati di 10 gradi quindi nel nostro nella nostra situazione questo td e quella 10 c'è uno riscaldamento uniforme per tutta la trave 10 gradi e qui potete calcolarvi lo spostamento lo spostamento indotto da questa distribuzione di temperatura e potete verificare con la soluzione numerica che avete ottenuto se questa concido con cd chiaramente con un'ottima accuratezza, chiaramente la temperatura è continua, non c'è discontinuità, scusate lo spostamento è continuo, non c'è discontinuità, abbiamo visto che le discontinuità sono sugli sforzi. Bene, questo è una verifica della bontà dell'approccio del codice, però la cosa più importante è tenere bene a mente che se abbiamo diversi materiali, anche se il vincolamento di questo è molto debole perché vi ricordo che abbiamo solo imposto che non ci siano spostamenti in questa direzione qua e che questo nodo non si sposti verso l'alto quindi uno potrebbe dire, va bene ma si espande che problema c'è? perché mai si dovrebbero generare degli sforzi? invece a seguito del fatto che i due materiali hanno un coefficiente di dilatazione termica diverso si generano degli sforzi importantissimi e siccome i MEMS voi sapete che hanno degli strati soprattutto i MEMS piezo hanno lo strato di silicio, lo strato di ossida, lo strato di LN, lo strato di molideno per l'elettrodo sotto, lo strato di molideno per l'elettrodo sopra o alluminio, quello che è, c'è uno stack di 6 o 7 materiali diversi, quindi questa è la situazione standard, materiali diversi, coefficienti di irratazione diverse, sforzi importanti che possono anche indurre la gelaminazione di questi strati. per cui valutare bene secondo voi come si deformerà questo oggetto qui noi siamo a temperature uniforme alfa maggiore qua e alfa minore sopra cosa farà? certamente sarà una flessione qui rimane fermo e va su a banana così verso l'alto potete vederlo vi ricordate come si fa? Dovete usare il plugin scalture e creare il campo avventorial sempre quello facciamo così giusto per ricordarci come si può vedere una deformata una deformata che può essere sempre utile faccio andare la soluzione ok quindi qui non c'è la deformata sono solo le due componenti quindi si va su plugins tools plugins e si va in fondo a cercare quello che si chiama scal to vec dagli scalari ai vettori bisogna selezionare qui a destra le componenti del vettore quindi la componente 0 che è un 1 la componente 1 che è un 2 e poi se clicca run per far andare il basso a destra per far andare il plugin ora qui al momento non si capisce niente uno perché ci sono due viste sovrapposte quindi tolgo la vista del campo 1 purtroppo per default lui usa le frecce per visualizzare, già si capisce un pochettino, le frecce vi dicono come si deforma ma se lo vogliamo vedere nella maniera classica dobbiamo andare su options e modificare l'azione di vista dell'ultima di view 5 che abbiamo appena creato bisogna fare due cose uno andare sull'aspetto l'aspect e a metà c'è questo vector display invece di utilizzare le frecce 3D selezioniamo spostamento e già migliora poi però siamo abituati di più a usare le isoline quindi andiamo ancora in tools options sempre 5 in general nel primo quindi intervals type selezioniamo invece che continuous map field, iso values e ad esempio io tipicamente uso una ventina di livelli e vediamo che questa è la visualizzazione quindi qui è abbastanza evidente che chiaramente abbiamo usato dei parametri che è tutto fittizio per cui la deformata è esagerata però vedete che bene che effettivamente esiste ed è anche molto importante. Quindi, MERS, effetti termici tremendi, brutta roba, bisogna sempre calcolare. Quindi la termoelasticità, un'analisi termoelastica in un MERS all'ordine del giorno. Quindi abbiamo visto due cose che si fanno in un MERS, l'analisi modale e analisi termoelastiche per vedere gli effetti di una variazione di temperatura ovvio le strutture sono più complesse però il concetto è questo la complessità la gestisce il codice di elementi finiti lo vedrete nel ciclo di esercitazioni che farete più avanti però quello che ci interessa oggi e soprattutto migrare alla dissipazione termoelastica. La dissipazione termoelastica è un concetto importantissimo soprattutto in tante classi di MEPS, ma in particolare nei risonatori. I risonatori, tantissimi MEMS risuonano, quindi uno potrebbe dire che non è un nome azzeccato, perché i giroscopi risuonano, i microspecchi risuonano, tutti insomma, a parte gli accelerometri che sono ad alta pressione e quindi hanno un comportamento quasi statico, quasi tutti i MEMS che noi trattiamo, i MEMS inerziali che noi trattiamo sono risonanti, però quando si dice risamatore, si parla di MEMS che sono proprio creati solo perché risuonino. E risuonino ad una frequenza precisa da utilizzare, ad esempio, come clock nei protocolli di sincronizzazione. Nel 5G ormai molti MEMS vengono utilizzati, le alte frequenze. Perché vengono utilizzate? Perché gli si fanno fuoco. E sono molto più, hanno un Q molto più alto dei circuiti elettronici, perché anche i circuiti elettronici possono essere risonatori un RLC è un risonatore ma i MEMS sono infinitamente più accurati e dissipano di più quindi però dobbiamo cominciare a capire oggi è una prima lezione sulla dissipazione però poi insisteremo anche se non abbiamo la possibilità di analizzare troppe cose però durante il ciclo di esercitazioni con console prenderete proprio un risonatore attualmente è effettivamente un progetto grosso che stiamo facendo con ST il progetto è l'analisi di una nuova famiglia di risonatori con un nuovo materiale dopo vi dico che cosa deve fare questo materiale e quindi in questo caso per i risonatori le dissipazioni sono di due tipi lo vedremo dopo la dissipazione termoelastica e dissipazione agli ancoraggi con console cercherete di stimare tutti e due la dissipazione agli ancoraggi è troppo difficile per introdurla in questo corso se volete può essere un argomento da investigare come progetto d'esame, ci sono degli articoli che bisogna leggere ne abbiamo fatti tanti però spiegare come funzionano esagerato però ormai sono stati implementati anche dei codici console ormai ce l'ha anche lui e quindi con console si riesce a stimare molto accuratamente la dissipazione di un MAPS di un risonatore che è stato appena progettato e che poi verrà prodotto quindi in questo caso sono due ma partiamo dalla base e prima diamo alcune definizioni perché oggi introduciamo tanti concetti probabilmente se avete sentito Gorgigliano avrete già una certa familiarità però adesso alcune equazioni hanno più senso prendiamo un risonatore massa molla smorzatore che è il primo là sopra ma mu2.cu.ku uguale a ft dove ft lo prendiamo con una forza una forza anche sinusoidale e voi direte ma come? siamo ritornati a trattare di un caso 1D un DOF solo ormai noi abbiamo imparato che una struttura è data dalla somma di tanti risonatori uno di basta proiettare una base modale, quello che abbiamo fatto settimana scorsa quindi noi stiamo semplicemente dicendo che abbiamo preso la struttura, abbiamo calcolato tutti i modi che ci interessavano e abbiamo proiettato le equazioni su questi modi e stiamo analizzando un modo particolare l'equazione di un modo particolare quindi ha un significato molto più generale non è un semplice esempietto di una massa mola smoltatore è l'equazione di movimento di un mesh proiettata sul modo di interesse ok allora noi sappiamo che la soluzione di quella di quell'equazione esprimibile in questa forma come si fa? lo scrivo perché è come si fa a ottenere in modo rapido lo scopo che poi dovremmo fare un ragionamento su questa cosa qui innanzitutto la riscrivo in due forme diverse quindi il primo è m2t più c1 più kv uguale a f per il seno di omega t che è quella che c'è qui in palazzo poi molto spesso si divide per la massa e allora cominciamo ad avere U2 punti più C diviso M per U punto più K diviso M per U uguale quindi F diviso M insieme a A negativo va bene questo viene chiamato omega 0 quadro omega 0 quadro uguale a k diviso m perché si dà questo nome? perché se è assente l'anticipazione questa è la sequenza propria di questo modellino di questo modellino perché se U ha una risposta sinusoidale qui avremo meno omega quadro più K per U non ci sono forzanti quindi diventa omega zero quadro val K diviso M quindi omega zero è la frequenza di risonanza in assenza di dissipazione in assenza di dissipazione quindi questo lo chiamiamo omega zero quadro lo chiamiamo, scusate, omega perché non so più come si può usare omega zero quadro per U questo lo lasciamo scritto così con due punti e poi questo invece uso scriverlo in questa maniera probabilmente l'avete visto meccanica nazionale non so se è fatto più meccanica nazionale ma in questa forma cioè invece di scrivere c diviso m invece di scrivere c diviso m si introduce questo coefficiente nu che è una misura di disfazione e quindi c diviso m diventa essere uguale a 2 nu per omega 0 questa è la definizione stessa di questo coefficiente va, lo accettiamo poi andiamo a destra e troviamo f diviso m però f diviso m può essere scritto in questa forma facciamo apparire 2 volte k quindi avremo k diviso m per f diviso k K si sente lì chiaramente quindi è uguale a quello che ho scritto sopra allora quindi K diviso L lo ritroviamo ancora una volta è omega 0 quadro non ho cambiato niente per adesso devo ancora scrivere stendio omega T scriviamo tutto così è uguale, ho solo diviso per K e rimonturato per K questo è omega 0 quadro quindi scrivo più omega 0 quadro e poi F diviso K attribuiscono un significato fisico molto semplice se siamo in statica cioè se non ci sono accelerazioni e velocità F diviso K è quella che chiamo soluzione quasi statica a parte la presenza della omega ma se ho una forzante statica F diviso K è lo scostamento stato per questo viene chiamato come c'è scritto la U statico U st è una definizione quindi U statico è stato definito U statico è stato definito come F diviso K e questa è la forma più frequente per presentare questo oggetto innanzitutto abbiamo già introdotto vari coefficienti vedremo tra poco che questo 2NU non ho ancora spiegato ma sarà uguale a 1 diviso q che è il coefficiente più importante per noi oggi il coefficiente di qualità daremo varie definizioni a questo coefficiente di qualità allora va bene per trovare la soluzione di questo oggetto si passa sempre in variabili complesse cioè si aggiunge a questo si moltiplica per l'unità immaginaria questa e si somma l'equivalente con il coseno e si ottiene 2 nu omega 0 z z punto più omega 0 quadro z uguale omega 0 quadro statico per e alla i omega t questo passaggio è classico quindi non è molto applicato, cioè una combinazione di due equazioni in cui la prima è quella che ho scritto l'altra è quella con il coseno questa la moltiplico per i la sommo e quindi z sarà la variabile complessa associata e di conseguenza se guardo la forma la soluzione di questo sarà certamente z uguale a una z grande per e alla i omega t assumo che la soluzione sia di questa forma dove z è una linea complessa quindi z grande risulterà essere uguale a 1 diviso meno omega quadro più 2 nu i omega omega 0 più omega 0 quadro per omega 0 quadro uno statico ok, questo è la soluzione sostituendo questa ipotesi all'interno della nostra equazione poi non si tiene mai omega 0 quadro al numeratore quindi si divide tutto per omega 0 quadro e quindi avremo che questo diventa 1 e qui abbiamo meno omega 0 quadro e qui avremo omega diviso omega 0 e quindi questa è la soluzione per z complesso da cui, prendendo parte reale e parte immaginaria, si arriva a queste due espressioni. Magari era del tutto inutile fare questo sviluppo, però almeno ci ritornerò a questa equazione tra poco, a questa equazione tra poco, quando cercherò di calcolare la frequenza propria di questo oggetto che non è più omega0 quadro perché c'è anche la dissipazione quindi per effetto di questa dissipazione di questo termine che è un dissipatore come cambia la frequenza propria va bene però lo vediamo dopo allora memorizziamo solo questo che la soluzione in termini di U viene spesso espressa in questa maniera come n che è un fattore di amplificazione per u statica che sarebbe la risposta in presenza di una forzante statica per il seno di omega t meno phi phi è la fase phi è la fase e omega è la frequenza cioè u è una funzione sinusoidale alla stessa identica frequenza della forzante forzo il sistema con una frequenza u risponde con la stessa frequenza e la differenza è che però è una sinusoide con una fase che presenta una fase, cioè uno shift della risposta n è il coefficiente che moltiplicato per la risposta statica dà la risposta dinamica e i due diagrammi che tipicamente importano sono questi che presentano il coefficiente di amplificazione n e la fase p questi sono i diagrammi tipici che abbiamo visto tante volte però quelli che non interessano sono i diagrammi per nu piccolo perché nei MEMS nu è sempre molto piccolo quindi tutto quello che faremo sarà sotto l'ipotesi che nu sia piccolissimo il fattore di qualità nei MEMS a parte gli accelerometri che hanno delle disfrazioni alte ma volontariamente qualità sono almeno 1000 nei MEMS quindi il nu è piccolissimo visto che 2 nu uguale 1 su q siamo all'ordine di 10 a almeno 10 a meno 2, 10 a meno 3 quindi il nu è sempre molto piccolo allora questo dunque plotta n in funzione del coefficiente del coefficiente nu e quindi a noi interessano le curve molto alte e questa invece è la fase la fase parte sempre da 0 cioè quando la frequenza omega è molto piccola rispetto a omega 0 la risposta è quasi statica e quindi non c'è differenza di frase fra la forzante e la risposta io applico una forzante con una frequenza bassissima e il sistema risponde immediatamente senza avere uno shift poi invece quando siamo in questa condizione omega diviso omega 0 uguale a 1 cioè la omega è uguale alla frequenza di risonanza del sistema non smorzato, allora arriviamo a picco e mezzi e poi dopo se continuiamo arriviamo a 180, la fase. E per i nostri problemi con un U molto piccolo, questa spezzata è molto simile alla realtà, cioè si mantiene basta basta basta basta basta passa a 180 rapidissimamente e resta 180 allora alcune considerazioni cosa sono omega max e del max? sono la omega in corrispondenza della quale la risposta risulta essere massima cioè omega max da queste espressioni si ottiene in funzione di omega 0 non è esattamente omega 0 ma viene variata, viene modulata dalla presenza della dissipazione solo che se la dissipazione è molto piccola omega max è molto simile a omega 0 cioè, cosa vuol dire? il picco, vedete che qui questo picco non è esattamente omega 0 è un po' citato a sinistra ma più queste diventano alte più lo diventa piccolo più il picco capita in corrispondenza di omega uguale a omega 0 questo è il senso possiamo fare anche un'espansione in serie ma vedremo che omega 0 è uguale a omega più un termine quadratico quindi se nu è dell'ordine 10 a meno 2 abbiamo detto quadrato 10 a meno 4 quindi veramente correzioni molto piccole per noi e poi in quanto vale n max cioè n in corrispondenza del picco il fattore di amplificazione è dato in maniera esatta da questa espressione ma se nu è piccolo lo possiamo approssimare come 1 diviso 2 nu andiamo a vedere se è vero allora, supponiamo io ho detto che questo n è il rapporto fra il valore massimo e il valore statico il valore statico è questo a frequenze bassissime n è sempre uguale a 1 a frequenze bassissime andiamo a vedere qua, nu uguale a 0,2 fate il conto voi con la calcolatrice quanto viene questo n se nu uguale a 0,2 1 su 0,4 quindi verrà 2 qualcosa 10 diviso 4 uguale 2,25 e questo è un pochino più alto perché nu è ancora grande quindi non vale questa approssimazione in maniera esatta ma se andiamo a prendere ad esempio questo che è nu uguale a 0,1 quindi questo deve essere esattamente 5 effettivamente il picco accade per 5 quindi più nu è piccolo più effettivamente questo n ha il significato di coefficiente che prendiamo per moltiplicare la risposta statica e ottenere la risposta dinamica e poi notate che questo n max è ovvio anticipato che verrà già chiamato Q quindi Q è il fattore di qualità che avrà un significato legato alla dissipazione ma è anche quel coefficiente che moltiplicato per la risposta statica ti dà il massimo della risposta dinamica quindi ha tanti significati Q ne abbiamo tutti due va bene, questo è un po' un preliminare, forse un po' troppo lungo ma la dissipazione qua dove entra questo sistema ha un dissipatore quindi è la massa con la rigidezza e col dissipatore e questa è la forzante che sto applicando per tenerlo in movimento oscillante devo mettere un po' di energia se non ci fosse il dissipatore abbiamo imparato l'altra volta che questo movimento alla sua frequenza propria potrebbe continuare all'infinito senza mettere energia. Ma siccome c'è dissipazione, devo metterci energia. Quindi come faccio? Io ho messo la soluzione, quindi cosa so? So che la soluzione è una sinusoide. Quindi faccio così. Calcolo la dissipazione che devo mettere nel sistema come l'integrale da 0 a t grande, che è il periodo quindi l'integrale della forza che devo applicare questa qui questa quantità qui la chiamo globalmente f piccolo è proprio f piccolo uguale a f sessantino di t per la velocità io sto applicando una forza quindi la potenza è forza per velocità l'integrale di questa potenza nel ciclo è la dissipazione perché se non dissipasse di questa integrale sarebbe 0. Va bene, allora cosa si fa per fare questa integrale? Si sostituisce ad f il termine di sinistra, quindi l'equazione mu2.c1.qk1 per un punto e si sviluppa questa integrale, che per fortuna non è troppo difficile da fare, visto che so che la risposta è sinusoidale. guardiamo ad esempio il primo termine ho m2 punti per un punto ma questa è la derivata del termine m1.4 metteteci davanti l'un mezzo e si fa la derivata rispetto al tempo che cosa ottiene? il 2 va a eliminare il denominatore quindi avrò m1. per 2 punti quindi esattamente questo termine qua quindi l'integrale di questo termine è questo che dovrò valutare all'inizio e alla fine del periodo. Poi la stessa cosa vale qua, KU per un punto è la derivata di questo termine, con un mezzo davanti. Quindi questo integrale basta valutarlo all'inizio e alla fine del periodo. Ma siccome la risposta è periodica, un punto e U assumono sempre lo stesso valore all'inizio e alla fine di un periodo, per cui questo termine sparisce. cioè infatti non dipende da c se non ci fosse il dissipatore la potenza sarebbe la dissipazione sarebbe zero come ci aspettiamo se non ci fosse dissipatore se ci fosse solo la massa e la cappa la potenza sarebbe zero l'unica cosa che mi si fa è effettivamente il dissipatore c un punto al quadrato ma ci un quadrato utilizzando l'espressione di un lo si può valutare e viene fuori questa quantità credetemi ma è facilissimo da fare d'altronde se prendiamo un seno abbiamo anche la luce non devo farne la derivata quindi apparirà la omega ma è al quadrato quindi appare la omega al quadrato l'integrale di seno al quadrato su di un periodo siccome se non ho dato un cos quadrato da 1 e quindi i due sono simmetrici su un periodo l'integrale è un mezzo per il periodo, quindi è t per un mezzo allora un mezzo viene da qui e questo è il periodo t moltiplicato per l'ampiezza del movimento uguale quindi, vabbè, insomma non è difficile da calcolare questo è la dissipazione che dipende effettivamente solo e esclusivamente dalla dalla presenza del dissipatore. Quindi per questo è fondamentale che quel termine C grande sia piccolo, perché più piccolo è, meno energia devo mettere nel sistema per tenerlo in vibrazione. Poi, questo è dunque il ΔW, la quantità di energia che devo mettere nel sistema. poi l'altra quantità che mi interessa è W, la massima energia elastica immagazzinata nel sistema essendo un oscillatore sarà dato da un mezzo k per quadro dove k è la molla quindi l'energia elastica all'interno della molla questo l'avete visto già alle superiori era il classico del cito delle molle un mezzo k per quadro è l'energia immagazzinata nella molla ok cos'è il fattore di qualità? per definizione è questo cioè la definizione vera del fattore di qualità è questo oggetto, è 2π per W diviso ΔW cioè è 2π per la massima energia elastica immagazzinata nel sistema diviso la dissipazione per ciclo questa è la definizione quindi quando si chiede cos'è il fattore di qualità, è questo poi troveremo altre 10 definizioni di Q ma che sono conseguenza di questa cosa questa è la vera definizione per cui la prima conseguenza è che se prendiamo W e ΔW che abbiamo ottenuto qui la prima cosa è questa Q è K diviso Cω poi però siccome noi sappiamo che nei nostri problemi torniamo qua, nei nostri problemi U è molto molto piccolo quindi quello che dico adesso mentre la definizione del dato è generica quello che dirò adesso è specifico per i nostri sistemi dove il nu è molto piccolo quindi omega è uguale a omega 0 quindi mettiamoci in condizione di omega uguale a omega 0 e allora utilizzando il fatto che omega 0 è questo omega 0 è questa cosa qua otteniamo anche questa seconda definizione omega qui dovreste vedere nelle dispense è corretto ma nella slide è rimasto dove scrivere omega 0 per m diviso c e beh, allora se andiamo a vedere qua se andiamo a vedere queste definizioni troviamo che effettivamente q è quello che ho detto prima q è uguale a 1 diviso 2 nu ora q è uguale a 1 diviso 2 nu se nu è piccolo, cioè vi ho detto attenzione, quello che dico adesso è valido solamente se la dissipazione è piccola quindi in generale questa equazione nel mondo membro non la troverete mai scritta così ma la troverete sempre scritta così omega 0 diviso q per un punto si spone già che la dissipazione sia piccola e se qualcuno riscrive quella relazione con dissipazione grande è perché non sa da dove arriva il q non sa che cos'è q perché questo è vero sempre ma questo passaggio è vero solamente se mu è piccolino ok? allora quindi questa è una cosa importante abbiamo quasi finito la parte la parte noiosa poi quindi q è circa uguale a 1 quindi q capite adesso perché scrivo circa uguale a 1 diviso 2 nu perché questo è valido solamente sotto l'ipotesi di nu piccolo e poi anche questo capite perché n max è uguale a q circa sempre perché è una successione di identità quindi l'amplificazione il rapporto tra risposta dinamica e risposta statica è q e questo è rigoroso nei sistemi lineari ma vale abbastanza bene anche se avessimo delle non linearità di cui non abbiamo mai parlato ma magari verso la fine del corso ne parliamo quindi abbiamo visto queste due o tre definizioni di Q e quella che interessa più di tutte queste sono importantissime nel senso che ogni giorno parlando di MEMS parlate di queste cose Q come massima amplificazione Q come dissipazione tutti questi concetti sono importanti ma ai nostri fini per calcolare la dissipazione termoelastica ce ne serve un'altra ce ne serve un'altra e ve la spiego in maniera un po' qualitativa però va molto bene, credete voglio calcolare a partire da questa equazione le nuove frequenze proprie del sistema omega 0 è la frequenza senza dissipazione ma io voglio calcolare la frequenza propria con la dissipazione quindi voi sapete che per calcolare le frequenze proprie devo assumere che non ci siano forzanti quindi tolgo la forzante e trovo la soluzione una soluzione sinusoidale di questa equazione in assenza di forzante allora assumo dunque che la risposta sia di questo tipo esattamente come prima assumo di avere quindi una risposta di questo tipo in assenza di forzante sostituisco lì dentro e che cosa ottengo? ottengo meno omega quadro lo z apparirà lo z grande apparirà ovviamente dappertutto quindi non lo scrivo nemmeno ottengo meno omega quadro più 2i nu omega 0 omega più omega 0 quadro uguale a 0 quindi questa è l'equazione che risolta mi dà la frequenza la nuova frequenza propria allora lasciatemi fare un po' una porcheria io so che omega è circa omega 0 questo lo so perché so che il mio sistema è altissima poco allora lo riscrivo così omega quadro è uguale a 2 I nu omega 0 omega 0 omega meno no scusate più omega 0 quadro questo ho solo spostato a destra alcuni termini adesso faccio la porcheria siccome so che omega è circa uguale a 1 a 0 sostituisco qui sulla destra a omega sostituisco omega 0 e quindi ottengo omega 0 quadro qua e anzi raccolgo omega 0 ma noto che c'è una parte reale e una parte immaginaria quindi la nuova omega quadro è diventata complessa e il fatto che ci sia un termine complesso riflette la presenza di dissipazione la parte immaginaria non ci sarebbe se non avesse dissipazione e vediamo un po' come faccio quindi a stimare il Q dove lo faccio apparire il Q? lo vedo nascosto qua 2NU è legato al Q l'abbiamo imparato quindi 2NU è 1 diviso Q quindi lo riscrivo in questa maniera avrò I per 1 diviso Q omega 0 quadro più omega 0 quadro allora adesso facciamo il conto di questa quantità prendo la parte reale di ω² e la divido per la parte immaginaria di ω² la parte reale è questa la parte immaginaria è questa questa è una identità molto molto utile nei metodi numerici perché i metodi numerici cosa fanno? calcolano la nuova frequenza propria del sistema e prendono la parte reale dividono per la parte immaginaria e dai già la stima di Q automaticamente è quello che faremo noi è quello che faremo noi tra poco ecco vediamo questo dopo che ore sono? ho perso il controllo 4 e mezza facciamo una piccola pausa e facciamo la parte finale però vi illustro abbiamo capito che la dissipazione c'è che il fattore di qualità Q è un concetto importantissimo nel MAPS e questa è un'immagine che usiamo spessissimo per spiegare come mai non parte questo ci dovrebbe essere una piccola animazione ma non parte, pazienza questo vediamo se parte da qua come la sono persa copiando mi sa che mi sono perso l'animazione vabbè, faceva vedere che la massa bianca e come funziona questo oggetto oscilla così questi sono ancoraggi queste sono delle molle che vanno qua sopra e poi tornano giù da una parte e dall'altra quindi questo oggetto ha la possibilità di oscillare solamente in questa direzione qui ci sono dei capacitori a comb finger si chiamano che generano una forza in direzione orizzontale pressoché costante almeno diciamo in prima approssimazione costante, quindi questo oscilla questo oscilla e questo oggetto viene preso viene messo all'interno di una camera vuota e viene testato con pressione a pressione variabile viene messo in movimento e si misura sperimentalmente il fattore di qualità a pressione variabile e si traccia un diagramma sperimentali quindi del fattore di qualità che ha sempre questo aspetto qui siamo siamo in millibar quindi 10 alla terza siamo a un'atmosfera qui la dissipazione tipicamente dovuta all'interazione col gas ad alta pressione si usano le equazioni di Nadi Stokes in questo caso con delle particolari condizioni al contorno che si chiamano sleep ma ne torneremo su un'altra volta sleep boundary conditions, poi la pressione diminuisce, un accelerometro non lavora mai a un bar, lavora a qualche frazione di bar, quindi qua in mezzo, quindi se volete stimare la dissipazione fluida in un accelerometro prendete un codice negli stokes, ci mettete le sleep boundary conditions e li calcolate la dissipazione. Poi vedete che questo aumento aumenta, aumenta e poi c'è una lunga regione in cui questo diventa praticamente lineare. e che cosa succede all'interno di questo? è il regime dove lavorano i giroscopi i magnetometri, i risonatori c'è un alto livello di vuoto perché qui siamo arrivati intorno al millibar il millibar è una pressione piuttosto bassa e difficile da tenere se non in una camera a vuoto il pacchettino in cui è messo il risonatore deve avere i famosi getter dentro per riuscire a mantenere un livello di pressione i getter fatti dalla size getter cioè quelli che mangiano le molecole non riescono solo a mangiare gli inerti ma tutto il resto lo mangia e quindi tengono una pressione abbastanza bassa intorno al millibar o frazioni di millibar e qui ahimè 10 anni fa no anzi purtroppo più di dieci anni fa era quasi vent'anni fa e sentivo e dicevo ma come facciamo a stimare la dissipazione in questo regime e beh, vedi che lì c'è ci sono effetti di rarefazione vi racconterò questa cosa qui, lì si usa l'equazione di Boltzmann e a Politecnico c'era Ciorcignani che era il guru mondiale dell'equazione di Boltzmann e studiava il rientro delle navicelle spaziali quando impattano l'atmosfera e usava le stesse identiche tecniche che poi abbiamo usato insieme a lui per stimare la dissipazione dei meccanismi. Gas dinamica era rifatta si chiama free molecule flow poi vedete che però la cosa interessante è che i fattori di qualità sono perfettamente lineari in quest'ampio sono quasi due ordini di grandezza di pressione. Tantissimi MEMS stanno qua dentro. Poi se continuate a diminuire la pressione, diminuire la pressione, diminuire la pressione, ma lì ci vuole veramente una pompa vuota, perché nessun MEMS riuscirà a garantire 10 alla meno 2 mili barogenti. Però se uno è proprio interessato a vedere come va, a un certo punto arriva un plateau. Cioè vuol dire che il fattore di qualità, la dissipazione del MEMS, non dipende dalla pressione. Questo vuol dire che non dipende dal filo. non dipende più dal gas, ma che ci sono altri fenomeni che diventano importanti, che diventano esclusivamente sorgenti di dissipazione, che sono ad esempio la dissipazione termoelastica, le anchor losses, cioè le anchor losses, giusto per non lasciare questo nell'alone di mistero, sono energia elastica che migra dagli ancoraggi nel sultrato e poi non ritorna fin dietro. Quindi è una perdita di energia elastica, di onde elastiche che vanno nel sultrato e siccome il soltanto è grande le anchor sono piccole non torna più indietro si dissipa prima di tornare indietro queste sono le anchor losses e tipicamente nei risonatori la somma di termolastico e di anchor losses spiega cosa succede poi vi dirò qualcosa anche su quelle altre che sono un po' più esotiche però questo diagramma viene fatto questo plot è importante conoscerlo perché se dovete sapere qual è la sorgente dissipazione è fondamentale sapere a che pressione lavora il vostro oggetto allora oggi termoelastico termoelastico è copiato termoelastico è copiato vuol dire una cosa molto semplice l'equazione che abbiamo visto prima l'equazione termoelastica in cui sigma dipende non solo da y ma anche dalla temperatura è quello che abbiamo appena visto no? allora è accoppiato con l'equazione della diffusione questo è il capitolo 6 delle dispense come sempre non c'è niente che non ci sia a parte quelle immagini di prima accoppiato con questa è l'equazione della diffusione che abbiamo analizzato prima non c'è la sorgente di volume ma è stata sostituita da questo oggetto qua allora guardiamo un po' che cos'è T0 è la temperatura di riposo dell'oggetto cioè stiamo lavorando a 25° beh, allora sarà 273 più 25 t0 α è lo stesso coefficiente di prima è il modulo di Young nu è il coefficiente di Poisson e questa è la traccia di y punto cioè la traccia della derivata rispetto al tempo di y fisicamente che cos'è questo oggetto? la traccia di y è legata alla variazione di volume l'abbiamo visto velocissimamente nel ripasso della seconda lezione non ci torniamo più, andate a vedere quel ripasso la traccia di y vi dice quanto sta variando il volume elementare è legato alla j punto allo j al derivato di j è solo che non è la variazione di volume ma è la velocità di variazione di volume questo vi dice che se la traccia di y è positiva cioè se il volume sta aumentando se il reticolo cristallino si sta espandendo questo agisce come una sorgente negativa di calore cioè un sink lo si chiama un heat sink si espande e raffredda se invece si comprime riscalda e questo è valido in tantissime situazioni per cui non è strano non è strano solo che quel termine lì è in genere così piccolo perché è così piccolo? perché alfa è nell'ordine 10 a meno 6 e quindi rispetto agli altri coefficienti tipicamente questa cosa qui gioca un ruolo completamente trascurabile, diventa importante nei nostri casi, perché? Perché le frequenze sono alte, quindi quando io derivo rispetto al tempo qui appare la frequenza, omega può essere l'ordine del megahertz, quindi viene in ballo un fattore di 10 alla 6 e capite che il 10 alla meno 6 di alfa viene compensato dalla 6 della frequenza, per cui nei nostri memori si è questo termine molto importante. queste due equazioni ovviamente sono completamente accoppiate tra di loro è un fully coupled quindi il fully coupled termine elasticity a differenza di prima in cui prima questo termine non lo prendevamo in considerazione perché? perché non stavamo considerando dei meno oscillanti abbiamo detto che l'evoluzione era molto lenta il package che si sta riscaldando non è un package si riscalda per l'effetto R e quindi si deforma è chiaro che però la sua frequenza non sta vibrando il package è il MEMS che vibra, il package è FEVER per cui in quel caso non era importante questo termine, ma se invece vogliamo studiare proprio il MEMS che vibra questo è molto importante dopo stabiliamo come se noi siamo capaci di studiare gli autovalori di questo problema accoppiato e poi se facciamo la parte, come vedete qui, la parte reale di omega quadro diviso la parte immaginaria di omega quadro, otteniamo la stima di Q. Dopo impariamo come facciamo, ma immaginiamo di averlo fatto. Dove è importante questa dissipazione? In tutti i NEMS che hanno dissipazione a flessione. Vedete, se avete un cantilever, cioè un oggetto che è incastrato qui ed è libero, questa è la curva, è un paragone. Questi sono sperimentali pallini bianchi i pallini bianchi sono stati sperimentali e questa linea e la teoria la formula di zener che trovate spiegata nella dispensa non abbiamo tempo di farla però istruttivo bisogna ritornare non lo facciamo soprattutto perché bisogna ripassare la teoria del travi prima perché la teoria di zener e si basa sulla teoria del travi non ne abbiamo mai parlato io... ho deciso che... finito. Con l'altro. Semplicemente si è messo il tau. Ok. Quindi abbiamo questa... per il cantilever abbiamo questa soluzione di riferimento che dipende dalla densità caustica termica modulo di Yang alfa T0 omega, che è la frequenza di vibrazione questo tau z che qui non ho scritto ma dispensa dato è una combinazione di alcoli per procenti quindi quella la linea è questa formula e questi diamanti neri sono la simulazione di elementi finti che poi vedremo poi vedete che gli esperimenti la soluzione analitica soluzione numerica matchano molto bene effettivamente per questo cantilever la dissipazione termoelastica è la sorgente di dissipazione evidentemente non ci sono altre sorgenti di dissipazione sensibili e questi sono i parametri utilizzati in questo esempio per la prima volta utilizzeremo i parametri fisici non fittici ecco attenzione voi direte va bene ma il cantilever che serve? questo è un risonatore di 7-8 anni fa fatto da Seatime a 500 kHz puramente capacitivo e questo è un reverse engineering di ST è più facile capire come funziona qui c'è l'ancoraggio e ci sono questi quattro braccetti che vibrano questo si chiama double tuning fork come un diapason doppio e voi sapete perché il diapason è a forma di diapason? perché non dissipa perché i due braccietti inducono un movimento della base ma i due movimenti opposti fanno sì che alla radice lo spostamento trasmesso sia veramente trascurabile quindi quello che dicevo prima la dissipazione di ancoraggio è praticamente trascurabile perché non c'è energia elastica che migra nel sultrato questo è un doppio tuning fork perché ce l'ha di qua e ce l'ha di qua quindi noi sappiamo che senza neanche farse un'azione che la dissipazione di ancoraggio qua è trascurabile e che l'unica sorgente è la dissipazione termoelastica in questo caso quindi abbiamo assolutamente bisogno di stimare di stimare ecco questo non lo scriverei più in questo forse il tempo è un po' ottimista qui 10-2 mili insomma è un po' esagerato comunque questo qui si comporta in questa maniera è il primo risonatore di grande successo in silicio perché fino al 2016 io ero dentro un progetto europeo e si studiava la possibilità di fare risonatori in stilicio perché il problema fondamentale era sostituire il quarzo i risonatori a 32K sono in quarzo gli orologi che negli anni 80 sono esplosi erano tutti risonatori in quarzo il quarzo è stabilissimo in temperatura cioè anche se varia la temperatura tanto il quarzo risona sempre alla stessa frequenza cioè varia sì ma pochissimo il silicio invece varia molto di più e si cercava di capire come effettivamente questa vedete è una zoom del risonatore il silicio invece dritta molto in temperatura quindi solo allora si era scoperto che se si dopa fosforo tantissimo 10 alla 20 cioè praticamente c'è più fosforo di silicio 10 alla 20 diventa molto meno sensibile la temperatura. Necchè attualmente anche l'ESETI sta producendo risonatori con un altissimo doping, cioè un silicio monopristallino dopato con solo 10 alla 20, che fa anche di più. Questa è una cosa, però il fatto è che la dissipazione termolastica, studiavamo la dissipazione termolastica di questo aggetto e avevamo capito che effettivamente inserendo questi slot queste fessure in maniera opportuna blocchi la trasmissione del calore cioè impedisce che il calore fluisca in questa direzione e la dissipazione termoelastica è legato a quello e quindi riesce ad aumentare molto il fattore di qualità legato alla dissipazione termoelastica quindi si può giocare un po' sulla geometria per diminuire la dissipazione termoelastica e quindi andiamo sull'ultimo un tool di ottimizzazione per arrivare a capire che effettivamente la geometria prodotta da Seatime era una geometria ottimale per questo. Quindi ci sono veramente tantissimi MEMS che lavorano, che hanno bisogno che sia nota con estrema coerentezza la distribuzione del terramatico. Questo per motivarvi a studiarla. è anche vero che non sempre coincide qui è un esempio molto simile, un cantilever ma a differenza di prima quello di prima era spesso 30 micron questi qui invece sono spesti 2,3, 1,2 micron allora vedete questa è la curva di condizione di Zener per il cantilever che è valida per qualunque espressore, questi invece sono dei dati sperimentali, vedete che nella picco inferiore anche se non ci sono tanti dati si può intuire che arrivano a dare la stessa condizione, ma poi quando la frequenza, cioè la lunghezza del cantilever varia e si allontana da quella che dà 1, F diviso F0, c'è una fortissima deviazione dalla condizione analitica rispetto alle misure sperimentali. Poi ci sono anche dei modi di vibrare torsionali che avrebbero deformazione traccia di Y nulla quindi secondo le equazioni la dissipazione termolattica dovrebbe essere 0 e che invece predicono una dissipazione diversa e quindi ancora adesso non mi interesso più di queste cose qui ma credo che il problema sia ancora aperto cioè quando il device diventa veramente sottile ci sono degli effetti di superficie ossidazione superficiale eccetera eccetera che producono una dissipazione che a un certo punto diventa maggiore di quella termoelastica, però al momento che io sappia nessuno è veramente in grado di stimare numericamente, di predire non solo numericamente ma anche con formulette semplificate questa dissipazione di superficie. Quindi sappiate che non tutto è noto nel mondo dei MAPS, tantissime cose sono state capite negli ultimi 10 anni, 15 anni, ma molte altre rimangono. proprio di questi effetti di superficie, l'ossidazione superficiale, come impatta, come predirla, nessuno lo sa. Per fortuna i MEMS sono ancora un po' più grandi, sono più spessi di un micron, quindi siamo lontani da questo limite. Ok, adesso dobbiamo cercare di stimare la dissipazione termoelastica. queste sono le 2.13 e le 2.14 non so perché si chiamano 2 perché ho compilato solo con quei due capitoli lì in realtà sono i capitoli 5 questo è il capitolo 6 quindi dovrò rifare questi numeri sono sbagliati il numero dell'equazione è giusto ma il capitolo è sbagliato e dobbiamo implementare queste equazioni qua non è tremendamente difficile però non ci pensiamo nemmeno vogliamo ridurre al massimo le cose nuove che dobbiamo fare e riutilizzare tutto quello che siamo capaci di fare prima di tutto voglio convincervi del fatto che anche se fossimo pazienti e volessimo implementare questa equazione se poi calcolare i modi propri per utilizzare omega quadro parteriali sarebbe difficile si prende la forma forte ormai siamo abituatissimi a calcolare la forma debole la forma debole chiaramente può essere discretizzata e si otterrebbe due equazioni discretizzate di questo genere dove abbiamo dalla dinamica abbiamo la classica massa per l'accelerazione più k per lo spostamento queste sono le teste identiche che calcolavamo la settimana scorsa quindi matrici di massa, matrici di k se abbiamo bisogno di questi oggetti andiamo nel codice di Enemis che copiamo tal quale quello che facevamo allora poi abbiamo però la novità, questo d per la temperatura che è la discretizzazione di questo oggetto No, l'abbiamo visto prima, l'abbiamo commentato un pochettino, qui abbiamo la temperatura che finisce nella lista delle temperature nodali incognite, poi questo è il campo test, la traccia di y era data da quella bd, famosa bd, quindi l'integrale di questo viene assemblato e lo chiamiamo D nell'assemblaggio quindi questo termine qua integrato genererà una matrice D che moltiplica la lista delle temperature incognite poi passiamo all'equazione della diffusione prima abbiamo visto che questo termine è quello che genera la capacità termica per la derivata prima della temperatura è lo stesso di prima quindi non c'è motivo che cambi questo termine qua genera la matrice di conducibilità per la questa parentesi non so che cosa perché ci sia comunque la matrice di conducibilità per la temperatura e poi c'è quell'altro termine che adesso dobbiamo giustificare ok vedete che qui ho usato la stessa lettera che abbiamo qua sopra la stessa matrice però trasposta è vero che la stessa matrice trasposta analizziamo più da vicino a quei termini, a parte il coefficiente alfa che reggono 1 meno 2 è il lente che è lo stesso. Poi il t0 è stato fattorizzato, quindi questo t0 non lo considero, questo coefficiente non lo considero, devo paragonare i due termini che sono l'integrale di traccia di epsilon valutato sul campo test U moltiplicato per la temperatura, questo è un integrale e poi avrò l'altro integrale che è la parte del segno è la derivata della traccia di U rispetto al tempo per la temperatura la derivata della traccia di E si è andata sulla soluzione U rispetto al tempo per la temperatura t, che è la temperatura testa perché questi giudicano la stessa matrice trasposta? proviamo a pensare al generico elemento per cui facciamo l'integrale sull'elemento numero E questo l'abbiamo sviluppato prima che cos'era questo? era l'integrale questo era l'abbiamo sviluppato così che è per l'integrale di quella matrice bv trasposta, vi ricordate? la matrice bv che era la somma che era la somma delle due righe di b la somma delle due righe di b per la temperatura quindi questo T e trasposto genera questa traccia e poi abbiamo la temperatura che è n integrato di omega integrato sull'elemento ok e poi fuori qui c'è la lista di temperatura T scusate ho sbagliato qui non è T ma è è la U e Ilde andate a vedere le espressioni di prima c'era la U e Ilde perché questo è di V per U e dalla traccia adesso guardiamo qua sotto allora innanzitutto la traccia di Epsilon e sarà ancora la traccia di y valutata su u e ancora la stessa matrice bv per la lista ue. Se valeva prima per il campo test, vale la stessa identica relazione per il campo soluzione. Però qui abbiamo la derivata rispetto al tempo. Ma la derivata rispetto al tempo nella semidescredizzazione degli elementi finiti finisce sul valore modale. Questa matrice non cambia. e poi abbiamo la temperatura T tilde ma la temperatura T tilde è ancora una volta la vista N per la lista delle temperature T tilde quindi sostituiamo all'interno questa cosa però il capotest è TE quindi scrivo dapprima la temperatura T tilde faccio il trasporto e ho TE tilde che moltiplica l'integrale sull'elemento E di N trasposto perché ho scritto prima E e per poter scrivere T e trasposto fuori ho fatto il trasposto di questo oggetto e poi devo scrivere la traccia di U che è uguale a BD poi lo integro sul volume e lo moltiplico a destra per U punto ecco, allora paragonate queste due matrici e trovate che sono trasposte una dell'altra ed è per questo che questo è vero a livello elementare cioè su un generico elemento T6 o qualunque elemento isoparaventrico le due matrici sono trasposte una dell'altra rimane vero anche dopo aver fatto l'assemblaggio per cui queste due matrici sono una alla trasposta dell'altra ok? questo per spiegare questa notazione allora vedete che è un sistema omogeneo perché non ci sono forze non stiamo forzando con niente perché vogliamo risolvere il problema agli autovalori il problema agli autovalori vuol dire che assumiamo che gli spostamenti hanno questa forma e le temperature hanno questa forma e inseriamo questa ipotesi qua dentro allora il problema è che se inseriamo questo guadagno avremo un omega quadro e se invece inseriamo questo qui avremo un omega quindi è un problema molto diverso da quello che eravamo abituati nella dinamica in cui avevamo solo omega quadro qui omega quadro, qui omega e qui omega quindi diventa quello che viene chiamato il problema agli ottovalori quadratico perché contiene tutte le potenze di omega omega quadro e omega ed è difficilissimo da risolvere mentre il problema classico ha degli algoritmi efficientissimi per risolverlo questo problema è straordinariamente difficile e l'unico modo per risolvere in maniera efficiente e ridurre tutto al primo ordine ridurre tutta questa equazione al primo ordine cioè aggiungere un'altra equazione che vi dice ok chiamiamo v la derivata prima di u e questa è un'equazione aggiunta e poi questo qui allora sarà ad esempio diventerà m per v punto quindi vedete che trasformo tutto e così via trasformo tutto in un sistema di equazioni differenziali del primo ordine solo che lo espando diventa un sistema molto più grosso perché sto introducendo anche tutte queste relazioni che prima non c'erano e tipicamente la dimensione di questi problemi diventa talmente grande da essere proibitiva e poi quello che non va bene è che comunque le matrici finali che ottengo non hanno le belle proprietà di essere simmetriche definite positive come era in dinamica quindi comunque lo si metta risolvere questo problema gli autore valori quadratico è una grossa sfida e nei codici alimenti finiti che lo fanno è facile mettere in crisi perché basta nelle nostre strutture che sono sottili lunghe che hanno insomma tante peculiarità frequenze alte risolvere il problema lì molto spesso genera solo delle genera solo delle grandi forcherie e quindi bisogna trovare una abbiamo sviluppato un po di tempo fa una soluzione alternativa che super efficiente e vale solo esclusivamente se la dissipazione bassa come sempre però nei nostri nei nostri mess a parte gli accelerometri ma gli accelerometri sono talmente facili da analizzare che non ci sono paranoie di questo genere questa sarebbe la procedura globale però c'è una procedura molto semplificata che ci obbliga a mettere assieme un po' tutti i tasselli che abbiamo visto settimana scorsa e oggi prima assunzione la dissipazione termoelastica è così piccola che il modo meccanico sarà modificato pochissimo dalla presenza dell'accoppiamento termoelastico. Quindi io affermo che questa U, il mio modo meccanico, sarà ancora proporzionale a questa quantità UB, che io chiamo UB, che è il modo calcolato col vecchio codice di dinamica modale. cioè io calcolo più in modo del cantilever e assumo che il movimento meccanico sarà ancora proporzionale a questo perché la dissipazione meccanica la dissipazione termoelastica è così bassa che non è in grado di modificarlo l'unica cosa che modifica la dissipazione termoelastica è l'aggiunta di una parte immaginaria alla sequenza ma quella la calcoliamo in maniera diversa quindi passo 1 che adesso dobbiamo fare bisogna andare a prendere il nostro cantilever e calcolare il modo che ci interessa con il vecchio codice lo facciamo perché è cambiato un pochettino e vi dico in che cosa è cambiato quindi andiamo dovrebbe funzionare senza i vostri problemi apriamo l'editor facciamo edit model punto m e me lo apre qui il vecchio codice si chiamava model è rimasto tal quale facciamolo girare sul file input che si chiama modal cantilever modal cantilever fa l'analisi e dovrebbe far vedere i modi solo che siccome i parametri qui sono reali il movimento è così piccolo che devo andare a amplificare molto la visualizzazione cambiando sulla view del primo modo cambio il fattore di amplificazione ma devo cambiarlo a meno di un fattore mi ricordo 10.000 ecco, 10.000 è troppo, 1000 ecco, quindi si è calcolato tutti i modi e il primo modo è il classico modo di flessione di un cantilever vi è andato? ok allora adesso lo chiudiamo però in realtà il codice ha fatto una cosa alla fine che la volta scorsa non faceva si è salvato esattamente come prima salvavamo la soluzione della diffusione termica così adesso dobbiamo salvare il modo il modo di vibrare in questo caso il primo modo a frequenza più bassa per poterlo poi riutilizzare successivamente. Allora vedete qui alla fine, nelle ultime linee, definisco qual è il modo che voglio salvare, crea questa matrice che ha la lunghezza dei nodi e due componenti perché ha il componente x e y, fa il loop su tutti i nodi per tutte e due le direzioni, riempie ub con il valore di V, di grande l'output della routine di matlab che calcola gli autovettori quindi di grande contiene tutti quanti gli autovettori prendiamo l'autovettore che ci interessa in questo caso il primo e poi riempiamo quindi ub con tutte le componenti nodali del modulo quindi ub è una matrice con tante righe quanti nodi e due colonne componente x componente y dell'autovettore ok? vi ricordate quali erano le proprietà degli autolettori nella meccanica e ci serviranno erano orto erano normalizzate in massa quindi se io prendo un B per farci un B trasposto M un B che cosa ottengo? lo ricordate? che cosa ottengo? invece se faccio uv trasforzo k uv che cosa ottengo? va bene lo ricordate? m ortonormali e k ortogonali? no allora questo è giusto qui viene 1 e qui viene omega 0 cioè normalizzato in massa poi sono normali anche rispetto a K cioè se prendo due autovettori diversi viene 0 se invece uso lo stesso autovettore mi viene la frequenza proprio al quadrato quindi queste due relazioni sono molto importanti le useremo tra poco ecco, tutti i valori nodali di questo autovettore sono stati salvati in questa matrice che poi li scrive e poi li scrive in questo file che quindi dovreste trovare a questo punto nella vostra directory e salva due cose, il valore di ω0, che è la frequenza associata all'autovettore, e UB, che contiene tutti i valori nodali dell'autovettore. Quindi è il vecchio codice, puramente meccanico, senza nessun effetto termolastico, ha calcolato il modulo e l'ha salvato. Quindi noi adesso ritorniamo alla nostra ipotesi, vediamo prima la procedura che abbiamo implementato e poi cerchiamo appunto di metterla a posto quindi abbiamo fatto il primo passaggio i quadri rossi sono le cose che bisogna fare bisogna fare tre passaggi il primo è stato fatto, abbiamo calcolato il modo meccanico ed ora e poi sottorremo sempre che l'evoluzione dei valori nodali sia proporzionale a questo modo meccanico moltiplicato per questo e alla i omega t dove omega è la frequenza che stiamo cercando omega non sarà omega 0 ma sarà una nuova frequenza che cerchiamo poi questa è l'equazione meccanica poi entriamo nell'equazione della diffusione cioè la seconda, entriamo in questa seconda equazione qua assumiamo che di conoscere questo termine e in questa assumiamo che la omega sia uguale alla omega zero perché omega è molto vicino alla omega zero a risonanza calcoleremo la variazione successivamente a partire dalla prima equazione quindi facciamo questa ipotesi e assumiamo che invece di voler calcolare la distribuzione di temperatura in questa forma cioè il prodotto di un vettore che vogliamo stimare per e alla i omega t no, ho detto una fesseria effettivamente non stiamo no, no, no no, questo è io ho fatto una risposta forse insomma pensavo di averlo cambiato sì, no, la risposta è giusta c'è anche l'omega 0 facciamo quindi lo faccio alla lavagna prendiamo l'equazione della diffusione e inseriamo quindi le nostre due ipotesi le nostre due ipotesi sono uno, che lo spostamento si comporta così ma in questa equazione assumeremo l'omega uguale a omega 0 per il comportamento della parte meccanica e poi assumiamo che la temperatura abbia questa forma il prodotto di un vettore che dobbiamo determinare per l'esponenziale complesso. Quindi se lo inseriamo là dentro, otteniamo chiaramente iω per il vettore incognita td, anzi m, scusate, m per td. Poi il secondo termine rimane identico. tutto ovviamente verrà molto indicato per e alla i omega t ok e poi abbiamo il c per c per t b c per t b e poi metto dentro anche il termine no facciamo così questo è e alla i omega t e questo sarà uguale alla derivata il corretto quindi sarà uguale a meno a t0 non c'era il meno lì vabbè comunque non c'era il meno in una frazione non vi fate qualche pasticcio nella scrittura delle vabbè usiamo meno a t0 per i omega i omega 0 cioè qui sostituisco omega 0 e poi la provo B per E alla I omega T. Allora, questo cancello è I omega T a destra e a sinistra e quindi questo mi genera un sistema lineare in cui l'incognita ITB, risolvendo questo sistema lineare, ottengo la stima della distribuzione di temperatura all'interno della trave, quando la trave sta vibrando con il modo meccanico UB. Questa è comunque la fase 2. Dopo di questo ho ottenuto la miglior stima possibile TB, con certe approssimazioni, vedete che sto facendo alcune approssimazioni, che però si basano tutte sull'idea che la dissipazione è piccola, quindi ω5, ω0, eccetera. poi adesso prendo l'equazione meccanica che era la prima delle due che ho che è questa che non ho ancora utilizzato perché quello che ho utilizzato in realtà è semplicemente la parte puramente meccanica per calcolare in modo non l'ho mai accoppiato con l'effetto della temperatura adesso prendo questa equazione con l'effetto della temperatura con l'effetto della temperatura e sostituisco all'interno di questa equazione le due distribuzioni che abbiamo trovato. Ok? Allora facciamolo e vediamo che cosa otteniamo. Ottengo, allora abbiamo l'equazione che è m due punti, quindi meno omega quadro, meno omega quadro per m che moltiplica la u b e poi dappertutto e alla i omega t ma non lo scrivo tutti i termini sono moltiplicati per questo termine per questo esponenziale complesso e quindi non lo scrivo da nessuna parte allora poi ho k u k u b k u b e poi ho questo termine qua che è meno d per la temperatura e quindi sarà quindi questo è uguale a d per la temperatura Tb che ho calcolato prima tutti quanti sono moltiplicati per l'esponenziale complesso che quindi si semplifica quindi non siamo ancora arrivati a questa espressione ma abbiamo solo fatto una prima sostituzione all'interno di questa equazione noi vogliamo utilizzare questa equazione per stimare la nuova frequenza omega quadro e questo è un nostro obiettivo però qui non riusciamo perché abbiamo un sistema quindi bisogna passare da un sistema a un'equazione scalare e il trucco che si usa correttamente in queste cose è proiettare questa equazione equazione sul modo meccanico che ci interessa. Cosa vuol dire proiettare sul modo meccanico? Vuol dire pre-multiplicare tutti quanti questi termini per un B trasposto. Cioè vuol dire andare nel principio delle potenze virtuali e imporre il principio delle potenze virtuali non per tutte le scelte possibili del campo teste ma per una scelta molto precisa in cui il campo teste coincide con l'autovettore meccanico coincide con l'autovettore meccanico quindi se noi moltiplichiamo a sinistra il ruolo di trasposto se dopo abbiamo tempo ritorno su questo punto e adesso lo accettiamo semplicemente e basta cioè ho premoltivicato a sinistra rubì, trasposto e ottengo quindi un'equazione scalare ma un'equazione scalare che ha delle bellissime proprietà perché adesso si semplifica notevolmente e si semplifica perché l'abbiamo ricordato prima che queste quantità qua hanno dei valori precisi questa vale 1 e questa vale ω0 quadro. Per cui, mettendo assieme le cose, ottengo che omega quadro è uguale a omega zero quadro meno u b trasposto per d per t b. Questa è la relazione che state cercando. quindi la nuova frequenza omega quadro sarà data dalla somma di omega zero quadro, la vecchia frequenza omega zero quadro più questo termine qua che sarà certamente immaginario perché tb sarà una soluzione immaginaria, cioè sarà certamente complesso perché il vettore tb sarà un vettore complesso che viene dalla soluzione di un sistema complesso io dopo prenderò la parte immaginaria di questo oggetto della parte reale e con la formula che abbiamo visto prima, Q uguale alla parte reale di omega quadro diviso la parte immaginare di omega quadro, stimo il fattore di qualità. Stimo il fattore di qualità. Quindi sono dei passaggi molto semplici che richiedono di calcolare un problema agli otto valori standard, risolvere un sistema lineare e fare questo calcolo finale che però è un calcolo rapidissimo. Tutto questo è stato implementato, questa procedura è stata implementata all'interno di un codicillo che voi avete che si chiama ah no, questo l'abbiamo già visto, l'applicazione al cantilever è la prima parte, abbiamo calcolato il modo e l'abbiamo salvato, questa è la prima parte e l'abbiamo già vista. Adesso invece quello che dobbiamo fare è analizzare e poi dovremo modificare leggermente questo file che si chiama TED.m. Apriamo il file TED.m. Ecco. Cosa fa il file TED.m? Per prima cosa carica il file che abbiamo salvato prima. Quindi legge omega0 e ub. Omega0 è la frequenza propria del modo meccanico e ub era la matrice che contiene tutti i valori nodali del modo. E poi si propone di implementare questo codice. Questo codice implementa le due fasi che abbiamo visto prima utilizza come file di input questo cantilever.m il cantilever.m che voi avete se guardate ha effettivamente le proprietà del materiale realistiche il bodo di yang è 150.000 cosa alla 150 giga pascal evidentemente sta usando il pascal e il micron comunità di misura quindi invece che 150 per 10 alla 9 è diventato 150 no, non può essere pascal non può essere pascal certamente usa il micron ma non può essere pascal probabilmente è pascal di 10 alla meno 6 quindi micron newton comunità della forza il micron newton comunità della forza e il micron comunità di lunghezza allora 150 giga pascal quindi 150 per 10 alla 9 su micron dovrebbe quindi essere diviso da 10 alla 12 però siccome è micron-newton sì, probabilmente è l'unità giusta poi la densità è 2300 kg al metro cubo quindi probabilmente usa i picochili i picochili come unità di misura e sono le unità tipiche di NEMS il picochilo come massa e il micron come lunghezza e così via è un esercizio non così facile è certamente la cosa più complicata all'interno di questo capire quali sono le unità che sono state utilizzate per assegnare Yang e Poisson è 0 densità questa è la C la capacità termica che unità è stata usata per la capacità termica la K e l'alfa, beh quell'alfa è facile perché tanto è addimensionale quindi in questo caso è sempre 10 alla meno 6 tipica del silicio quindi questo è effettivamente un materiale realistico poi solido è definito come sempre e gli spostamenti sono imposti a sinistra del cantilever li blocca, quindi il cantilever è effettivamente incastrato se no non si chiamerebbe cantilever e poi definisce la temperatura T0 che è la temperatura iniziale che è la 298, quindi 25° insomma siamo a 25° temperatura ambiente, quindi niente di che dunque legge il file legge queste informazioni legge anche questo vettore UB come abbiamo visto prima e comincia ad implementare questo, deve arrivare a risolvere questo sistema deve arrivare a risolvere questo sistema e quindi assembla la matrice di capacità e la matrice di conducibilità esattamente come nei problemi di diffusione però deve aggiungere questo oggetto qua deve aggiungere questo integrale che oggi abbiamo già visto tante volte per fare questo integrale adesso se provate a far andare il codice vi dovrebbe dare errore vi dovrebbe impiantarsi a un certo punto se sono sbagliato a caricarmi io ho caricato quello giusto gli dai errore andiamo a vedere perché gli dai errore alla linea 78 alla linea 78 dai errore alla linea 78 e dove, no, com'è alla linea 78 come è possibile no, allora ho sbagliato un'altra cosa il nostro errore forse è giusto perché il nostro errore è giusto perché reindirizza al completamento delle forze nodali vediamo un secondo come mai ma forse ho dato solo ho scelto male il file di input sì, eccolo dice che la linea 117 questa funzione che calcola che calcola l'integrale di prima che calcola per un elemento questo integrale non è completa certo perché ho levato alcune righe non ho legato qualche riga allora andiamo a vedere effettivamente se è così l'apro T6 tende FE eh sì, mi manca vedete che ho messo puntini puntini puntini puntini cioè quello che chiedo non so se riuscite a farlo però almeno una decina di minuti cercare di sbattere la testa per capire cosa bisogna fare lì dentro sapendo che dei termini simili ci sono nella termoelasticità che abbiamo visto prima quindi sapendo che tutta una parte di cose è stata fatta, qui mancano un paio di righe non mancano tante righe forse una, adesso non mi ricordo bene dovete riempire quelle righe lì per riuscire a calcolare questo integrale questo integrale attenzione voi t0 è esterno vedete che qui la funzione viene chiamata con la x che è la geometria modulo di Y con il vicente di ragione termica e Y che cos'è we? we è la lista degli spostamenti nodali per l'elemento quindi sui sei nodi dell'elemento le due colonne del modo del modo del modo meccanico perché ricordatevi che io sto implementando questa equazione l'equazione termica supponendo che lo spostamento sia quello del modo meccanico quindi devo implementare questa questa quantità supponendo che U sia noto e che sia in modo meccanico quindi se U è noto ed è in modo meccanico come faccio a implementare questo oggetto all'interno di questa routine integrazione numerica Gauss è già tutto inserito devo solamente dire come calcolo T e come calcolo questo oggetto non l'abbiamo detto prima provate a insistere una decina di minuti per capire come si deve completare come si devono completare queste linee ricordatevi che la derivata rispetto al tempo è come sono la derivata rispetto al tempo diventa questa I omega 0 perché stiamo assumendo questo stiamo assumendo che lo spostamento sia proporzionale al modo moltiplicato per l'esponenziale complesso quindi quando fate la derivata viene io omega 0 quindi questa derivata non vi deve preoccupare l'unica cosa che dovete fare è la traccia di y valutata su questo su questo b Grazie. vedete pure tutto quello che volete allora più o meno scoppiazzando l'estra manca dovreste riuscire a non scrivere neanche una riga però bisogna sapere bene dove andare a copiare cercare di rispondere a tutte le domande ti dite come lo scrivete sono molte cose che ho scritto prima quando parlavo del termo elastico ti dite come lo scrivete n che moltiplica t e t tilde n moltiplicato per t tilde e attenzione che quando fate quell'integrale su un elemento ok quando provo a riscrivervi vi riscrivo esattamente quello che avevo scritto prima quindi togliendo t0 alfa e diviso 1 lascio solamente scrivo t tilde per primo e poi la traccia la derivata rispetto al tempo della traccia di epsilon valutato su u che però u è in questo caso il campo uv questo è il noto e questo deve diventare qualcosa scritto così t e tilde trasposto per l'integrale di qualcosa e questo è il contributo che state cercando il contributo delle forze nodali che cercate perché questo viene noto quindi voi dovete fare il prodotto di trasposto per Perfetto, andiamo. Ma perché se è carico? Per su, così. Quindi questo è la prima cosa. Errore dimensionale... e poi un'altra cosa di cui avete bisogno è la seguente magari non ve lo ricordate più ma qui avete tanti elementi che vi servono questa è la cos'era D, la derivata delle funzioni di forma rispetto alle coordinate e poi però per calcolare il gradiente di U vi serve la matrice del gradiente la matrice G come si calcolava la matrice G? questo è stato fatto in tantissime routine G, eccola. Nocchi. e adesso il nostro primo obiettivo è calcolare la matrice G che abbia come componente GJ le derivate di MI rispetto a XJ perché poi questo vi serve per calcolare il gradiente del campo 1 però voi non avete questo avete la derivata di n rispetto ad a che è in d questo lo trovate nel capitolo 3 oppure nelle routine k6, kt6 qualcosa del genere derivata composta questa è DIM e questo è dove potete andare a cercare la derivata di A rispetto a X è l'inverso della derivata di X rispetto ad A ok ok, trasposto come faccio a calcolare questa g ho tutto ho trovato su una video in precedente sì assolutamente allora se volete il gradino interno cominciamo a scriverlo cominciamo a dire i elementi di cui abbiamo bisogno allora innanzitutto dovremo calcolare la matrice gradiente che è quindi g uguale a d uguale a d per inv j inv inv j quindi calcolerà g che è la derivata di ni rispetto a xj quello lo uso per calcolare il gradiente di lo devo usare per calcolare il gradiente di U andiamo un po' insieme come si calcola il gradiente di U allora il gradiente di U è la derivata è quella matrice che ha la derivata di Ui facciamo così la derivata di UM rispetto a XK uso degli indici diversi perché se no ci confondiamo comunque questa è la matrice del gradiente se io riesco a calcolare questo poi la traccia è facile perché la traccia sarà la somma della derivata di 1 rispetto a x1 più la derivata di 2 rispetto a x2 quindi io prima voglio calcolare questo ma u è interpolato quindi è la sommatoria di n i moltiplicato per le componenti nodali quindi i nella componente m è questa cioè questa è la componente um è la sommatoria su tutte le funzioni su tutti i nodi della funzione di forma per i valori nodali e poi devo fare la derivata allora avrò la derivata di ni rispetto a xk questo è il gradiente che ho appena calcolato e questo lo trovo nella matrice che contiene i valori nodali nella matrice Ue che contiene i valori notari questa qua Ue è la matrice che ha 6 righe e 2 colonne quindi quella è Ue Ue tolgo la E per non confondermi con gli indici i nodi sono la riga I e la la colonna e la componente N e devo riuscire a fare questo come faccio a scrivere questo per tutte le componenti ho la G ho la U e questo è il gradiente di U che quindi io chiamo GU la matrice che ha proprio quelle componenti lì non può essere G per U perché l'indice I è il primo per tutte e due quindi devo fare la trasposta di una quindi devo fare g trasposto per ue quindi la traccia poi è facile calcolarla se ho quella matrice perché la traccia la chiamo treps sarà data da che cosa? dalla somma di gu 1 1 più g g 2 2 questa è la traccia siamo pronti dovete solo compilare ultimanti marino mali 20 elementi se lo comprate e lo fate andare con me dovrebbe aprirsi Gmesh con la distribuzione di temperatura parte reale e parte immaginaria adesso taglio di trasposto per 1 e abbiamo questo abbiamo n andiamo a frequenziare sicuro per frequenziare ecco se ti riusciti a compilare questo forse aspetta fa vedere ti aspetta quello lì aspetta la temperatura e quella lì si perché la temperatura fluisce dalle parti però io sono spallatissimo no infatti è spallatissimo nel bending provate a pensare al bending c'è una parte compressa che la parte se la pala stava su la parte compressa e su si sta tirando verso il super si mette che dipende l'autovettore mentre va verso l'alto c'è la parte compressa su quindi che la parte conversa genera una sorgente di calore mentre quella che si estende un sink di calore no per cui poi calore e fluisce in una certa direzione fluisce in direzione toglie alla trave per questo che la variazione di temperatura è ha un gradiente anche lungo x però l'ingrediente più importante è lungo y e questo quindi è se andrete a studiare la soluzione di Zener però vediamo di correggere insieme questa prima parte cosa manca? per integrare questa cosa qui bisogna mettere l trasposto e poi dobbiamo mettere la traccia quindi la traccia io l'ho scritta qui treps uguale a g1 quindi dovrò mettere treps treps e quella è la traccia e poi cosa dovrò mettere? dovrò mettere lo jacobiano per integrare e il peso il peso di WG per completare l'integrazione e qua non bisogna dimenticare i parametri che io non ho più scritto che sono e diviso alfa per e diviso 1 quindi avremo n trasposto copio il coefficiente che era alfa per e diviso 1 meno 2 per 1 questo era il coefficiente e poi devo moltiplicare per 3ps la traccia di epsilon che è moltiplicato per lo jacobiano che si chiama determinante di j per il peso di gauss che si chiama w gauss ovviamente del punto di gauss corrente che è il peso di Gauss G perché sto facendo il loop su G se questo è giusto se questo è giusto per me è facile perché sto copiando è giusto e lo faccio andare ecco fino a qui arriva c'è questa distribuzione di temperatura e la distribuzione di temperatura che è una distribuzione con un gradiente molto debole in direzione x ma un gradiente molto forte in direzione y perché le parti compresse e le parti tese sono sopra e sotto, quindi il flusso di calore avviene in questa direzione. Poi c'è la parte reale e la parte immaginaria, vedete che la parte reale, qui sotto avete le due scale, la parte reale è piccolina 10 alla meno 5 la parte immaginaria è pur sempre piccola ma ha un ordine di grandezza superiore 10 alla meno 4 quindi la componente principale è la parte immaginaria se voi chiudete Gmesh lui va avanti e vi dice ancora ciao, no, devi finire un altro calcolo andiamo a vedere che cosa rimane da fare l'ultimo sforzo abbiamo finito vi ricordo il riassunto delle fasi le fasi da fare non ho tre le fasi da fare queste tre 1 calcolo del modo in questo l'avevamo fatto 2 il calcolo il calcolo di queste di questo sistema e questo è quello che abbiamo appena fatto solo che qui manca un A omega 0, l'ho detto prima, che devo aggiungere nella slide, ma nella dispensa giusta, e poi adesso dobbiamo fare l'ultimo, l'ultima fase, questa. Cioè noi dobbiamo calcolare questo prodotto qua, un B trasposto di TB, non lo faremo in questa maniera formale, cioè calcolandoci la matrice e poi facendo il prodotto U trasposto per TB, ma lo calcoleremo direttamente all'interno della routine, cioè cos'è che facciamo? Dobbiamo andare a fare questo integrale qua, perché da dove deriva questo prodotto? Deriva da questo integrale, cioè sempre lo stesso coefficiente, alfa e diviso 1 meno nu, che moltiplica la traccia di y valutata sul campo test, ma qui forse è un pochino più chiaro, io prendo in questo caso il campo test come modo meccanico cioè non lo impongo per un generico campo test non è vero che lo impongo per un generico campo test ma lo impongo per la scelta di campo test utile e uguale a modo meccanico a unito quindi qua dentro io dovrò prendere il modo meccanico dovrò calcolare la traccia di y sul modo meccanico e caspita, l'abbiamo appena fatto quindi bisogna copiare questo e poi bisogna però valutare la temperatura la temperatura TH è una temperatura variabile da punto a punto quindi dobbiamo interpolare anche la temperatura però questo è un po' più semplice ormai di prima secondo me riuscite a modificarlo voi andiamo a aprire questa funzione che cosa manca all'interno di questa funzione e se hai fatto questo hai fatto tutto che miseria no, questo è il vero vecchio, scusate questo, cut head manca anche qui una piccola parte però alcune cose sono molto simili a quelle che avete appena scritto il calcolo di Treps è lo stesso esattamente scusate il calcolo di Treps è esattamente lo stesso quindi fate un ultimo sforzo un ultimo sforzo e vi manca solo da valutare la temperatura da interpolare la temperatura per conto che vedete è che entrate in questa funzione con non solamente il campo di spostamento su tutti i sei nodi ma anche il campo di temperatura calcolato prima quindi avete 6 per 2 Ue 6 per 1 Tie perché Tie è un campo scalare il resto è lo stesso la geometria e i nuovi ... ... ... ... ... ... ... ... ... ... ... ... ... ... se lo combinate bene il codice va avanti e paragona la vostra soluzione con la soluzione analitica e vi scrive alla fine nell'ultima riga del codice trovate già la formula di Zener implementata lui scrive il Q stimato e il Q analitico Grazie. Grazie. Grazie. si potrebbe rimettere la formula all'integrale quella di prima? sì no, quella sulle slide sì ok anche se ho imbrogliato un po' perché era sul vecchio codice c'era? una volta scorsa, sì ah, che miseria! però non è leggendo roba e quindi quanto ti viene? i due Q quanto ti vengono? uno viene il Q-analytic viene 2.3 per 10 alla quarta e l'altro non ho capito, vi do 1 per 10 alla quarta per una matrice eh no, allora forse non era ancora completo un consiglio è che quello che manca lì è praticamente per alcune righe identico a quello che avete appena scritto nell'altra funzione no decisamente no no 2 per 10 era 4 2,3 per 10 alla 4 no ma quello analitico quantomeno dovrebbe calcolarlo giusto no è giusto quello sì l'altro no eh no assolutamente no allora vediamo qual è quello sbagliato il peso c'è allora noi dobbiamo fare due cose allora dobbiamo calcolare la traccia di epsilon valutata su b che è esattamente la stessa cosa di prima e poi th andiamo a fare allora così visto che è la stessa cosa di prima io copio la parte le tre righe che abbiamo scritto prima in cui calcoliamo la traccia questa qui no? queste tre righe le preso e le copio tal quali nel nuovo lo metto qui via però mi manca vedete qua la temperatura però i valori nodali e abbiamo detto che la temperatura è n per t n per t per cui t la temperatura siccome nella lista che c'è qua sopra la lista delle c'è quindi nt uguale a n per t e dobbiamo solo fare l'integrale adesso no però l'integrale si fa sempre la stessa maniera perché sarà dato da quelli coefficienti ci sono i coefficienti elastici che sono quindi alfa però questi li vado a copiare anche loro da prima non voglio sbagliarmi alfa per e diviso nu è sempre lo stesso coefficiente quindi lo copio dalla routine di prima bene poi che cosa devo aggiungere? devo aggiungere la traccia di eps che moltiplica traccia di eps che moltiplica la temperatura che moltiplica che cosa ho combinato? che moltiplica il peso di Gauss quindi sempre lui che moltiplica lo jacobiano scusate, det j det j che moltiplica il peso di Gauss quindi v Gauss g allora questo è giusto vediamo un po' se funziona se no non sapevo se funziona mi chiede sempre lui termo no ok flotta questo bam si funziona 239 ok quindi la predizione analitica dice che è 23.300 e questa semplice procedura un po' approssimata ma mica tanto dice che è 23.600 allora è la produzione analitica che sbaglia o è la nostra produzione approssimata io sono convinto che è molto più vicino la verità alla nostra soluzione che la produzione analitica la produzione analitica fa un sacco di se andrete a vedere fa un sacco di ipotesi se sei, se sei, se sei, se sei eccetera eccetera noi ne abbiamo fatte effettivamente molto molto meno con un fattore di qualità così alto quindi nu è piccolissimo e tutte le ipotesi che abbiamo fatto di mischiare una volta un po' la omega con la omega 0 un po' di queste porcherie sono assolutamente giustificate perché facciamo errori almeno del nu quadro quindi nell'ordine dell'inverso di questo al quadrato quindi nell'ordine di 10 alla meno 8 quindi assolutamente completamente trascurato quindi se volete se volete giocare con questo ma immagino che non sia la vostra la vostra priorità, ritornate a questo diagramma e vi costruite questo diagramma per varie lunghezze delle travi. Noi siamo in una trave, che lunghezza ha la nostra trave? Sapete che non me lo ricordo, bisogna andare a vedere un angeo, quella che stiamo cercando di simulare, dove è il cantilever modal cantilever che poi è esattamente la stessa 1000 1000 e altezza 30, quindi siamo qua più o meno, siamo in questo punto qua abbiamo simulato questo qui ci sta in scala logaritmica 23.000 ci sta ok per oggi direi che vi ho torturato abbastanza.