Utilizzando un elemento isoparametrico a 4 nodi, conviene introdurre un sistema di coordinate (o,ξ,η) detto sistema di coordinate naturali. Definite le funzioni di mappatura come:

$$\left\{\begin{matrix} x(\xi \eta ) =\sum_{i}^{4}N_{i}(\xi \eta )x_{i}\\ y(\xi \eta ) =\sum_{i}^{4}N_{i}(\xi \eta )y_{i} \end{matrix}\right. \qquad(1)$$

dove (x,y) sono le coordinate del sistema fisico, $N_{i}(ξ,η)$ le funzione di forma e $x_{i} y_{i}$ le posizioni dei nodi nel sistema di riferimento fisico. Le funzioni di mappatura possono essere esplicitate in modo da ottenere una espressione del tipo:

$$\left\{\begin{matrix} x=\alpha \xi +\beta \eta +\gamma \xi \eta +\delta\\ y=\alpha '\xi +\beta' \eta +\gamma' \xi \eta +\delta' \end{matrix}\right. \qquad(2)$$

in cui $ \alpha \beta \gamma $ sono costanti.

Tali espressioni ci suggeriscono che tali funzioni seguono un comportamento lineare lungo sezioni ortogonali agli assi (ξ,η) del sistema naturale. Infatti, se considero il bordo tra i nodi 1 e 2 nel riferimento naturale, essendo un segmento ad η costante, viene trasformato linearmente ed i segmenti mostrati in figura sono perciò proporzionali.

$$\frac{a}{b}=\frac{a'}{b'} \qquad(3)$$

Lo stesso si può dire per tutti i bordi dell'elemento quadrato.

Un punto medio del bordo rimarrà dunque un punto medio del bordo anche nel riferimento fisico. Per le stesse ragioni anche gli assi ζ ed η vengono trasformati linearmente e conservano dunque la loro rettilineità. Sapendo inoltre che questi intersecano i bordi dell'elemento deformato nei punti medi dei bordi diventa immediato tracciare gli assi delle coordinate naturali anche nel sistema fisico.

$$Figura 1$$

Ragionando in modo analogo è ancora possibile ricavare per via geometrica la posizione di un punto P dal sistema di riferimento naturale a quello fisico. Individuato P infatti come intersezione tra due rette parallele agli assi ζ ed η nel riferimento naturale esso corrisponderà all'intersezione delle stesse rette trasformate linearmente nel riferimento fisico. Tali rette sono individuate come rette passanti per due punti appartenenti al bordo che vengono ricavati tramite una relazione di proporzionalità analoga alla precedente.

$$\frac{c}{d}=\frac{c'}{d'}$$

$$Figura 2$$

Si considera un triangolo rettangolo sul piano naturale con cateti paralleli agli assi $(\xi ,\eta)$ di lunghezza infinitesima $d\xi , d\eta $come mostrato in Fig.2. L’area del triangolo risulta:

$$A_{\xi \eta }=\frac{1}{2}d\xi d\eta \qquad(4)$$

I nodi del triangolo nel sistema naturale hanno coordinate :

$Nodo_{1}=\xi ,\eta$

$Nodo_{2}=\xi+d\xi ,\eta \qquad(5)$

$Nodo_{3}=\xi ,\eta +d\eta $

Nel piano fisico , con le funzioni di mappatura le coordinate dei nodi del triangolo sono :

$Nodo_{1}=x(\xi,\eta),y(\xi,\eta)$

$Nodo_{2}=x(\xi+d\xi,\eta),y(\xi+d\xi,\eta)\qquad(6)$

$Nodo_{3}=x(\xi,\eta+d\eta),y(\xi,\eta+d\eta)$

Consideriamo come esempio :

$x(\xi ,\eta )=\sum_{i}^{3}N_{1}(\xi ,\eta )x_{i}$ $$\qquad(7)$$ $x(\xi+d\xi ,\eta )=\sum_{i}^{3}N_{1}(\xi+d\xi ,\eta )x_{i}$

Con l’ipotesi che i lati del triangolo nel piano naturale siano piccoli si ha:

$d\xi ,d\eta < < 1$

Sviluppando con Taylor le coordinate dei punti nel piano fisico (6), queste risultano essere

$Nodo_{1}=x(\xi ,\eta ),(\xi ,\eta )$

$Nodo_{2}=x(\xi ,\eta )+\frac{\delta x}{\delta \xi }|(\xi \eta )d\xi ,y(\xi \eta )+\frac{\delta y}{\delta \xi }|(\xi \eta )d\xi\qquad(8)$

$Nodo_{3}=x(\xi ,\eta )+\frac{\delta x}{\delta\eta}|(\xi \eta)d\eta,y(\xi \eta)+\frac{\delta y}{\delta\eta}|(\xi\eta)d\eta$

Con queste coordinate è possibile calcolare l’area del triangolo nel piano fisico :

$$A_{xy}=\frac{1}{2!} \begin{bmatrix} 1 & 1 & 1\\ x & x+\frac{\delta x}{\delta \xi }d\xi & x+\frac{\delta x}{\delta \xi }d\eta\\ y & y+\frac{\delta y}{\delta \xi }d\xi & y+\frac{\delta y}{\delta \eta }d\eta \end{bmatrix}\qquad(9)$$

Sottraendo alla seconda riga la prima moltiplicata per x [ R2 / (R2 - x*R1)] e alla terza riga la prima moltiplicata per y [ R3 / (R3 - y*R1)] si ottiene: $$A_{xy}=\frac{1}{2}\begin{vmatrix} 1 &1 &1 \\ 0&\frac{\delta x}{\delta \xi }d\xi &\frac{\delta x}{\delta \eta }d\eta \\ 0& \frac{\delta y}{\delta \xi }d\xi &\frac{\delta y}{\delta \eta }d\eta \end{vmatrix}=\frac{1}{2}(-1)^{2}\begin{bmatrix} \frac{\delta x}{\delta \xi }d\xi &\frac{\delta x}{\delta \eta }d\eta \\ \frac{\delta y}{\delta \xi }d\xi & \frac{\delta y}{\delta \eta }d\eta \end{bmatrix}\qquad(10)$$

Essendo $(d\xi , d\eta)$due scalari , ed essendo presenti indistintamente su tutti i termini di una colonna, è possibile scrivere:

$$A_{xy}=\frac{1}{2}d\xi d\eta\begin{bmatrix} \frac{\delta x}{\delta \xi } &\frac{\delta x}{\delta \eta } \\ \frac{\delta y}{\delta \xi } & \frac{\delta y}{\delta \eta } \end{bmatrix}\qquad(11)$$

Dove:

$\frac{1}{2}d\xi d\eta=A_{\xi \eta }$

e

$$\begin{bmatrix} \frac{\delta x}{\delta \xi } &\frac{\delta x}{\delta \eta } \\ \frac{\delta y}{\delta \xi } & \frac{\delta y}{\delta \eta } \end{bmatrix}\qquad(12)$$

rappresenta la matrice Jacobiana trasposta , il cui determinante chiameremo Jacobiano

Per quanto riguarda le deformazioni, queste possono essere scritte come

$$\underline{\varepsilon }\begin{bmatrix} \varepsilon _{x}\\ \varepsilon _{y}\\ \varepsilon _{xy} \end{bmatrix} =\begin{bmatrix} 1 &0 &0 &0 \\ 0&0 &0 &1 \\ 0 &1 &1 &0 \end{bmatrix} \begin{bmatrix} \frac{\delta u}{\delta x}\\ \frac{\delta u}{\delta y}\\ \frac{\delta v}{\delta x}\\ \frac{\delta v}{\delta y} \end{bmatrix}\qquad(13)$$

La matrice 3×4 viene indicata come $\underline{H}$,mentre ci si concentra sulla scrittura del vettore a 4 componenti, a destra nella (13) , in funzione del sistema naturale. Come già visto,considerando momentaneamente, solamente lo spostamento $\underline{u}$:

$$\begin{bmatrix} \frac{\delta u}{\delta \xi }\\ \frac{\delta u}{\delta \eta } \end{bmatrix} =\underline{J}\begin{bmatrix} \frac{\delta u}{\delta x}\\\ \frac{\delta u}{\delta y }\ \end{bmatrix}\qquad(14)$$

Supponendo $\underline{J}$ invertibile , ossia$\left | \underline{J} \right |\neq 0$ per ogni valore , è lecito scrivere

$$\begin{bmatrix} \frac{\delta u}{\delta x}\\ \frac{\delta u}{\delta y} \end{bmatrix} =\underline{J}^{-1} \begin{bmatrix} \frac{\delta u}{\delta \xi }\\ \frac{\delta u}{\delta \eta } \end{bmatrix}\qquad(15)$$

$$ \begin{bmatrix} \frac{\delta v}{\delta x}\\ \frac{\delta v}{\delta y} \end{bmatrix} =\underline{J}^{-1} \begin{bmatrix} \frac{\delta v}{\delta \xi }\\ \frac{\delta v}{\delta \eta } \end{bmatrix}\qquad(16)$$

considerando che:

$$\underline{J}^{-1} =\frac{1}{\left | \underline{J} \right |} \begin{bmatrix} J_{22} & -J_{12}\\ -J_{21} &J_{11} \end{bmatrix} $$

Utilizzando la matrice $\underline{J}^{-1}$ si può inserire il vettore a 4 componenti trovato nella (13)

$$\begin{bmatrix} \frac{\delta u}{\delta x}\\ \frac{\delta u}{\delta y}\\ \frac{\delta v}{\delta x}\\ \frac{\delta v}{\delta y} \end{bmatrix} = \begin{bmatrix} \underline{J_{11}^{-1}} &\underline{J_{12}^{-1}} &0 &0 \\ \underline{J_{21}^{-1}} &\underline{J_{22}^{-1}} &0 &0 \\ 0& 0 & \underline{J_{11}^{-1}} &\underline{J_{12}^{-1}} \\ 0 & 0& \underline{J_{21}^{-1}} &\underline{J_{22}^{-1}} \end{bmatrix} \begin{bmatrix} \frac{\delta u}{\delta \xi}\\ \frac{\delta u}{\delta \eta}\\ \frac{\delta v}{\delta \xi}\\ \frac{\delta v}{\delta \eta} \end{bmatrix}\qquad(17)$$

La matrice 4×4 trovata nell'equazione (17) si indica con $J_{inv}^{*}$

Per esplicitare il vettore a 4 componenti a destra nulla (17) attraverso le funzioni di forma, si ricava che:

$\frac{\delta u}{\delta\xi}=\sum_{i}^{4}\frac{\delta N_{i}}{\delta \xi }u_{i}$

$\frac{\delta u}{\delta\eta}=\sum_{i}^{4}\frac{\delta N_{i}}{\delta \eta }u_{i}$

$\frac{\delta v}{\delta\xi}=\sum_{i}^{4}\frac{\delta N_{i}}{\delta \xi }v_{i}$

$\frac{\delta v}{\delta\eta}=\sum_{i}^{4}\frac{\delta N_{i}}{\delta \eta }v_{i}$

con: $$\underline{N}=\begin{bmatrix} N_{1} &0 &N_{2} &0 &N_{3} &0 &N_{4} &0 \\ 0 &N_{1} &0 &N_{2} &0 &N_{3} &0 &N_{4} \end{bmatrix}$$

Il vettore sopra citato nella (17) può essere scritto in forma matriciale come:

$$\begin{bmatrix} \frac{\delta u}{\delta \xi }\\ \frac{\delta u}{\delta \eta }\\ \frac{\delta v}{\delta \xi }\\ \frac{\delta v}{\delta \eta } \end{bmatrix} = \underline{Q}(\xi ,\eta) \begin{bmatrix} u_{1}\\ v_{1}\\ u_{2}\\ v_{2}\\ u_{3}\\ v_{3}\\ u_{4}\\ v_{4} \end{bmatrix}\qquad(18)$$

Il vettore a destra sulla (18) si indica come $\underline{\delta }$ e $\underline{Q}(\xi ,\eta )$

ed è una matrice 4×8 così descritta:

$$\underline{Q}(\xi ,\eta )= \begin{bmatrix} \frac{\delta N_{1}}{\delta \xi } &0 &\frac{\delta N_{2}}{\delta \xi } &0 &\frac{\delta N_{3}}{\delta \xi } &0 &\frac{\delta N_{4}}{\delta \xi } &0 \\ \frac{\delta N_{1}}{\delta \eta } &0 &\frac{\delta N_{2}}{\delta \eta } &0 &\frac{\delta N_{3}}{\delta \eta } &0 &\frac{\delta N_{4}}{\delta \eta } &0 \\ 0&\frac{\delta N_{1}}{\delta \xi }& 0 & \frac{\delta N_{2}}{\delta \xi } &0 &\frac{\delta N_{3}}{\delta \xi } &0 &\frac{\delta N_{4}}{\delta \xi } \\ 0& \frac{\delta N_{1}}{\delta \eta }& 0 & \frac{\delta N_{2}}{\delta \eta } &0 &\frac{\delta N_{3}}{\delta \eta } &0 &\frac{\delta N_{4}}{\delta \eta } \end{bmatrix}$$

Possono essere descritte a questo punto le deformazioni nel sistema naturale come:

$\underline{\varepsilon }=\underline{H}\cdot \underline{J_{inv}^{*}}(\xi ,\eta )\cdot \underline{Q}(\xi ,\eta )\cdot \underline{\delta }\qquad(19)$

L’energia potenziale elastica dell’elemento isoparametrico a 4 nodi (U) può essere scritta come:

$U=\iint_{a_{xy}}^{ }\frac{1}{2}\underline{\varepsilon }^{T}\cdot \underline{\sigma }\cdot da_{xy}\qquad(20)$

Ricordando che $a_{xy}=a_{\xi \eta }\left | J \right |$ la (20) può essere scritta come:

$U=\iint_{a_{\xi \eta }}^{ }\frac{1}{2}\underline{\varepsilon }^{T}\cdot \underline{\sigma }\cdot \left | J(\xi \eta ) \right |da_{\xi \eta }\qquad(21)$

Utilizzando la (19) ed introducendo la matrice D del legame costitutivo che lega e tensioni alle deformazioni ,tali che:

$\underline{\sigma }=\underline{D}\cdot \underline{\varepsilon }\qquad(22)$

La (21) può essere scritta come:

$U=\frac{1}{2}\int_{-1}^{1}\int_{-1}^{1}\underline{\delta }^{T}\cdot \underline{Q}^{T}\cdot \underline{J_{inv}^{*T}}\cdot \underline{H}^{T}\cdot \underline{D}\cdot \underline{H}\underline{J_{inv}^{*}}\cdot \underline{Q}\cdot \underline{\delta }\cdot \left | \underline{J} \right |d\xi d\eta \qquad(23)$

Ponendo

$\underline{H}\cdot \underline{J_{inv}^{*}}\cdot \underline{Q}\cdot \underline{\delta }=\underline{B}(\xi \eta )\qquad(24)$

la (23) può essere scritta come :

$U=\frac{1}{2}\underline{\delta }^{T}\left ( \int_{-1}^{1}\int_{-1}^{1} \underline{B}^{T}\cdot \underline{D}\cdot \underline{B}\left | \underline{J} \right |d\xi d\eta \right )\underline{\delta }\qquad(25)$

(Si è potuto portare fuori $\underline{\delta }^{T}$ ed anche $\underline{\delta }$ in quanto non dipendenti dalle variabili di integrazione $(\xi,\eta)$.)

Ponendo:

$\underline{K}=\int_{-1}^{1}\int_{-1}^{1} \underline{B}^{T}\cdot \underline{D}\cdot \underline{B}\left | \underline{J} \right |d\xi d\eta$

dove K rappresenta la matrice rigidezza dell’elemento, la (25) si scriverà come:

$U=\frac{1}{2}\underline{\delta }^{T}\cdot \underline{K}\cdot\underline{\delta }\qquad(26)$

A differenza dal caso relativo all'elemento triangolare, l'integrale contenuto nell'espressione (26), esplicitato nella (25), non si è in grado di risolverlo analiticamente. Viene perciò utilizzato il il metodo della quadratura di Gauss secondo il quale:

$\underline{K}=\sum_{l=1}^{4}W_{l}\underline{B}^{T}(\xi_{l} \eta_{l} )\cdot \underline{D_{\left [ l \right ]}}\cdot \underline{B}(\xi_{l} \eta_{l} )$

dove i punti candidati ad essere punti di integrazione di Gauss sono:

$P_{1}=(\frac{1}{\sqrt{3}};\frac{1}{\sqrt{3}})$

$P_{2}=(\frac{1}{\sqrt{3}};-\frac{1}{\sqrt{3}})$

$P_{3}=(-\frac{1}{\sqrt{3}};-\frac{1}{\sqrt{3}})$

$P_{4}=(-\frac{1}{\sqrt{3}};\frac{1}{\sqrt{3}})$

$$Figura 3$$

$W_{i}$ è il peso specifico dell'l-esimo punto di Gauss (in genere $W_{i}$=1). La matrice $\underline{D}$ è la matrice che esprime il legame costitutivo e può essere eventualmente funzione dei punti di Gauss; si può dunque pensare che un elemento isoparametrico a 4 nodi possa avere globalmente un comportamento elastoplastico espresso mediante i propri punti di Gauss che singolarmente possono avere differenti matrici $\underline{D}$. Si preferisce, tuttavia, evitare elementi parzialmente plasticizzati sfruttando elementi triangolari aventi un unico nodo e un'unica matrice $\underline{D}$ associata.

Il software Mark utilizza due modalità di integrazione gaussiana a seconda del numero di punti di Gauss sfruttati: Full Integration e Reduced Integration.

Full Integration: si lavora con 4 punti di Gauss, numero necessario e sufficiente per integrare esattamente la matrice rigidezza per elementi aventi Jacobiano costante (elementi indeformati).

Reduced Integration: non si lavora mai con punti di Gauss sufficienti per risolvere esattamente l'integrale ma compare sempre un errore di calcolo. Tuttavia a volte la matrice $\underline{K}$ inesatta risulta migliore di quella calcolata con il metodo Full Integration.

Un caso in cui si ha questa peculiarità è quello relativo al fenomeno dello Shear Locking analizzato di seguito.

Lo Shear Locking è un fenomeno che si manifesta nella teoria dell'elemento isoparametrico a 4 nodi, poiché si ammette che la rettilineità dei bordi venga mantenuta prima e dopo la trasformazione nell'elemento finito, a discapito dell'esatta deformata dell'elemento. Si prenda come esempio un elemento di trave soggetto a flessione pura. La fig.4 mostra la soluzione esatta delle deformata dell'elemento di trave.

$$Figura 4$$

D'altro canto, affinché vengano mantenuti i bordi deformati rettilinei, tramite la teoria degli elementi finiti isoparametrici a 4 nodi si manifesta una deformazione tangenziale che prima era assente.

$$Figura 5$$

La figura successiva mostra invece un concio in cui viene applicata la sola deformazione angolare comparsa nella deformazione a trapezio $\gamma _{xy}$

Da cui a livello energetico ne deriva:

$E_{anodi}=E_{esatta}+E_{\gamma }$

ossia: l'energia associata alla deformata dell'elemento isoparametrico a 4 nodi $E_{4nodi}$ è la somma dell'energia associata alla deformata esatta $E_{esatta}$ e quella associata alla pura deformazione angolare $E_{\gamma }$. In tal modo $E_{\gamma }$ risulta l'errore che fa divergere la soluzione esatta da quella ottenuta con la teoria degli elementi finiti. Accade infatti che all'elemento di trave sia applicata una coppia maggiore di quella associata al caso esatto affinché si ottenga la deformata a trapezio.

$$Figura 6$$

Ora, con una ipotetica Full Integration secondo 4 punti di Gauss come in Fig.6, vengono prese in considerazione inevitabilmente le $\tau _{xy}$ provocanti quell'incremento del 48% dell'energia potenziale elastica che si sarebbe voluto evitare. Utilizzando, invece, un unico punto di Gauss posto nel centroide ,con una Reduced Integration, è immediato intuire come nel calcolo dell'energia potenziale non venga aggiunto il surplus energetico proprio della deformazione a taglio. Tuttavia se da una parte l'utilizzo di un unico punto di Gauss elimina il fenomeno dello Share Locking, non si ottiene comunque un risultato che si possa dire accettabile, infatti, il problema è che il centroide non risulta immune solo alla deformazione a taglio ma anche all'effetto degli sforzi normali e dunque l'energia dell'elemento sarà paradossalmente nulla.