a cura di Ludovico Posa e Alessandro Rando

Partiamo ricordando che le funzioni peso per l'interpolazione di funzioni nodali sull'elemento sono di natura BILINEARE, lungo i lati però la funzione bilineare è semplicemente LINEARE.

L'interpolazione che utilizziamo per modulare una quantità da n1 a n2 (come nella figura sopra) prevede una modulazione lineare. Per ricavare la funzione peso nel punto considerato si sfrutterà quindi la semplice regola della leva.

Supponiamo ora di avere due elementi:

I nodi a cavallo (g132 e g3) sono disegnati separati per chiarezza ma, sono in realtà lo stesso nodo. Consideriamo un ipotetico spostamento e un ipotetica rotazione ,del nodo indicato con g3, che rappresentiamo in forma di vettore ( rosso in figura). Sono questi due i vettori che definiscono la rototraslazione di quel nodo. Ovviamente anche il nodo g3 visto dall'altro elemento avrà gli stessi spostamenti e le stesse rotazioni (essendo questi due nodi (disegnati separati per chiarezza) in realtà lo stesso nodo). Essendoci un unico grado di libertà di spostamento in x per tale nodo, non potrà avere due diversi spostamenti in x. Facciamo ora ragionamenti analoghi per il nodo indicato con g132 (spostamenti e rotazioni indicate in blu). Quindi ALMENO in corrispondenza di quei due nodi gli spostamenti sono indubbiamente continui. La funzione spostamento non vede salti passando da quei due nodi perchè all'interno dell'elemento l'interpolazione è una forma continua e non posso avere discontinuità e ai nodi a cavallo ho continuità (gli spostamenti visti dai due elementi sono uguali). La questione adesso è valutare se un punto preso sul lato compreso tra i nodi g3 e g132 (che ricordiamo essere lo stesso lato e sono disegnati separati solo per maggior chiarezza) si sposti in maniera diversa tra i due elementi.

Consideriamo quindi un generico punto preso su quel lato (che nella figura in alto è contrassegnato in verde). Come definiamo gli spostamenti di quel punto? Tali spostamenti si potranno definire a partire dagli spostamenti nodali, per interpolazione. L'interpolazione segue una legge lineare e gli spostamenti nodali sono uguali tra i due elementi (come stabilito precedentemente) quindi sarà semplice definire gli spostamenti di quel punto! Inoltre, l'interpolazione ha un'altra peculiarità, ovvero, il nodo indicato con g73 (per esempio) ha una funzione peso che ha influenza nulla sui due lati che non toccano quel nodo (quindi sui lati g132-g3 e g132-g8). Quindi in definitiva, non importa come si sposta il nodo g73 se vogliamo valutare lo spostamento del nostro punto,perchè il moto di quel nodo viene pesato con peso identicamente nullo. Se tale influenza non fosse nulla, sarebbe una potenziale fonte di discontinuità, perchè ovviamente, nulla lega il moto del punto g73 con quelli dell'altro elemento. Quindi, dato che nei nodi ho valori comuni e dato che nel trovare gli spostamenti del punto considerato utilizziamo gli stessi schemi di interpolazione lineare possiamo dire che il punto avrà lo stesso modo tra i due elementi. Osserviamo adesso la natura lineare dell'interpolazione DELL'ELEMENTO DI DESTRA. Poniamo il punto a una certa quota che indichiamo con λ che indica la distanza tra il nodo g3 e il punto contrassegnato (come si può vedere dallo “zoom” nella figura qui in basso del lato in comune).

Dalla figura si può notare che quando λ è uguale a 0,5 praticamente siamo a metà del lato. Ovviamente il tratto che va dal punto contrassegnato al nodo g132 sarà pari a (1-λ). Il vettore spostamento del punto contrassegnato si sposterà di una quantità definita dalla somma degli spostamenti dei due nodi (ug3 e ug132) moltiplicati per i loro corrispettivi pesi (rispettivamente (1-λ) e λ). Lo spostamento del punto considerato è definito sotto forma vettoriale nella figura qui in basso:

Questo è lo spostamento del punto interno al lato SECONDO L'ELEMENTO DI DESTRA ma, attraverso considerazioni analoghe si può ricavare che tale spostamento sarà uguale ANCHE PER L'ELEMENTO DI SINISTRA, perchè ovviamente si interpolano linearmente i medesimi valori. Stesso discorso vale per le rotazioni. Quindi in definitiva per tali motivi c'è continuità passando da un elemento all'altro. Ricapitolando la continuità è garantita da: -Le leggi di interpolazioni sono da entrambe le parti lineari -Non c'è influenza dei nodi non condivisi tra elemento e elemento -I valori estremali sono comuni.

