Indice
Analisi della subroutine per l’assemblaggio di una matrice di rigidezza relativa ad un elemento finito isoparametrico 4 nodi piano in una matrice di rigidezza globale bandata
Con AK(8,8)
viene passata la matrice di rigidezza dell’elemento: matrice quadrata con 8 righe (e 8 colonne) pari ai gdl complessivi dell’elemento (u e v per ognuno dei 4 nodi);
con STRUTK(2*nodes,mband)
viene passata la matrice di rigidezza globale bandata
Dove 2*nodes indica il numero di gdl complessivi ed mband indica la larghezza di banda;
IPOINT(8)
è un vettore (riga) fondamentale per la compilazione di STRUTK (piena), in quando riporta in sé la correlazione esistente tra gli indici dei gdl nella matrice di rigidezza di elemento e gli indici degli stessi nella matrice di rigidezza globale;
JBSTRK
= indice di colonna della matrice bandata.
ES: iso4 di nodi [3,7,11,5]
(1,2,3,4,5,6,7,8) | numerazione implicita gdl rispetto ad elemento iso4 | |
Ipoint | =(5,6,13,14,21,22,9,10) | numerazione gdl rispetto a struttura globale |
---|---|---|
(u3x,v3y,u7x,v7y,u11x,v11y,u5x,v5y) | identificazione implicita grado di libertà |
La parte esecutiva della subroutine è composta da due cicli do
nidificati con un indice “i” che scorre da 1 a 8 per le righe di AK e un indice “j” che scorre sempre da 1 a 8 per le colonne di AK, ciò significa che il blocco di operazioni centrale viene ripetuto per ogni elemento della matrice di rigidezza dell’elemento.
Per ogni elemento di AK (tramite il vettore ipoint) si riesce a risalire alla posizione esatta che esso occupa nella matrice di rigidezza globale piena ed è quindi possibile assemblarlo poi nella matrice di rigidezza globale bandata. Questo “doppio” passaggio si deve fare poiché non è possibile sapere a priori se un determinato elemento di AK sta sopra o sotto la diagonale della matrice di rigidezza globale piena, per cui, una volta determinati gli indici tramite le operazioni precedenti, occorre fare un controllo sulla sua posizione effettiva e se tale controllo è verificato allora è possibile procedere all’assemblaggio dell’elemento di AK nella matrice bandata globale.
Oss: un controllo alternativo potrebbe essere
if (JBSTRK.GE.1) then
Analisi della subroutine per il processo di vincolamento in una matrice di rigidezza bandata
idof
= gdl scelto per il vincolamento (idof-esima riga della matrice globale)
val
= valore imposto come vincolo
nsize
= 2*nodes = numero gdl totali
aband
= strutk (bandata)
b
= forces (termine noto)
Il processo da seguire per il vincolamento di un gdl è sempre quello visto nelle lezioni precendenti (Lezione 6), ovviamente però la subroutine deve essere riadattata per poter lavorare con una matrice non più piena ma bandata.
Operativamente parlando, nella matrice aband il primo elemento della riga idof-esima corrisponde a quello diagonale (nella matrice piena) al quale va sostituito “1”; gli elementi sopradiagonali (appartenenti alla banda nella matrice di rigidezza globale piena) devono essere presi in diagonale a 45°; quelli sottodiagonali (non stoccati direttamente), siccome è possibile sfruttare la proprietà di simmetria della matrice piena, sono tutti quelli (tranne il primo) posti alla riga idof-esima della matrice bandata; in questo modo si riesce a “ricostruire” tutta la idof-esima colonna della matrice piena (utile concettualmente per seguire passo-passo in processo di vincolamento già visto).
Sono implementati due cicli do
separati, uno per trattare gli elementi sopradiagonali e l’altro per quelli sottodiagonali in cui scorre un indice “j”; tali indici scorrono da 1 al minimo tra due grandezze che indicano sostanzialmente i limiti da non superare per rimanere all’interno della banda e all’interno della matrice piena, questo per non rischiare di andare a cercare termini che non sono dentro la banda, quindi non sono da considerare assolutamente o che addirittura non esistono!
Per seguire il processo di vincolamento, all’interno di ogni ciclo, scorrendo l’indice di riga, (muovendosi quindi lungo l’idof-esima colonna della matrice piena) viene selezionato un elemento, moltiplicato per val, cambiato di segno e aggiunto all’elemento ternime noto corrispondente (con stesso indice di riga), una volta fatta questa operazione il termine considerato in partenza nella matrice di banda viene azzerato.
Oss: le due operazioni non possono essere invertite, altrimenti se prima si azzera poi si moltiplica di fatto non si aggiunge nessun contributo!
Oss: la subroutine prende al suo interno la matrice STRUTK, la modifica e la restituisce al main
Successivamente nel programma principale si può proseguire nel calcolo della soluzione attraverso una ulteriore subroutine “solutore” sempre tenendo presente che deve essere adattata per lavorare con una matrice bandata; la particolarità di questa subroutine è che la soluzione viene sovrascritta direttamente sul vettore del termine noto, evitando così di dover allocare un altro vettore solo per le soluzioni; in oltre usa metodi di calcolo molto efficienti e veloci.
Approfondimento su singolarità tensionali
I casi qui riportati e discussi sono stati riassunti e presentati nella seguente pubblicazione:
sinclair_i.pdf
(le figure cui si fa riferimento sono riportate a pagina 4 del documento)
Caso A) il carico 2F risulta applicato come “concentrato lungo una linea” ⇒ ordine di singolarità= $1/r$
dove r è la distanza dal punto di singolarità, ovviamente tale relazione vale solo nell’intorno del punto;
nell’intorno del crack (spigolo vivo rientrante con angolo di apertura uguale a zero) in corrispondenza di P2 si ha ordine di singolarità= $1/\sqrt{r}$
Caso B) gomma – asfalto la buca offre uno spigolo alla gomma quindi si rientra nel caso di un oggetto elastico indentato da uno spigolo rigido ordine di singolarità= $(1/\sqrt{r})ln{(r)}$
In corrispondenza di P6 c’è un'altra singolarità poiché il pneumatico (schiacciato dal peso) tende ad allargarsi, per sui si instaura uno scorrimento tra la superficie della gomma e la superficie rigida del suolo, in conseguenza di ciò si creano delle tensioni laterali $\tau$, del tipo: $P_s*f_a$ (pressione schiacc.*coeff.attrito);
dalla teoria si sa che le $\tau$ vanno sempre in coppia (per equilibrio alle rotazioni del cubetto di materia), ma essendo il punto analizzato in corrispondenza del bordo del pneumatico la $\tau$ lato “aria” non c’è perché proprio non esiste, per cui si arriva ad un “nonsenso della teoria dell’elasticità” quindi si dice che c’è uno stato tensionale singolare, in cui la $\tau$ deve andare a zero immediatamente per cui assume un valore singloare.
Caso C) si hanno tre corpi elastici: cilindro, anello, pistone di cui uno (anello) con spigoli vivi che viene spinto ad indentare gli altri due con conseguente effetto lama , in realtà è un effetto un pò più debole poiché lo spigolo di apertura è di 90°.
Caso D) cava per chiavetta in albero, non c’è singolarità a spigoli convessi mentre c’è singolarità in corrispondenza degli spigoli rientranti ordine di singolarità= $1/\sqrt[3]{r}$;
si osserva che la condizione di singolarità dipende dalla geometria locale e non da quella complessiva e neanche dal materiale, vale a dire che lo spigolo rientrante a 90°, proprio per come è fatto, comporta una condizione di singolarità sempre, anche al di fuori di questo caso presentato.
Caso E) contatto multi-materiale: incollaggio di testa: acciaio_ epossidica_ acciaio ; l’ordine di singolarità è funzione di E e $\nu$
Caso F) problema simile al caso B in quanto si ha uno scorrimento relativo questa volta tra un oggetto rigido sopra e uno più cedevole sotto;
dal punto di spigolo verso destra c’è pressione e moto, quindi anche tensioni tangenziali $\tau$ (sempre del tipo: $P_s*f_a$), mentre dal punto di spigolo verso sinistra non c’è niente, il materiale sottostante è completamente scarico, per cui non ci sono nemmeno le $\tau$; questa situazione rappresenta un altro “nonsenso”.
Caso G) effetto lama quindi singolarità dovuta al tendere della pressione all’infinito.
caso H) caricamento concentrato su un singolo nodo quindi pressione “matematicamente” infinita.
Esempio di cattiva modellazione FEM
Modellazione di cuscinetti a sfera di auto da competizione: caso 3D con iso8
situazione: albero su due appoggi caricato da una forza ad una certa distanza da uno dei due cuscinetti
scopo: ottenere freccia albero date certe condizioni al contorno tra cui carico applicato e caratteristiche dei cuscinetti (precarico)
L’albero è modellato tramite una linea che è intesa come asse di una trave che possiede tutte le caratteristiche e le proprietà geometriche dell’albero stesso;
modellazione dei gusci (piste) esterno ed interno come elementi deformabili, con guscio esterno incastrato in supporto (si potrebbe fare anche con lieve interferenza) e guscio interno modellato in modo che si muova solidale all’albero su cui è montato; questo vincolo si realizza sostanzialmente collegando i nodi della parete interna del guscio interno con l’asse rappresentate l’albero tramite elementi filiformi indeformabili in modo da avere moti di corpo rigido di questi elementi, di conseguenza il moto dell’albero si trasferisce interamente al guscio interno.
Secondo una prassi ormai ben definita e non messa in discussione le sferette all’interno dei cuscinetti sono modellate mediante delle travette che partono da un nodo sulla parete interna della pista esterna e arrivano ad un nodo della parete esterna della pista interna, tali travette hanno sezione data in funzione delle sferette e materiale dato;
le travette lavorano in realtà come puntoni, poiché essendo collegati ad elementi che non portano le rotazioni è come se avessero una cerniera sferica in testa.
Oss
differenza tra trave e puntone (cavo,tirante): la trave reagisce a sforzo normale, taglio e flessione (anche a torsione), mentre il puntone reagisce solo a sforzo normale (compressivo o trattivo)
Ciò che non va in questa modellazione è il fatto che le travette, essendo collegate in testa ad un solo nodo, esercitano sul corpo una forza che è concentrata, per cui si hanno inevitabilmente delle condizioni di singolarità tensionale.
Anche se si usa l’accorgimento di “spalmare” la forza concentrata su di un’area fittizia nell’intorno del nodo (considerando parte degli elementi che concorrono al nodo stesso), si nota che più fitta è la mesh, più piccoli sono gli elementi che concorrono a quel nodo quindi più piccola è l’area su cui si “spalma” la forza applicata, pertanto nel momento in cui si infittisce la mesh si ha che contestualmente si amplifica la singolarità tensionale, quindi si ha una deformazione compatibile con la caratteristica di forza concentrata (tipo: lenzuolo teso pinzato) in cui le tensioni teoricamente tenderebbero ad infinito (tangente verticale).
Si osserva che il solutore, per ragioni di calcolo, opera una troncatura della soluzione ad un valore finito, per cui restituisce una cedevolezza del cuscinetto finita che però non tiene conto assolutamente delle deformate in soluzione reale che, nei punti di applicazione di forza concentrata, sarebbero infinite;
per cui in definitiva si ha che la soluzione discreta porta a dei risultati della cedevolezza che, all’aumentare dell’infittimento della mesh, al posto che convergere andranno a divergere fino ad essere considerati addirittura non fisici!!!
Pertanto tale modellazione è molto poco robusta
Solitamente questo tipo di modellazioni sono utilizzate con alcune regole di taratura del modello da seguire, fissata la taglia e il tipo dell’elemento di mesh; ad esempio le travette posssono essere collegate in testa a più di un singolo nodo fissando così la grandezza dell’area su cui si pensa applicata la forza, infittendo la mesh ci saranno più elementi all’interno dell’area che però rimarrà invariata
Attenzione: questo tipo di modellazione potrebbe portare a grossi errori di valutazione!