Teoria 14

Da CdM_unimore.

PROBLEMA VINCOLAMENTO MATRICE BANDATA

Una volta ottenuta la matrice di rigidezza del sistema è necessario imporre le forze applicate ai nodi, prima considerando le forze esterne e, conseguentemente, le reazioni vincolari. Per queste ultime però non si è a conoscenza della reazione, ma del fatto che sia imposto un certo spostamento, in particolare nullo. Nell'analisi agli elementi finiti agli spostamenti questi risultano incogniti e quindi non è possibile vincolarli, è piuttosto necessario andare ad imporre una certa forza che tramite un'equazione di facile risoluzione, corrispondentemente al nodo interessato, fornisca un risultato intuitivo per lo spostamento.


Per la matrice intera è necessario annullare la riga corrispondente al nodo di interesse nella matrice di rigidezza del sistema, fatta eccezione per la componente diagonale, quindi è necessario annullare la colonna ancora corrispondente a quel nodo, ma non prima di aver portato a termine noto gli elementi qui presenti. Per come viene costruita la matrice di rigidezza del sistema in forma bandata questo procedimento risulta più complesso, in quanto gli elementi appartenenti ad una data riga della matrice piena sono disposti nella matrice bandata sia sulla riga che sulla diagonale e lo stesso per gli elementi sulla colonna della matrice piena, in particolare:

  • gli elementi sulla riga della bandata sono gli elementi della parte sopradiagonale della riga e della parte sottodiagonale della colonna della matrice piena,

la quale è simmetrica;

  • gli elementi sulla diagonale della bandata sono gli elementi della parte sottodiagonale della riga e della parte sopradiagonale della colonna della matrice

piena (che è simmetrica).


Si capisce quindi come sia impossibile annullare le componenti corrispondenti alla riga senza annullare quelle corrispondenti alla colonna. Perciò è necessario portare gli elementi della colonna al termine noto. Questo non è comunque l'unico problema, infatti per le particolarità della matrice bandata non è detto che sia possibile annullare tutti gli elementi della riga o della diagonale corrispondenti ad un certo g.d.l. per un dato nodo. In particolare è possibile che per una data riga, disposta troppo in alto, la diagonale incontri il limite superiore della matrice prima di raggiungere l'ultima colonna, o, al contrario, per una riga disposta troppo in basso la riga incontri la zona in basso a sinistra della matrice, riempita con gli zeri senza senso fisico, prima di raggiungere l'ultima colonna. Per questo motivo è necessario imporre un criterio di arresto in questi casi, come si ha nell'algoritmo sotto descritto, e che è intuibile dallo schema per una generica matrice bandata mostrato in Figura 1 e 2:

Figura1


Figura2

Uno schema analogo è utile per intuire in che modo va modificato il vettore dei termini noti, così come si vede nella subroutine sotto riportata. Per farlo è necessario andare a sottrarre al vettore dei termini noti pre-vincolamento un vettore contenente al termine corrispondente al g.d.l. del nodo da vincolare il valore dello spostamento imposto da vincolo, con parte superiore contenente gli elementi nella diagonale e con parte inferiore contenente i membri della riga della matrice bandata, come si vede chiaramente dallo schema qui riportato:

Figura3

Subroutine di vincolamento per matrice bandata

c234567================================================================!      
      Subroutine vincola (nsize,mband,aband,b,idof,val)
      real aband(nsize,mband), b(nsize)
c     Dove nsize indica il numero delle righe della matrice bandata, 
c     mband è il numero delle colonne e quindi della larghezza della 
c     banda, aband è la matrice di rigidezza, b è il vettore dei 
c     termini noti, idof è il grado di libertà specifico da vincolare e 
c     val è il valore a cui vincolare quel gradi di libertà.
c     Correzione termine diagonale a(idof,idof)=aband(idof,1) e b(idof)
       aband(idof,1)= 1.0       !viene così portato a valore unitario
       b(idof)= val