Cambiamo argomento. Voglio descrivere ora, lo schema di interpolazione delle componenti delle deformazioni fuoripiano (dette anche Traverse) che utilizza l'elemento 75 di MARC per evitare di incorrere nello SHEAR LOCKING. Consideriamo la seguente struttura:

Allora, di base le deformazioni vengono campionate a centro lato. Perchè a centrolato? Perchè è l'unico punto nel quale a fronte di una pura curvatura (la pura curvatura si nota dalla figura in basso, guardando l'elemento dal fianco lungo z del lato 12, infatti ho lo stesso angolo -alfa e + alfa, c'è quindi una semplice rotazione),non vedo nessuna deformazione tagliante spuria. Al contrario, ai punti estremali vedo delle deformazioni taglianti spurie. L'oggetto rappresentato in verde quindi, nasce rettangolare e rimane rettangolare anche dopo la curvatura(forse deforma un pò a trapezio ma per semplcità lo assumiamo rettangolare). Questa curvatura uniforme che si sviluppa parallelamente al lato lo va a flettere.

Stessi ragionamenti li possiamo fare sul lato compreso tra i nodi 1 e 4.

Quindi in sostanza il centrolato è l'unico punto dove non si vede deformazione spuria ed è proprio qui che campiono deformazione tagliante fuoripiano utilizzando la seguente formula già vista precedentemente:

Dove “d” è il vettore dei gradi di libertà dell'elemento e la matrice B_γ(ξ,η), è funzione di ξ e η i quali in questo caso vengono sostituiti dalle coordinate del centrolato.

Quindi quando campiono il centrolato ottengo due componenti di γ: γ_zx γ_zy

