$$ \newcommand{\vec}[1]{\smash{\underline{#1}}} \newcommand{\mat}[1]{\smash{\underline{\underline{#1}}}} $$
Prendendo in esame la stessa struttura si illustra un esempio di costruzione della “Matrice di rigidezza vincolata”. Nel caso specifico si considerano vincolati i nodi “1” e “4” rispettivamente con carrello e cerniera, corrispondenti ai gradi di libertà “2”, “7” e “8”.
La struttura è a elementi triangolari, caricata sul nodo 2 da $Px$ e $Py$. Possiamo associare a ogni elemento triangolare una matrice rigidezza 6×6; gli elementi della matrice dell'elemento 1,2 e 3 sono rispettivamente $a_{ij}$ , $b_{ij}$ , $c_{ij}$. I contributi delle matrici di rigidezza di ogni elemento possono essere assemblati nella matrice rigidezza dell'intera struttura; le matrici rigidezza di ogni elemento possono essere costruite mediante equazioni di equilibrio ad ogni nodo dell'elemento a cui appartiene. Vediamo il procedimento per determinare la matrice rigidezza dell'elemento 1. I nodi dell'elemento 1 sono 4,5,2 a cui associamo una nomenclatura locale rispettivamente $i,j,k$. Per ogni nodo scriviamo l'equazioni di equilibrio lungo l'asse $x$ e l'asse $y$.
La prima riga contiene i contributi alle reazioni elastiche sul nodo 4. Le righe della matrice dell'elemento sono legate alle equazioni di equilibrio, le colonne invece sono associate a specifiche incognite di spostamento $u$. Per trasportare gli elementi dalla matrice locale alla matrice globale possiamo rifarci a una regola pratica: Chiamando “l” il nodo locale, il nodo globale sarà:
Ad esempio la 1° riga della matrice di rigidezza dell'elemento 1 andrà a finire nella 7° riga della matrice di rigidezza globale, la 2° riga alla 8° della matrice globale e così via. Analogamente si procede per le colonne
Ripetendo il procedimento per gli altri 2 elementi, otteniamo la matrice di rigidezza globale andando a sostituire e quindi sommare i termini $a_{ij}, b_{ij}, c_{ij}$ nella loro posizione globale. Notiamo che gli elementi diagonali delle matrici locali rimangono sulla diagonale della matrice globale. La matrice di rigidezza globale è una matrice singolare e simmetrica. La matrice globale si presenta nella seguente forma:
Introduciamo i vincoli per il nodo 4, i quali devono essere considerati uno per volta.
Gli spostamenti imposti sono:$ u_4=0, v_4=0,v_1=0 $. Partiamo dallo spostamento in direzione x, procedo cancellando il vincolo e inserendo un vincolo cinematico ossia: $u_4=\bar{u}_4$ Tolto così il vincolo ho imposto un vincolo cinematico e quindi considerando nello specifico il grado di libertà “7”, si annullano tutti gli elementi della settima riga tranne quello diagonale posto uguale a “1”. In questo modo senza togliere il carattere di incognita ad $u_4$, risulta essere pari a zero indipendentemente dagli altri valori.
Si esegue la medesima operazione sulla colonna “7” (essendo tutti termini noti) collocando gli elementi della stessa nel vettore dei termini noti, in questo modo si recupera la simmetria della matrice.
Si ottiene così il sistema matriciale vincolato relativo al grado di libertà “7”.
Si effettuano le stesse operazioni per i gradi di libertà “2” e “8” relativi ai restanti vincoli scelti ( $v_4=\bar{v}_4$, $v_1=\bar{v}_1$)
La matrice ottenuta è la “Matrice rigidezza vincolata”. Al termine noto è presente una quota parte delle reazioni elastiche degli elementi associati a uno spostamento. La matrice di rigidezza si può risolvere mediante un qualsiasi metodo di calcolo numerico.
Prendiamo un sistema previncolamento
$$ \mat{K} \vec{\delta} = \vec{P} $$ Consideriamo il sistema vincolato del tipo: $$ \mat{K}^v \vec{\delta}^* = \vec{P}^v $$
La soluzione avrà forma del tipo $ \vec{δ}^* = \mat{[K^v]}^{-1} \vec{P}^v $. Perciò essendo $$ \mat{K}\vec{\delta}^* \neq \vec{P} $$ possiamo scrivere che $$ \mat{K}\vec{\delta}^* = \vec{P} + \vec{R} $$ dove $\vec{R}$ è vettore reazione vincolare. La differenza tra le reazioni elastiche e le azioni esterne non sarà mai nulla poichè raccolgono l'errore numerico di calcolo. Dagli spostamenti nodali della struttura otteniamo gli spostamenti nodali elemento per elemento da cui ricaviamo le deformazioni e le tensioni.
Consideriamo questo elemento quadrilatero liberamente distorto nello spazio fisico $x,y$ (il corrispettivo nel 3D è il tetraedro a 8 nodi).
Costruisco le funzioni di interpolazioni sullo spazio fisico come: $$ u(x,y)=ax+by+cxy+d $$ dove $u$ è lo spostamento
La funzione di interpolazione può essere espressa anche in forma lineare del tipo $ u(x,y)=ax+by+d+cx^2 $ oppure $ u(x,y)=ax+by+d+cy^2 $
Dato che i termini non possono essere più di quattro devo scegliere una delle due forme, per cui non compaiono entrambi i termini al quadrato. Avrò un sistema del tipo $$ \left\{\begin{matrix} u(x_i,y_i)= u_i \\u(x_j,y_j)= u_j \\ u(x_k,y_k)= u_k \\ u(x_l,y_l)= u_l \end{matrix}\right. $$
Queste 4 equazioni ci permettono di definire i coefficienti $ a,b,c,d $. Procedendo in questa maniera tuttavia riusciamo a definire correttamente solo elementi rettangolari non distorti. Infatti non si rispetta la continuità tra lato e lato e si ha la generazione di tagli e sovrappozione di materiale.
Abbandoniamo questo tipo di approccio. Procederemo con un altra procedura più generale e generalizzabile. Per costruire la teoria dell'elemento isoparametrico a 4 nodi occorre affiancare al piano fisico un sistema di coordinate naturali e perciò un sistema naturale $ \xi $ e $ \eta $
Nel sistema di coordinate naturali l'elemento è sempre un quadrato che si estende su $ \xi $ e $ \eta $ da -1 a +1. Per ogni punto del piano naturale esiste un corrispondente sul piano fisico.
Definisco perciò le funzioni di mappatura: $ x(\xi , \eta) $ $ y(\xi ,\eta) $
Otteniamo così l'associazione di ogni punto del piano naturale alla sua immagine sul piano fisico. Le relazioni scritte sono biunivoche perciò possiamo definire anche una mappatura inversa.
Definisco le funzioni come combinazione lineare delle coordinate dei nodi(x,y) sul piano fisico pesate mediante le funzioni di forma. Scrivo la funzione di forma al nodo 1: $$ N_1(\xi,\eta)=a\xi + b\eta +c\xi\eta+d $$
L'elemento risulta non distorcersi mai sul paino naturale. Imponendo dei valori sui 4 nodi ottengo i coefficienti $ a,b,c,d $.
Consideriamo il nodo 1. Sapendo che le funzioni di forma valgono 1 sul nodo stesso e 0 sui restanti ottengo che: $$ \left\{\begin{matrix} N_1(-1,-1)=1 \\ N_1(1,-1)=0 \\ N_1(1,1)=0 \\ N_1(-1,1)=0 \end{matrix}\right. $$
Risolvendo il sistema ottengo $N_1$: $$ N_1= \frac{1}{4} (1-\xi)(1-\eta) $$
Possiamo perciò definire così le funzioni di forma: $$ N_{1,2,3,4}= \frac{1}{4} (1\pm \xi)(1\pm \eta) $$
Correggendo il segno in base alla caratteristica della funzione di essere unitaria sul nodo stesso e nulla sugli altri.
Consideriamo la funzione di forma al nodo 1
Nel sistema di coordinate naturali almeno uno tra \xi e \eta rimane costante. In generale la funzione di forma è bilineare ma fissando uno tra $\xi$ e $\eta$ si può considerare lineare sui 4 lati nel sistema di coordinate naturali. Rappresentiamo perciò l'andamento della funzione di forma $ N_1 $:
Possiamo trovare il valore della funzione di forma in un punto generico:
Si può notare come ognuno dei nodi viene mappato in se stesso.
$$ \left\{\begin{matrix} x(\xi,\eta)= N_1(\xi, \eta)x_1 +N_2(\xi, \eta)x_2 +N_3(\xi, \eta)x_3 +N_4(\xi, \eta)x_4 \\ y(\xi,\eta)= N_1(\xi, \eta)y_1 +N_2(\xi, \eta)y_2 +N_3(\xi, \eta)y_3 +N_4(\xi, \eta)y_4 \end{matrix}\right. $$
In generale risulta che la mappatura è di semplice stesura solo dal piano naturale a quello fisico, il contrario invece è più complesso. Le funzioni di forma sono anche utili per mappare gli spostamenti $$ \left\{\begin{matrix} u(\xi,\eta)= N_1(\xi, \eta)u_1 +N_2(\xi, \eta)u_2 +N_3(\xi, \eta)u_3 +N_4(\xi, \eta)u_4 \\ v(\xi,\eta)= N_1(\xi, \eta)v_1 +N_2(\xi, \eta)v_2 +N_3(\xi, \eta)v_3 +N_4(\xi, \eta)v_4 \end{matrix}\right. $$ dove $ u_i $ è lo spostamento lungo $x$ del nodo i-esimo, e $ v_i $ è lo spostamento lungo $y$ del nodo i-esimo.
In generale possiamo dire che se uno più funzioni di forma per definire le coordinate $(x,y)$ rispetto a quelle usate per gli spostamenti $(u,v)$ ho elementi detti superparametrici, nel caso inverso invece elementi subparametrici. Gli elementi isoparametrici sono elementi per cui l'interpolazione delle coordinate nodali per passare dal sistema naturale(dove l'elemento è non distorto) al sistema fisico (dove l'elemento è geometricamente distorto) e l'interpolazione degli spostamenti da sistema naturale a sistema fisico si ottengono tramite le stesse funzioni di forma.
Ottenuti gli spostamenti posso ricavare le deformazioni, tramite queste le tensioni. Attraverso tensioni e deformazioni posso ricavare l'energia potenziale elastica interna e la matrice rigidezza. Tuttavia le deformazioni sono definite sul piano naturale e non su quello fisico. Si vogliono calcolare le deformazioni generalizzate $\epsilon_x , \epsilon_y, \gamma_{x,y}$ a partire dagli spostamenti globali $u$ e $v$:
$ \epsilon_x=\frac{\partial u(\xi,\eta)}{\partial x} $
$ \epsilon_y=\frac{\partial u(\xi,\eta)}{\partial y} $
$ \gamma_{xy}=\frac{ \partial u(\xi,\eta) }{ \partial x } +\frac{\partial v(\xi,\eta)}{ \partial y } $
La deformazione in direzione $ x $ è: $$ \epsilon_x=\frac{\partial u(\xi,\eta)}{\partial x}=\frac{\partial u}{\partial \xi}\frac{\partial \xi}{\partial x}+\frac{\partial u}{\partial \eta}\frac{\partial \eta}{\partial x} $$
Il problema è che le derivate $\frac{\partial \xi}{\partial x}, \frac{\partial \eta}{\partial x} $ non sono facili da calcolare in quanto è disponibile solo x in funzione di $\xi$ e $\eta$, non il viceversa. Si decide quindi di utilizzare il seguente approccio solo sugli spostamenti $u$, quelli per $v$ sono analoghi.
$$ \left\{\begin{matrix} \frac{\partial u}{\partial \xi}=\frac{\partial u}{\partial x}\frac{\partial x}{\partial \xi}+\frac{\partial u}{\partial y}\frac{\partial y}{\partial \xi} \\ \frac{\partial u}{\partial \eta}=\frac{\partial u}{\partial x}\frac{\partial x}{\partial \eta}+\frac{\partial u}{\partial y}\frac{\partial y}{\partial \eta} \end{matrix}\right. $$
dove i termini $\frac{\partial u}{\partial \xi} ,\frac{\partial x}{\partial \xi}, \frac{\partial y}{\partial \xi},\frac{\partial u}{\partial \eta} ,\frac{\partial x}{\partial \eta}, \frac{\partial y}{\partial \eta} $ sono di facile soluzione. Il sistema è di due equazioni in due incognite e si può scrivere in forma matriciale: $$ \begin{bmatrix} \frac{\partial u}{\partial \xi} \\\frac{\partial u }{\partial \eta} \end{bmatrix} = \begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} \\\frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} \end{bmatrix} \cdot \begin{bmatrix} \frac{\partial u}{\partial x} \\\frac{\partial u}{\partial y} \end{bmatrix} $$
La matrice 2×2 è la matrice Jacobiana della trasformazione $\mat{J}(\xi,\eta)$.