c     spostamento termini sopradiagonali di a(*,idof) a termine noto e
c     annullamento degli stessi.
c     L'asterisco indica "per ogni riga"
      do isd=1,min(mband-1,idof-1)   !evito di sforare la banda e il
                                     !limite superiore della matrice
        b(idof-isd)= b(idof-isd) - aband(idof-isd,1+isd) * val
        aband(idof-isd,1+isd)= 0.0  !così facendo viene annullato 
      enddo                         !il termine
c     spostamento termini sottodiagonali di a(*,idof) a termine noto e 
c     annullamento degli stessi
       do isd= 1,min(mband-1,nsize-idof)
         b(idof+isd)= b(idof+isd) - aband(idof,1+isd) * val
        aband(idof,1+isd)= 0.0
      enddo
      end subroutine


Subroutine.jpg


Analizziamo gli elementi della subroutine:

  • Nsize: sono le righe della matrice di rigidezza bandata (numero dei gradi di libertà del problema)
  • Mband: sono le colonne della matrice di rigidezza bandata (larghezza di banda)
  • Aband: è la matrice di rigidezza in forma bandata
  • b: è il vettore dei termini noti
  • idof: è l'elemento (grado di libertà) della matrice che scelgo di vincolare
  • val: è il valore che assegno a quell'elemento, cioè gli assegno il corrispondente valore di δi (delta segnato i)

Dichiariamo aband come una matrice con nsize righe e mband colonne e dichiariamo b come un vettore di nsize elementi. Volendo, al posto di mband, potrei mettere *, in quanto è necessario dare tutte le dimensioni, tranne l'ultima, che può essere lunga a piacere, a seconda di dove si va a leggere.

Successivamente effettuiamo la correzione del termine diagonale:

  • Aband (idof,1) = 1.0 --> qui sto assegnando un valore alla riga idof, colonna 1, della matrice, cioè al termine diagonale
  • b(idof) = val --> qui sto assegnando un valore al termine idof-simo del termine noto

Se modifichiamo subito i termini diagonali, è buona cosa che non li ritocchiamo più in seguito, cioè se facciamo qualcosa all'inizio che porta un oggetto al suo stato definitivo, dobbiamo stare attenti a non modificarlo ulteriormente.

Vediamo ora lo spostamento dei termini sopra-diagonali di A(*,idof) a termine noto ed annullamento degli stessi:

  • innanzi tutto quando dico i termini sopra-diagonali di A(*,idof), nelle parentesi di A indico con:
  * "*"   --> ogni riga della matrice A
  * idof  --> solo quel valore idof di quella colonna

Successivamente:

  • isd=1 --> indica la variabile di scorrimento in questo caso 1, tale indice scorre sul termine appunto sopra-diagonale. Cioè


..e così via..
idof 3
idof 2
idof 1


Successivamente viene introdotto il ciclo do:

 do isd=1, min (Mband-1, idof-1)

Ora:

  • isd rappresenta l'indice di scorrimento
  • Min: su min dobbiamo porre molta attenzione in quanto è un oggetto molto costoso per il processore perchè va a richiamare l'if, tuttavia lo introduciamo per "abbreviazione didattica", cioè lasciamo che sia il processore a scrivere più righe di codice, risparmiando a noi il lavoro. Successivamente:
  b(idof-isd) = b(idof-isd) - A(idof-isd, 1+isd)*val
  Aband(idof-isd, 1+isd) = 0.0
  • idof-isd: indica la riga che sale
  • 1+isd: colonna che si sposta verso destra
  • val: δi (delta segnato i)

Osservazione:SOLO DOPO CHE LO LEGGO E LO PORTO A TERMINE NOTO POSSO ANNULLARLO

Nota: isd parte da 1, quindi non lavoriamo mai sul termine idof stesso, ma sempre da idof-1 in su.

Ora posso dedicarmi allo spostamento dei termini sotto-diagonali di A(*,idof) a termine noto e l'annullarsi degli stessi. Eseguo direttamente il ciclo do:

do isd=1, min(Mband-1, NSIZE-idof)
b(idof+isd) = b(idof+isd) - A(idof, 1+isd)*val
  Aband(idof, 1+isd) = 0.0
  • Mband-1: evita lo sfondamento della matrice A alla sua destra
  • Nsize: evita lo sfondamento dell'eventuale "triangolino inferiore di zeri della matrice A"
    |             |
    |             |
    |         _ |0|
    |      _ |0  0|
    |     |0  0  0|

Nota : La subroutine di cui sopra non coincide con la subroutine di vincolamento della Dispensa Progettazione Assistita, in quanto quest'ultima applica in sequenza tutti i vincoli prescritti, mentre la subroutine sopra illustrata è una specie di sotto-subroutine che vincola un solo grado di libertà alla volta. Pertanto, qualora si debbano vincolare n gradi di libertà, si procede con un ciclo FOR sugli elementi della lista e, per ogni passo del ciclo FOR, si effettua una chiamata alla subroutine.

PROBLEMI DI CONTATTO

I codici di analisi agli elementi finiti commerciali dispongono generalmente di specifici algoritmi per la gestione dei contatti tra parti diverse all'interno del modello, in modo da garantire che non ci sia compenetrazione tra di esse. La risoluzione dei modelli in tal caso non è infatti semplice, perchè se si ha presenza di possibili contatti si ha una discontinuità della risposta del sistema meccanico alle reazioni esterne. Le non linearità da contatto sono tra le peggiori che si possono avere in un modello agli elementi finiti, perché sono non linearità che possono dare discontinuità sulla tangente della risposta di un sistema. La risposta di un sistema meccanico ad una sollecitazione è sperabilmente continua. Al variare continuo delle sollecitazioni (forze applicate o spostamenti imposti), di solito si ha una variazione continua della risposta. (Di solito, poiché esistono controesempi, indotti dalla presenza dell'attrito). Per capire, vediamo un esempio semplice:

trave incastrata con possibilità di appoggio in mezzeria

Figura 1
Consideriamo l'esempio di figura 1, in cui si ha una trave di lunghezza l incastrata e caricata all'estremità libera da un carico P(t) lineare nel tempo. Supponiamo inoltre che sia presente un appoggio, a distanza "a" dall'incastro ed a distanza "s" dalla trave. Qual è la risposta del sistema? La trave viene caricata in modo variabile nel tempo con un carico P(t); nei primi istanti, dove P è ancora piccolo e perciò ancora non si è verificato il contatto con l'elemento di appoggio, l'andamento di δ (spostamento verticale del punto della trave sottostante al carico P) è lineare nel carico e quindi anche nel tempo; in pratica, il sistema si comporta come se quell'appoggio non ci fosse. All'aumentare del tempo, poi, la trave arriverà a contatto con l'appoggio; continuando a procedere nel tempo con questo tipo di deformazione si violerebbe il vincolo di non compenetrazione con l'appoggio. Da quel momento in poi la struttura si comporterà come appoggiata, quindi lo spostamento imposto del punto della trave che ora giace sull'appoggio sarà pari alla distanza s. Tale imposizione è ottenuta introducendo una reazione vincolare: infatti, l'appoggio restituirà una reazione vincolare verso l'alto. Questa reazione vincolare che si è venuta a creare sarà causa della discontinuità nella derivata di δ; infatti, come si vede dal grafico dello spostamento (δ) in funzione del tempo, il punto di discontinuità sarà presente nell'istante in cui la trave tocca l'appoggio, e successivamente δ avrà un andamento ancora lineare ma di pendenza minore, come mostrato in figura 2. Da quell'istante la natura del sistema viene a essere modificata; in particolare, la trave si comporterà con la rigidezza propria della trave incastrata e appoggiata: cioè avrà comportamento più rigido.
Figura 2


Per cui abbiamo in definitiva due risposte della trave in base a P(t):

  1. le configurazioni che essa assume prima del contatto con l'appoggio, che creano uno spostamento lineare
  2. le configurazioni che la trave assume dopo il contatto con l'appoggio: dove gli spostamenti sono ancora lineari, ma di pendenza diversa perchè dopo il contatto varierà anche la rigidezza (K) della trave, in particolare avrà maggior rigidezza in quanto la pendenza della curva diminuisce, per cui a parità di ΔP imposto (tra le due tipologie di configurazioni assunte) il Δδ sarà minore.