[P.S. Ricordiamo che stiamo parlando ovviamente di deformazioni fuoripiano, perchè il fenomeno dello SHEAR LOCKING associato a deformazioni entropiano non è compensato nell'elemento 75, ovvero nei modi elementari 1 e 2 è compensato ma non per i modi 8 e 10 (che sono presenti nella figura in basso)]

Comunque,dovrei quindi avere due componenti x e y di deformazione ma, preferisco proiettare le due componenti lungo il lato tirando fuori la componenti taglianti parallele e ortogonali al lato (in pratica scompongo una componente γ allineata al lato 12 (che è il lato preso in esame) e una ortogonale al lato 12). Quella ortogonale non la considero perchè non ha nessuna proprietà particolare che ci interessi e tengo solo quella allineata al lato. Quindi in poche parole considero solo l'angolo deformante a taglio che vedo sul fianco laterale e non l'angolo che si vede nel piano ortogonale.

Quindi ricapitolando, campionando a centrolato non si risente della deformazione tagliante spuria dovuta alla pura curvatura. Ma, se ci fosse una deformazione tagliante legittima si riuscirebbe a vederla campionando a centrolato? Per scoprirlo aggiungiamo una deformazione tagliante legittima. Manteniamo gli stessi angoli dei segmenti ai nodi (per esempio nella figura in basso si vede che avremo sempre gli angoli alfa visti precedentemente) però è stato aggiunto uno spostamento fuori piano differenziale:

Tale spostamento fuori piano differenziale fa in modo che i due nodi si spostino uno verso l'alto e uno verso il basso inducendo sul fianco dell'elemento (insieme alla curvatura pura dovuta alle rotazioni alfa che erano già presenti) una deformazione a mazzo di carte tagliante che osservo sempre sul piano Z - lato12. L'angolo sul rettangolino verde che si viene a formare (come si vede nella figura in basso,sempre guardando dal fianco del lato 12) lo chiamo γ_z12.

Tale angolo è opposto alla convenzione per i tagli positivi, cioè questo angolo è in realtà -γ_z12.

Quindi nel caso di presenza di una deformazione tagliante legittima,nonostante abbia campionato al centrolato, si rileva un deformazione tagliante non nulla, per cui in questo caso sono in grado ci campionare tale deformazione tagliante.

Ora, tornando al discorso di prima, campiono la proiezione di γ_z12 (componente parallela al lato) sul lato corrispondente e con ragionamenti analoghi campiono la γ_z41 sul lato 41. γ_z12 e γ_z41 posso trattarle come vettori nel piano normale a z.

In particolare posso fare una costruzione geometrica,ovvero prendo l'angolo che definisce la pendenza puramente tagliante (precedentemente definito come -γ_z12) e ne faccio la tangente che sarà uguale al γ vettore, orientato su questo piano, come il lato in esame. Stesso procedimento faccio per γ_z41 sul lato 41. Quindi ottengo così le due rappresentazioni vettoriali che appaiono nella figura in alto.

Dopo aver campionato le tensioni taglianti a centro lato, considero ciò che accade come rappresentativo di ciò che accade nell'intero lato in esame. Al nodo in comune tra i due lati (quello cerchiato in figura) sono associate ambo le deformazioni taglianti. Quindi la deformazione tagliante complessiva in quel nodo dovrà avere la seguente peculiarità: se proiettata sul lato 12 sarà uguale al vettore γ_z12 e se proiettata sul lato 41 sarà uguale al vettore γ_z41. Da questo vincolo appena citato, riesco univocamente a ricavare il vettore di deformazione fuoripiano associato a quel nodo (che sarà ricavato tramite semplice somma vettoriale delle due componenti parallele ai due lati).

Dopodichè di tale vettore posso estrarre le due componenti rispetto agli assi locali fisici x e y, ottenendo così il valore nodale di -γ_yz e -γ_zx al nodo 1:

Tale valore è diverso da quello che si otterrebbe andando ad applicare la formula scritta precedentemente (quella con la matrice B_γ(ξ,η)). A questo punto si ripete la procedura per ognuno dei 4 nodi ottenendo i valori di γ_yz e γ_zx ad ogni nodo. E una volta che si hanno i valori di ogni nodo tramite interpolazione posso ottenere i valori di γ_yz e γ_zx di ogni punto dell'elemento e in particolare ai 4 punti di integrazione che uso per le altre componenti di deformazione e tensione (le epsilon x e y e le sigma x e y). Il vantaggio di questo metodo è che il campo di deformazioni fuoripiano, che ottengo in tal modo, è filtrato rispetto alle deformazioni spurie che danno SHEAR LOCKING (perchè parto dai centro lato che come detto prima non risentono delle deformazioni spurie dovute alla pura curvatura). Una volta che si hanno i valori di ogni nodo tramite interpolazione posso ottenere i valori di γ_yz e γ_zx ai punti di Gauss, da quelle γ posso estrarre le componenti di tensione τ_yz e τ_zx, per poter quindi comporre, insieme alle altre tensioni (calcolate separatamente secondo quadratura gaussiana), il tensore delle tensioni completo ai punti di Gauss e poter infine calcolare la tensione secondo Von Mises e vedere se questa sia al di sotto della tensione di snervamento per poter verificare l ipotesi elastica (o al di sopra ipotesi plastica). Il vantaggio principale di utilizzare gli stessi punti di Gauss per l'integrazione di tutti i contributi energetici è proprio quella di sancire qual è il legame costitutivo (matrice D) da utilizzare tra quello elastico o quello plastico.

In definitiva con questo sistema di “campionamento delle deformazioni taglianti fuoripiano a centrolato” descritto dalla guida di Marc, riesco a togliere il problema dello shearlocking utilizzando sempre 2 punti di campionamento per asse secondo le regole di quadratura gaussiana e senza generare modi incompatibili alterando le funzioni di forma. Comunque, il Campionamento che utilizziamo filtra il termine spurio che genera shearlocking, ma ovviamente non sparisce, semplicemente si fa in modo di non osservarlo.

La formula che è stata scritta prima con la matrice B_γ(ξ,η) vale ancora. Cioè è ancora vero che per ogni punto dell'elemento γ_yz e γ_zx sono uguali alla matrice B_γ(ξ,η) moltiplicato per il vettore dei gradi di libertà nodali d però si deve aggiungere che la matrice B_γ è stata formulata in modo da contenere i contributi di tutti i passaggi appena visti. Quindi è una matrice leggermente diversa e per identificarla si può aggiungere un asterisco (indicante proprio l'assenza del contributo di taglio spurio che da problemi di shearlockinig) :

B_γ *≠ B_γ

CHIUDIAMO QUI CON L'ANALISI A LIVELLO DI SINGOLO ELEMENTO

Andiamo ad analizzare ora:

La struttura di riferimento è questo oggetto nella figura in alto a forma di “tetto” che è composto da 8 nodi e 4 elementi. Da cosa è descritta questa struttura? Ogni nodo è posizionato nello spazio mediante le sue coordinate in un sistema di riferimento globale. Questi elementi sono degli elementi quadrato di lato 4 (4 unità di lunghezza) con spessore 1 (quindi elementi abbastanza tozzi per non vedere sparire la deformata tagliante rispetto a quella flessionale) e con un angolo rispetto al piano orizzontale pari a 30° (quindi c'è un semplice disallineamento tra gli elementi, in modo da avere una differenziazione tra sistemi di riferimento locale e globale, cioè gli assi z degli elementi non sono allineati tra loro (tutti i termini trigonometrici saranno dei seni o coseni di 30°)). C'è un sistema di riferimento globale definito dai versori di x,y,z : i*_g,j*_g,k*_g. Essendo 3 versori non è importante dove sia l'origine che si può posizionare arbitrariamente. Su ogni elemento vi saranno degli assi locali rappresentati tramite i versori i,j,k di cui i e j sono orientati proprio come i lati di ciascun elemento (mentre ovviamente k sarà ad essi ortogonale). Quindi ogni elemento ha un sistema di riferimento globale diverso. Una volta definite le coordinate nodali [noi non le abbiamo definite ma è abbastanza semplice, basta definire un origine (per esempio al nodo 5 centrale) e tramite esso ricaviamo tutte le altre coordinate nodali dato che sono tutti elementi quadrati e il disallineamento tra loro è definito da angoli di 30°] occorre definire la matrice che lega i nodi di elemento con i nodi globali di struttura. Tale matrice si chiama matrice di CONNETTIVITA' (che fa parte degli input del codice Marc per descrivere la mesh) che in realtà, più che una matrice, è una semplice tabella. In pratica la struttura è descritta da una parte dalle coordinate nodali, mentre dall'altra dalla seguente tabella :

Questa tabella vista in figura quindi, non è altro che la tabella di corrispondenza tra numerazione locale e globale. Ad esempio se si vuol sapere qual è il terzo nodo locale dell'elemento 2 basta guardare la tabella e si vedrà che corrisponde a : e₂n₃ ≡ g6. Praticamente dalla tabella si può facilmente vedere quali sono i nodi di ogni elemento e come sono ordinati.

Nel tentar di far interagire e elementi , per esempio, l'elemento 1 con l'elemento 2 occorre definire un comune sistema di gradi di libertà. Ogni nodo a livello globale può avere un suo sistema di riferimento, ad esempio se consideriamo il nodo 5, posso definire per quel nodo una direzione i_g5 che definisce il primo asse, una direzione j_g5 che definisce il secondo asse e una direzione k_g5 che definisce il terzo. Quindi, quando definisco 3 componenti di spostamento, il vettore spostamento finale avrà le 3 componenti secondo gli assi del sistema di riferimento proprio del nodo. Questo sistema di riferimento è una proprietà del nodo, per cui è un oggetto che è visto in maniera comune da tutti gli elementi che toccano quel nodo (funge da “stele di rosetta” per comunicare tra elemento e elemento). Ora, per semplicità, tutti i sistemi di riferimento nodali sono allineati col sistema di riferimento globale citato prima [i*_g,j*_g,k*_g]. L'asterisco su quei versori del globale sta ad indicare che si riferisce a tutti i sistemi di riferimento nodali. (Sostanzialmente per il sistema di riferimento nodale ho una totale arbitrarietà che nel caso in esame ci porta ad agire così per semplificare le cose).

E' possibile anche associare ad ogni nodo un sistema di riferimento diverso, sembra che potrebbe complicare le cose, ma in alcuni casi è molto utile. Ad esempio, se si volesse lavorare in una struttura di riferimento a geometria anulare, per esempio un mozzo che viene calettato su un albero:

Può essere comodo descrivere il moto dei nodi sulla superficie interna non secondo un sistema cartesiano, bensì secondo un sistema polare che ha origine sull'asse dell'albero su cui il mozzo viene calettato. Per cui il nodo g123 (che si vede dalla figura in alto) avrà un primo asse locale (i*_g123) radiale, un secondo asse (j*_g123) tangenziale e un terzo asse (k*_g123)in direzione assiale. In tal modo vengono definite localmente le direzioni radiali, tangenziali e assiali del nodo g123. Il nodo g332 avrà altre 3 direzioni locali che saranno differenti da quelle del nodo g132 (apparte la k ovviamente). Quindi in questo modo, essendoci un sistema di riferimento nodale diverso per ogni nodo e definito secondo queste direzioni, è possibile (per esempio imporre) spostamenti radiali nulli a tutti nodi (invece che imporre spostamenti lungo x o y che in questo caso sarebbero inutili) o tangenziali nulli. Quindi solo in casi di necessità come questo, può essere comodo far si che il sistema di riferimento sia diverso per ogni nodo.

Tornando alla nostra struttura di riferimento, nella quale tutti i sistemi di riferimento nodali sono allineati col sistema di riferimento globale. Adesso voglio prendere tutti i gradi di libertà locali (spostamenti e rotazioni) e trovare per ognuno di essi una controparte globale. Concentriamoci sull'elemento 1 (e poi il procedimento che stiamo per descrivere si potrà ripetere per tutti gli elementi).

L'elemento 1 (così come anche gli altri elementi) ha 24 gradi di libertà associati (20 se non consideriamo i moti di drilling), di questi 24 G.D.L. ne consideriamo solo 2: lo spostamento in direzione locale z del nodo locale 2 (w_e1n2) e la rotazione attorno all'asse locale x del nodo locale 1 (θ_e1n1).

w_e1n2 è uno scalare per cui lo trasformo in una quantità fisicamente vettoriale associata (ovvero la freccia nera in figura in basso) che è semplicemente definita dal prodotto di questo scalare (w_e1n2) per il versore associato dell'elemento, k_e1.

(Ovviamente il moto di quel nodo locale sarà definito dalla somma dei contributi degli altri spostamenti v_e1n2·j_e1 e u_e1n2·i_e1, ma noi stiamo considerando lo spostamento del nodo solo di un asse locale dell'elemento (quello in direzione k))

Questo vettore che abbiamo ottenuto è effettivamente un vettore fisico, per cui posso scomporlo nelle 3 componenti x,y,z secondo il sistema di riferimento nodale (in questo caso del nodo g2). Poichè dalla tabella di connettività notiamo che e₁n₂ ≡ g2 , considero le direzioni i_g2,j_g2 e k_g2 (che fanno da ponte tra elemento e elemento) per cui, rappresento il vettore w_e1n2·k_e1 dall'equivalente nel sistema globale associato al nodo 2. Per cui rappresento il vettore w_e1n2·k_e1 attraverso la seguente equivalenza:

Quindi questo stesso vettore nello spazio posso rappresentarla attraverso queste due notazioni equivalenti, da cui posso facilmente ricavare la trasformazione di coordinate. In particolare ricaveremo che i termini di quella equivalenza corrispondono alle componenti del vettore w_e1n2·k_e1 (frecce rosse nella figura in alto), orientate secondo gli assi del sistema del riferimento nodale.

[P.S. ci sarebbe anche una terza freccia rossa dovuta al contributo w_g2·k_g2 ma è stato evitato per maggiore chiarezza nel disegno]

Ragionamenti analoghi li possiamo ripetere per la rotazione attorno all'asse locale x del nodo locale 1, θ_e1n1 (freccia nera nella figura in basso).

Quindi per questo G.d.L., l'equivalenza tra notazione locale di elemento e notazione globale del nodo è esprimibile dalla formula:

Quindi abbiamo:

Le due equivalenze mostrare sono solo 2 di esempio, ma ricordiamo che questo procedimento deve essere ripetuto per ogni GdL di ogni elemento.

Adesso, come opero la trasformazione di coordinate? Allora considero la prima equazione (quella con w) in cui la sua LHS (Left Hand Side) sarà ovviamente uguale alla RHS (Right Hand Side): <LHS>= <RHS>. Proietto, ora, ambo i membri su una direzione specifica che è k_e1 stessa,per cui moltiplico scalarmente ambo i membri per k_e1: <LHS · k_e1 > = <RHS · k_e1 > per cui ottengo:

In pratica ottengo che il grado di libertà locale w_e1n2 è definito come la combinazione lineare dei 3 gradi di libertà globali spostamento del nodo g2 per dei coefficenti che sono i coseni degli angoli compresi tra il versore locale k_e1 e i 3 versori degli assi globali del nodo g2. Ottengo così un banale cambio di coordinate.

Quindi abbiamo visto che è possibile scrivere una generica relazione tra gradi di libertà locali e gradi di libertà globali. Cerchiamo adesso di scrivere il tutto in forma compatta. Costruisco un vettore colonna con tutti i gradi di libertà globali della struttura, e lo chiamo vettore d_g:

Tale vettore colonna (che ricordiamo non essere un vettore fisico, bensì una matrice colonna) contiene spostamenti e rotazioni di ogni nodo, scomposti rispetto ai sistemi di riferimento nodali. Tale vettore è composto da 54 termini (9 nodi per 6 GdL a nodo). I termini sono stati raggruppati secondo un ordine preferenziale che è quello di mettere vicini i GdL associati a un singolo nodo (raggruppamento Base-Nodo).

In maniera analoga costruisco dei vettori colonna con tutti i gradi di libertà locali specifici di ogni elemento. Per esempio, considero l'elemento 1 e il vettore ad esso associato lo chiamo vettore d_e1:

Tale vettore colonna contiene spostamenti e rotazioni di ogni nodo dell'elemento, scomposti rispetto al sistema di riferimento locale di elemento. Tale vettore è composto da 24 termini (4 nodi per 6 GdL a nodo). Qui però il raggruppamento è stato fatto per componente, perchè in seguito semplificherà i passaggi algebrici (anche se solitamente il raggruppamento si fa mettendo vicini i GdL associati a un singolo nodo come è stato fatto per il vettore d_g, però per semplicità noi raggruppiamo come in figura).

Adesso devo trovare una sorta di equivalenza tra questi due vettori appena descritti (e poi si dovrà ottenere ovviamente anche per tutti gli altri elementi, intanto partiamo dall'elemento 1 rappresentato da d_e1 ). Tale equivalenza si otterrà tramite una trasformazione lineare che sarà ottenuta tramite una matrice che moltiplicata per i GdL globali ci da quelli locali. (figura in basso):

Tale matrice ha tante righe quanti i termini del vettore di GdL locali (24 nel nostro caso) e tante colonne quanti i termini del vettore di GdL globali (54 nel nostro caso). Tale matrice rappresenta un generico legame lineare che adesso andiamo a definire nello specifico. Consideriamo per esempio uno dei termini che abbiamo trattato : w_e1n2 (cerchiato in rosso nel vettore locale nella figura in alto). A questo termine nel vettore dei GdL locali, corrisponde una riga nella matrice di trasformazione (riga in rosso in figura). Poi, nel vettore dei GdL globali vado a isolare i termini che mi interessano, ovvero u_g2, v_g2 e w_g2 (cerchiati in blu), che tra l'altro sono gli unici che appaiono nella formula scritta prima. Questi 3 termini del vettore d_g occupano la settima, ottava e nona posizione, di conseguenza sulla matrice di trasformazione, ci posizione alla settima ottava e nona colonna (blu in figura) e incrociando la riga (rossa) citata prima, trovo ,in corrispondenza di quelle 3 caselle della matrice, proprio quei coefficienti ricavati anche prima (che sono i coseni degli angoli compresi tra il versore locale k_e1 e i 3 versori degli assi globali del nodo g2, ovvero i prodotti scalari tra quei versori) che moltiplicati per i vettori u_g2, v_g2 e w_g2 ci fanno ottenere w_e1n2.

Tutti gli altri contributi sono 0, ovvero in quella riga, avremo solo zeri tranne per i termini nella settima, ottava e nona colonna (e quei tre termini sono termini trigonometrici). Quindi ottengo, proprio come ottenuto prima che:

Quindi procedo così per tutti i GdL locali.

La matrice di trasformazione la chiamiamo MATRICE DI MAPPATURA da GdL globale a GdL locale, per l'elemento 1: $\underline{\underline{P}}$_e1.

Ovviamente è difficile rappresentare una matrice di 24 righe e 54 colonne, per cui la rappresentiamo nella seguente forma semplificata “a scala di grigi”:

Tale scala mostra il modulo di ognuno dei termini normalizzato da 0 a 1 (dato che sono termini trigonometrici). Il bianco è associato allo 0, il nero è associato all'1, grigio scuro è associato al coseno di 30° e grigio chiaro al seno di 30°. Tale scala lega i GdL locali per l'elemento 1 ai GdL globali della struttura. Si nota che è piena di zeri e i termini non nulli sono una frazione ridotta. Inoltre si nota che gli spostamenti vengono trasformati in spostamenti e le rotazioni in rotazioni. Inoltre ho termini non nulli solo in corrispondenza delle colonne associate ai nodi su cui effettivamente l'elemento 1 va a insistere (ovvero nella nostra struttura l'elemento 1 tocca i nodi globali 1 2 4 e 5, per cui i blocchi non nulli sono solo quelli associati a quei nodi globali).

Matrici analoghe si ottengono per ogni elemento, quindi in totale avremo 4 MATRICI di MAPPATURA (nel nostro caso).

Quindi se voglio scrivere il tutto in forma compatta, per ogni elemento posso scrivere:

Adesso, consideriamo una configurazione deformata a livello di struttura che è nota (vettore d_g) e tramite essa andiamo a derivare (tramite la matrice di mappatura) la configurazione deformata di un elemento ,per esempio dell'elemento 1 (vettore d_e1). L'elemento deformato reagisce elasticamente per cui, per mantenere l'elemento 1 nella configurazione d_e1, devo applicare delle forze che sono descritte da:

La matrice di rigidezza dell'elemento (che tra l'altro sarà la stessa per tutti gli elementi), $\underline{\underline{K}}$_e1 può essere rappresentata anch'essa in scala di grigi:

Bianco = 0 e Nero= termini con valore maggiore. Questa è una matrice lineare e simmetrica, in cui stiamo escludendo i moti di drilling, quindi possiamo escludere l'ultima riga e l'ultima colonna della tabella. Comunque, notiamo che, le forze entropiano sono funzione solo dei GdL di spostamento entropiano, perchè c'è disaccoppiamento tra moti membranali e moti flessotorsionali taglianti ( perchè c'è simmetria del materiale e ricordiamo che la simmetria porta a un disaccoppiamento dei due moti). Inoltre questa matrice è semidefinita positiva e si occupa di valutare il contributo energetico dei vari modi. In particolare esistono modi a contributo energetico nullo che sono (oltre a eventualmente i modi di drilling) i 6 moti di corpo rigido che si possono ottenere senza applicare nessuna forza. Quindi tale matrice trasforma i 6 moti di corpo rigido (che sono 6 vettori di modulo non nullo) in vettori forza nulli. Trasforma quindi, 6 vettori indipendenti in 6 zero, avrà quindi rango pari a n-6 (senza i moti di drilling n-10),quindi rango pari a 24-6=18 (14 escludendo il drilling). Quindi è una matrice ampiamente SINGOLARE.

Quindi in sostanza F_e1 sono le forze descritte secondo il sistema di riferimento locale all'elemento 1, che quell'elemento deve ricevere dall'esterno per far si che rimanga nella configurazione deformata d_e1 ed equilibrare quindi le proprie reazioni elastiche interne.

Adesso mi piacerebbe trasformare quell' F_e1 (che è in forma locale), in una equivalente globale. Quindi prendo $\underline{\underline{P}}$ che è la matrice di trasformazione da globale a locale e la inverto, ottenendo così la matrice inversa di trasformazione da locale a globale. Però se notiamo bene la matrice $\underline{\underline{P}}$ è una matrice rettangolare ,e ovviamente non esiste un inversa propriamente definita di una matrice rettangolare, ma esiste una PSEUDO-INVERSA. La matrice Pseudo-inversa più utilizzata è quella di Moore-Penrose. La matrice $\underline{\underline{P}}$ deriva da delle rotazione e delle permutazioni, e sia le rotazioni che le permutazioni sono delle trasformazioni ortogonali ,cioè l'inversa della matrice di trasformazione coincide con la trasposta. Quindi andiamo ad analizzare che proprietà ha la trasposta della matrice $\underline{\underline{P}}^{T}$_e1. Posso per esempio sviluppare i prodotti: 1) $\underline{\underline{P}}^{T}$_e1 $\underline{\underline{P}}$_e1 e 2) $\underline{\underline{P}}$_e1 $\underline{\underline{P}}^{T}$_e1.