Figura 3

Se l'appoggio fosse posizionato all'estremità (figura 3), sotto il punto di applicazione del carico P, la trave si deformerebbe fino al contatto con lo stesso; a tal punto, la P(t) sarebbe totalmente equilibrata dalla reazione vincolare verticale. In questo caso, dopo la discontinuità, la deformazione δ risulterebbe costante (figura 4), cioè un incremento di carico ulteriore non causerebbe un'ulteriore deformazione elastica della struttura. Anche in questo caso, risulta evidente il fatto che la struttura evolve con derivata della funzione spostamento nel tempo discontinua.

Figura 4



Le reazioni di contatto possono essere espresse tramite uguaglianze o disuguaglianze; nel primo caso parliamo di vincoli olonomi o bilateri, nel secondo caso siamo in presenza di vincoli anolonomi o monolateri. Nel caso di problemi di contatto, siamo in presenza di vincoli del secondo tipo quindi analizziamo questi ultimi.



Vincolo monolatero

Il vincolo monolatero si ha, considerando l'esempio precedente, quando l'appoggio tocca la trave in uno degli estremi, ed è determinato da disequazioni poichè il gioco deve essere sempre positivo, mentre la reazione deve sempre essere repulsiva. Devono inoltre valere le seguenti condizioni di Signorini:

  • condizioni di Signorini:
 se g>0  => R=0
 se R>0  => g=0

relazioni che vengono soddisfatte al verificarsi del seguente prodotto scalare:

<g,R>=0

in cui g e R sono rispettivamente il vettore dei giochi e il vettore delle reazioni vincolari di tutti i nodi o punti della struttura.

Dove si ha gioco non nullo, la reazione è nulla. Dove c'è reazione non nulla, il gioco è nullo. Gioco e reazione possono essere solo ristretti in segno. In particolare per quanto riguarda i segni è buon costume orientare il gioco e le reazioni in modo che siano positivi quando si ha separazione per il gioco e repulsione per la reazione.

Contatto nodo su faccia

Come gestisce un codice agli Elementi Finiti, come il Marc, il problema dell'imposizione di stati di contatto monolatero tra corpi modellati agli elementi finiti?

Si supponga di lavorare su un dato valore di carico. In particolare, supponiamo di avere un valore di caricamento (in senso lato, in quanto può essere anche uno spostamento imposto) che varia nel tempo. Si supponga che una sollecitazione ad uno step j-simo sia stata approvata dai test di convergenza e considerata buona. Si vuole calcolare lo stato del sistema con un caricamento incrementato al valore j+1. Ci si muove da un valore noto, in cui si è risolta la struttura, verso un valore successivo, procedendo con un'iterazione di tipo Newton-Raphson, in quanto il sistema è genericamente non lineare. Si noti come se il caso non lo fosse stato, introducendo il contatto, lo diventerebbe.

Consideriamo la seguente figura, che rappresenta due corpi modellati agli Elementi Finiti, con profili lineari a tratti, a spezzata, tipici di oggetti modellati agli Elementi Finiti. La natura più semplice del contatto è un nodo di un corpo che entra nella faccia di un altro corpo. Questo si definisce contatto "NODO SU FACCIA".

L'algoritmo storicamente utilizzato dai programmi FEM per la gestione dei contatti è quello "NODO SU FACCIA" (Figura 2), per il quale il contatto avviene tra il nodo del corpo detto "touching" e la faccia del corpo detta "touched"; questo funziona correttamente soltanto se il corpo del quale consideriamo i nodi è meshato più finemente di quello di cui consideriamo le facce.

Quindi, il Marc usa la denominazione, di seguito riassunta:

--> touching (toccante) --> corpo di cui si considerano i nodi
--> touched (toccato) --> corpo di cui si considerano le facce
Figura5

Nella pratica però questo algoritmo è implementato in una sua versione più complessa, che va a migliorare la rappresentazione dei bordi dell'oggetto considerato rispetto a quelli ottenibili tramite meshatura. Per inquadrare il problema, consideriamo infatti il problema del contatto tra albero e mozzo supponendo non ci sia nè interferenza nè gioco, caso in cui non ci può essere trasmissione di momento. Se andiamo però a considerare una possibile meshatura per tali elementi (particolarmente grossolana per capire il concetto), si vede che il problema si trasforma praticamente nel contatto di una chiave a brugola e del relativo bullone, e quindi a meno di un gioco molto elevato sarà sempre presente una certa coppia, dovuta al contatto dei componenti meshati. Per evitare problemi come questo, l'algoritmo che gestisce i contatti nei software commerciali va ad estrapolare tramite una cubica passante per i diversi nodi del corpo touched quello che poi interpreterà come il bordo del corpo reale. Andrà quindi a gestire ugualmente il contatto a quanto faceva nel caso semplice precedente, ma in questo caso la compenetrazione negata sarà già quella tra nodo del corpo touching e bordo fittizio del corpo touched, come si vede da figura:

Figura6


Nella pratica, non è sempre detto che questa spline tra i nodi limite di un oggetto nel modello sia effettivamente conveniente per gestire i problemi di contatto. Se facciamo riferimento agli oggetti in cui sono presenti degli spigoli non dovuti a meshatura ma realmente esistenti, la spline va a modificare il bordo dell'oggetto creando avvallamenti o protuberanze. Di conseguenza potrei erroneamente dichiarare in compenetrazione un nodo che, nella realtà, è ancora fisicamente distaccato dalla superficie del pezzo. Oltre al problema di precoce contatto, potrei ottenere una pressione localizzata priva di senso, causata da una fittizia interferenza tra un profilo corretto in maniera sbagliata ed il lato opposto. E' quindi intuitivo come il possibile contatto tra un oggetto così modificato ed un altro non rispecchia la realtà, ed è per questo motivo che nel caso in cui il corpo presenti nel bordo una discontinuità reale (spigolo) e non dovuto a meshatura si va a considerare il bordo definito dalla mesh, quindi le interpolanti non sono presenti qualora si sia in prossimità di uno spigolo reale dell'oggetto, mentre sono utilizzate in tutti gli altri casi, soprattutto per contorni circolari o sferici. Il Marc, infatti, ha un controllo specifico delle discontinuità. In particolare, prevede che si dichiarino esplicitamente come discontinui tutti i punti attraverso cui non si vuole che passino le interpolanti; così facendo, si vieta al Marc di costruire interpolanti a cavallo dei nodi dichiarati discontinui. In un modello Elementi Finiti modellato a elementi lineari, tutti i nodi sono punti di discontinuità della derivata, però solo alcuni sono rappresentativi di una vera discontinuità.

Indifferentemente dalla configurazione considerata nel modello tra le due sopra descritte, supponiamo ora di caricare gradualmente la struttura con diversi passi di carico (LOADSTEPS). Tra i diversi loadsteps viene analizzato il modello tramite un algoritmo di Newton-Raphson in più iterazioni fino al raggiungimento di una certa tolleranza, per cui si assume soddisfacente la soluzione ottenuta e si procede al loadstep successivo. Supponiamo che all'i-esima iterazione sia stata raggiunta la convergenza del metodo di Newton-Raphson, ma non ancora la tolleranza richiesta tra le due iterate i-1 e i, quindi si procede con l'iterata i+1. Supponiamo però che all'i+1-esima iterazione il nodo compenetri la faccia o la spline, come mostrato nella configurazione di destra delle figure precedenti. L'iterazione ora non può convergere, poichè non è risolto il vincolo di contatto (e non compenetrazione) tra i due corpi in esame. Tale vincolo viene introdotto quando dichiariamo che vogliamo avere vincoli di contatto monolatero tra due corpi. In generale, tra due elementi di un modello Elementi Finiti, ci può essere compenetrazione, salvo indicazione contraria. Il fatto di dichiarare che due pezzi di mesh sono in contatto vuol dire attivare il controllo di compenetrazione tra i due corpi. Quando troviamo un'avvenuta compenetrazione tra un nodo e una faccia, sia essa in forma originale, sia in forma corretta (correzione di tipo puramente cosmetico), cambiamo la condizione relativa tra il nodo e la faccia. Prima della compenetrazione, lavoravamo in ipotesi che un dato nodo fosse libero di muoversi come voleva rispetto ad una data superficie. Tuttavia, tale ipotesi ha portato a compenetrazione. Trovata la compenetrazione, non possiamo più permettere al nodo di muoversi in maniera libera rispetto alla superficie, ma viceversa, creiamo un vincolo di tipo cinematico per cui il nodo deve muoversi come la superficie, nel grado di libertà normale alla stessa. Raggiunto il contatto subentra così il vincolo cinematico per il quale il nodo è vincolato a muoversi seguendo la superficie senza compenetrarla, come mostrato in figura.