Il primo prodotto dà luogo a una matrice quadrato con dimesioni pari alla piu elevata delle due (54×54 nel nostro caso. Il secondo prodotto dà luogo a una matrice sempre quadrata ma, con dimensione pari all più piccola delle 2 (24×24 nel nostro caso). Inoltre questo prodotto dà una matrice identità, quindi potrei dire che $\underline{\underline{P}}^{T}$_e1 è l'inverso di $\underline{\underline{P}}$_e1 (o quantomeno ci assomiglia, in quanto tutti i termini extradiagonali sono nulli, e i termini diagonali sono o 1 o 0 in base se i GdL globali in esame insistono o meno sull'elemento considerato). [Nel primo prodotto purtroppo non si riesce a tirar fuori la matrice identità perchè il rango della matrice identità dovrebbe avere al più la dimensione minore (mentre qui è proprio la dimensione maggiore).]

Quindi in pratica dato che consideriamo $\underline{\underline{P}}^{T}$_e1 l'inversa, possiamo quindi ricavare le forze sul sistema globale a partire dalle forze locali dell'elemento 1:

Dove F_ge1 è un vettore di 54 elementi raggruppati attraverso un raggruppamento Base-Nodo. Non sono altro che le forze sul sistema globale che devono essere assorbite dall'elemento 1 per mantenersi in equilibrio nella configurazione deformata d_g. Per cui, andando a sostituire otteniamo:

Ovviamente procedimenti analoghi devono essere fatti per tutti e 4 gli elementi. Definisco ora la matrice di rigidezza $\underline{\underline{K}}$_g,e1 associata all'elemento 1 e ai GdL globali:

Tale matrice $\underline{\underline{K}}$_g,e1 ha la seguente forma rappresentata in scala di grigi:

I termini nulli si trovano ovunque ci fossero degli zeri nella matrice $\underline{\underline{P}}$_e1, mentre i termini non nulli sono in corrispondenza degli incroci associati ai nodi globali 1,2,4,5 (ovvero quelli che insistono sull'elemento 1).

Quindi:

Quindi $\underline{\underline{K}}$_g,e1 moltiplicato per d_g mi da il contributo di forze esterne da applicare ai nodi dell'elemento per mantenerlo nella configurazione deformata descritta d_g. Però analogamente si devono applicare delle forze per mantenere l'elemento 2 nella configurazione deformata descritta d_g. E per fare ciò servirà $\underline{\underline{K}}$_g,e2 che sarà diversa dalla precedente $\underline{\underline{K}}$_g,e1 perchè sarà diversa la matrice di mappatura $\underline{\underline{P}}$_e2 rispetto a $\underline{\underline{P}}$_e1. La somma di $\underline{\underline{K}}$_g,e1 e $\underline{\underline{K}}$_g,e2 darà la seguente matrice:

Notiamo che sono stati aggiunti altri termini, in particolare quelli relativi all'elemento 2. Ora con considerazioni analoghe devo considerare l'elemento 3 e quindi la $\underline{\underline{K}}$_g,e3. Per cui la somma di $\underline{\underline{K}}$_g,e1, $\underline{\underline{K}}$_g,e2 e $\underline{\underline{K}}$_g,e3 darà una matrice con ulteriori termini aggiuntivi:

E infine considero nella somma anche $\underline{\underline{K}}$_g,e4 ottenendo quindi la matrice completa. Quindi semplicemente ottengo le forze che (data una generica configurazione deformata) si devono applicare sui nodi per equilibrare le reazioni elastiche di tutti gli elementi,nonchè le forze che devo applicare alla struttura per mantenerla nella generica configurazione deformata in esame:

Una volta sommati tutti i contributi di $\underline{\underline{K}}$ ottengo la seguente matrice che chiamo $\underline{\underline{K}}$_g (solo per distinguerla da quella di elemento, ma quando si parla di struttura sarà semplicemente $\underline{\underline{K}}$ ), che è la MATRICE DI RIGIDEZZA dell'intera struttura:

Tale matrice è somma di matrici simmetriche (perchè le matrici di rigidezza degli elementi sono di per sè simmetriche e poi le abbiamo pre e post moltiplicate rispettivamente per una matrice trasposta e per la stessa matrice non trasposta, per cui le $\underline{\underline{K}}$_g,ei sono simmetriche) perciò è simmetrica. Quindi considero solo la parte triangolare superiore e poi la parte inferiore la derivo per simmetria. Cosa importante è che oltre una certa distanza dalla diagonale, ci sono aree della matrice nel quale ci sono solo zeri. Questo è dovuto al fatto che un nodo non ha un ponte elastico diretto con tutti i nodi della struttura, per esempio il nodo 1 ha un ponte elastico diretto (il ponte sarebbe l'elemento) con i nodi 2,4 e 5 ma non con il nodo 9, per cui la casella tra nodo 1 e 9 sarà bianca per esempio. Ma perchè non si vogliono dei termini non nulli su quegli angoli? E in particolare, perchè il Marc quando calcola fa in modo da rendere il maggiore possibile l'area piena di zeri (avvicinandola di fatto alla diagonale)? Perchè così posso avere una matrice bandata simmetrica. Chiamo b (banda) la distanza (in direzione riga o in direzione colonna) della più lontana cella non nulla dalla diagonale, e voglio fare in modo che sia il più ridotta possibile. A che cosa mi serve che sia ridotta? Per ridurre la quantità di memoria che mi serve per calcolare la matrice, cioè invece che calcolare $N^{2}$ (dove N sono i GdL) termini, o meglio $N^{2}$/2 (dato che è simmetrica), se la matrice è bandata dovrò calcolarne solamente N·b. E molto spesso risulta che b « N. Un'altra questione è che se quella $\underline{\underline{K}}$_g diventa la matrice di un sistema di equazioni lineari, il costo computazionale della matrice non bandata è proporzionale a $N^{3}$ (che sia simmetrica o meno): α·$N^{3}$. Se invece la matrice è bandata il costo computazionale è proporzionale con diverso coefficiente di proporzionalità a N·$b^{2}$: β·N·$b^{2}$. Quindi se b è molto minore di N il costo computazionale cala vertiginosamente nel passare da matrice piena a matrice bandata. Quindi non importa solo il numero di GdL nel costo computazionale ma anche la banda, che tra l'altro ha anche maggiore influenza perchè è al quadrato in quella relazione.

dispensa_2018_al_2018-04-26.pdf

---