Figura7


Questa condizione è possibile tramite la creazione di un SERVO LINK automatica da parte dell'algoritmo tra il nodo del corpo touching e la faccia dell'elemento touched, e sarà quindi riferito a 4 o 5 nodi di cui quello appartenente al corpo touching sarà dipendente dai 3 o 4 nodi (rispettivamente per elemento triangolare o quadrilatero) del corpo touched che definiscono la faccia (o la spline) considerata. Viene in questo caso definito un cambio di coordinate per il nodo "touching", del quale devo vincolare solo lo spostamento normale alla superficie "touched". Cioè il nodo che deve essere portato in contatto ad una superficie deve avere le incognite di spostamento così orientate: una normale alla superficie, le altre due tangenziali. La componente normale alla superficie deve perciò essere vincolata a muoversi in funzione della superficie stessa, quindi deve essere vincolata a muoversi in funzione dei nodi, più prossimi allo stesso, sulla superficie.

A questo punto, otteniamo un corretto iterato i+1, a cui succederanno gli iterati i+2, i+3, ..., fino a raggiungere un'eventuale convergenza. Ovviamente, tale vincolo è un vincolo di tipo bilatero, al quale sarà associata una reazione vincolare, che, se il nodo tende a compenetrare, sarà una reazione di tipo repulsivo.

Se la fase di carico fosse seguita da una successiva fase di scarico, durante la quale viene meno il contatto fra i corpi, l'individuazione del distacco avviene considerando le reazioni di contatto. Durante lo scarico vi sarà infatti un istante nel quale le reazioni di contatto diverranno inaccettabilmente trattive e non più repulsive come dovrebbe essere; dichiarando questa condizione inammissibile, quando la stessa si verificherà il Marc procederà automaticamente alla disattivazione del servolink fra nodo e superficie, ottenendo il distacco voluto fra i corpi.

Riassumendo, il carico prima cresce e poi cala. A un certo punto, durante la crescita del carico, attivo un servolink, il quale viene mantenuto per tutta la porzione in cui un dato nodo diversamente tenderebbe a compenetrare. Arriverà, prima o poi, un istante in cui il servolink giungerà a forza nulla (caso di sfioramento, realmente improbabile, in quanto condizione singolare e di fatto istantanea). Seguirà sicuramente uno step di scarico, in cui la tying force (forza legante) associata al link risulta non più repulsiva, ma tende a trattenere il nodo contro la superficie. A quel punto, l'iterato che genera una forza repulsiva deve essere dichiarato inammissibile. Il servolink va eliminato e si continua l'analisi, non più supponendo che quel nodo sia vincolato bilateralmente alla superficie con il servolink, ma considerando quel nodo indipendente dalla superficie.

Quindi, il Marc gestisce il contatto attivando e disattivando dinamicamente dei servolinks sul grado di libertà di spostamento normale alla superficie del corpo, con una lievissima variante, per cui i test di compenetrazione sono fatti durante l'iterazione Newton-Raphson, quindi tra l'iterato i e l'iterato i+1, mentre i test di forza trattiva vengono fatti esclusivamente a convergenza avvenuta.