Ingegneria degli acquiferi - Appunti

Proprietà delle funzioni di interpolazione lagrangiane

A partire dalla funzione N, si possono osservare le seguenti proprietà:

• La funzione N vale 1 sul nodo 1 e zero su tutti i punti del segmento 2-3.
• La funzione N vale 1 sul nodo 2 e zero su tutti i punti del segmento 1-3.
• La funzione N vale 1 sul nodo 3 e zero su tutti i punti del segmento 1-2.

Le proprietà delle funzioni di interpolazione lagrangiane associate ai tre nodi dell'elemento generico V sono:

Nelle equazioni (13.49) e (13.50) si può osservare che:

N + N + N = 1

Si ipotizza di prendere nel dominio dell'elemento generico V la seguente funzione interpolante:

ĥ(x, t) = h(t)N(x)

Come conseguenza della condizione espressa dalla (13.49), risulta che nei nodi dell'elemento la funzione interpolante ĥ assume gli stessi valori della funzione vera h, si ha cioè dalla (13.49):

ĥ(x, x, t) = h(x, x, t) = h(t)

ĥ(x, x, t) = h(x, x, t) = h(t)

ĥ(x, x, t) = h(x, x, t) = h(t)

x , t) = h (t)³ dove h (t), h (t), h (t) sono i valori della variabile dipendente nei tre nodi 1, 2, 3 dell'elemento V al generico tempo t. Inoltre, il fatto che ognuna delle funzioni di interpolazione di un nodo del triangolo sia identicamente nulla sul segmento opposto al nodo considerato, assicura che la soluzione approssimata sia continua lungo i segmenti di unione di elementi contigui, in quanto il valore di ĥ su tale segmento è una combinazione lineare dei valori nodali delle funzioni di interpolazione associate ai due soli nodi del segmento comune ai due elementi contigui.

(Se si indica con ĥ) il residuo di (13.41) che viene ortogonalizzato rispetto alle funzioni di interpolazione associate ai nodi 3 dell'elemento V, si ha:

∫ eR L(ĥ)N dv = 0 n = 1, 2, 3 (13.53)

E dalla (13.51) e dalla (13.41) si ha pertanto:

∂ ĥ / ∂ ĥ∂Z Z Ze e e−(T )N dv bqN dv = S N dv (13.54)ij n n n∂x ∂x ∂te ee V VV i j

Applicando il teorema di Gauss al

primo termine di (13.54) e utilizzando la(13.51) si ottiene: e e e∂N ∂N ∂NZ Zm m ne − −T h N α ds T h dv (13.55)ij m ni ij mn∂x ∂x ∂xe eS Vj j i∂hZ Z me e ebqN dv = SN dvNn m n∂te eV VInserendo la condizione di flusso prefissato (13.43) nella (13.55) si ottieneinfine:13.3. MODELLO DI FLUSSO 2D CON EF GALERKIN 229#"Z e e∂N ∂N Z m n e edv h + SN N dv ḣ = (13.56)T m mij m n∂x ∂x ee VV j iZ Ze e− −bq̄N ds bqN dv n, m = 1, 2, 3n ne eS VIn forma matriciale la (13.56) può essere scritta:enm e e eA h + B ḣ = Q̄ + Q n, m = 1, 2, 3 (13.57)m mnm n nenm e ×dove A e B sono le matrici elementari (3 3) dell’ elemento genericonmeV , che prendono, rispettivamente il nome di matrice delle conduttanzaen en ×e matrice capacitiva, mentre Q̄ e Q sono i vettori elementari (3 1)che mettono in conto rispettivamente, le condizioni di flusso imposte su S 1ed il termine sorgente.Nel caso di

elementi triangolari a 3 nodi, gli integrali per il calcolo dei termini delle matrici e dei vettori elementari possono essere valutati in forma chiusa, utilizzando le seguenti due formule per il calcolo degli integrali di superficie e di linea che compaiono nel calcolo delle matrici elementari:

I_ijk = (N_i) (N_j) (N_k) dv = 2A (13.58)

I_3e1 = 0 0 A /3
I_3e2 = 0 0 A /6
I_3e1 = 1 0 A /12
I_3e0 = 1 1 A /12

Tabella 13.3: Valori dell'integrale I (13.60)
i cui valori sono riportati nella tabella (13.3) e dalla formula:

I_ij = (N_i) (N_j) ds = (13.59)

Dalla (13.45) si ha inoltre:
∂N_b / ∂N_ci = I_ii = 1, 2, 3 (13.60)
∂e_e / ∂x = 2A / ∂x = 2A1 2

Tenendo presenti la (13.58), la (13.59), la (13.60) dalla (13.56) si ha:
Z_m = T dv = (13.61)
Z_ij / ∂x / ∂x_e = V_j i1 (T_bb + T_bc + T_cb + T_cc) n, m = 1, 2, 3
n m 12 n m 21 n m 22

n me4A ( eSA /6 per n = mZe e eB = SN N dv = (13.62)nm m n e 6SA /12 per n = meV 1Zen e enmQ̄ = bq̄N ds = bq̄L (13.63)n 2eS 1Z een ebqN dv =Q = bqA (13.64)n 3eV
Ottenute le matrici elementari di tutti gli elementi appartenenti al dominio
di flusso globale V è possibile pervenire ad una equazione matriciale formal-
mente analoga alla (13.57) dove compaiono le matrici globali A , B , Q̄nm nm ne Q , ottenute dalle corrispondenti matrici elementari mediante la proce-
ndura di assemblaggio che consiste nel sommare per ogni nodo le equazioni
elementari di tutti gli elementi che contengono detto nodo. A livello globale
si avrà pertanto:
A h + B (13.65)ḣ = Q̄ + Q n, m = 1, 2 . . . , nnm m nm m n n t
che costituisce un sistema di n equazioni differenziali alle derivate ordinarie
tcon variabili dipendenti i carichi idraulici ai nodi.
Ora discretizzando anche il dominio del tempo, ossia facendo t = k∆t, dalla
(13.65), applicando la tecnica delle DF per eprimere la derivata rispetto al
tempo di

Equazioni di flusso 2D con EF Galerkin

Si ottiene:

k+1 k−h hm mk+1A h + B = Q̄ + Q n, m = 1, 2 . . . , n (13.66)

nm nm n n tm k∆t che può essere messa nella forma:

k+1 kD h = P n, m = 1, 2, . . . n (13.67)

nm tm n dove: BBnm nmk kn kn k+ Q +D = A + ; P = Q̄ h i, n = 1, 2, . . . nnm nm tn mk k∆t ∆t (13.68)

La risoluzione del sistema di n equazioni algebriche, dopo aver modificato in modo opportuno le righe della matrice dei coefficienti D e del vettore P per fare in modo da rispettare la condizione al limite di carico imposto (13.44) sui nodi del confine S, permette di calcolare il valore dei carichi idraulici al tempo t = (k + 1)∆t. Bisogna inoltre notare che in questo problema l'applicazione del teorema di Gauss sulla condizione di Galerkin di ortogonalizzazione del residuo, ha di fatto utilizzato la cosiddetta formulazione debole di Galerkin che ha come conseguenza l'inserimento automatico delle condizioni di Neuman espresse dalla (13.43).

Per meglio comprendere la procedura di assemblaggio che permette di sommare ad ogni nodo del dominio il contributo di tutti gli elementi che contengono il nodo che si considera, cioè di passare dalle matrici elementari alle matrici globali, si fa l'esempio elementi e 5 nodi come in Fig.(13.7). In questo semplice caso la tabella delle connessioni è riportata nella tabella(13.4).

Per ognuno dei tre elementi vale la relazione:

h_ij = e_iA_h = Q_e = 1, 2, 3; i, j = 1, 2, 3 (13.69)

Scrivendo le 9 equazioni espresse dalla (13.69), utilizzando per le variabili h gli indici globali che si rilivano dalla tabella delle connessioni (13.4), si ha:

111 112 113 1a_h + a_h + a_h = q (13.70)

1 4 2

1121 122 123 1a_h + a_h + a_h = q (13.71)

1 4 2

2131 132 133 1a_h + a_h + a_h = q (13.72)

1 4 2

3211 212 213 2a_h +

a_h + a_h = q (13.73)⁴ 3 2 1221 222 223 2

a_h + a_h + a_h = q (13.74)⁴ 3 2 2231 232 233 2

a_h + a_h + a_h = q (13.75)⁴ 3 2 3311 312 313 3

a_h + a_h + a_h = q (13.76)⁵ 3 4 1321 322 323 3

a_h + a_h + a_h = q (13.77)⁵ 3 4 2331 332 333 3

a_h + a_h + a_h = q (13.78)⁵ 3 4 3

Da questo sistema di 9 equazioni si passa al sistema di 4 equazioni globali del caso in esame sommando le equazioni (13.72), (13.75) relative al nodo 2, e le equazioni (13.74), (13.77), relative al nodo 3, le equazioni (13.71), (13.73), e (13.78) relative al nodo 4, cosı̀ facendo si ottiene che le equazioni per i nodi 1, 2, 3, 4, 5, sono nell'ordine:

111 113 112 1a_h + a_h + a_h = q (13.79)1 2 4 1131 133 233 232 132 231 1 2a_h + (a + a_h)h + a_h + (a + a_h)h = q + q1 2 3 4 3 3223 222 322 221 323 321 2 3a_h + (a + a_h)h + (a + a_h)h + a_h = q + q2 3 4 5 2 2121 123 213 212 332 122 211 333 331 1 2 3a_h + (a + a_h)h + (a + a_h)h + (a + a + a_h)h + a_h = q + q + q1 2 3 4 5 2 1 3312 313 311 3a_h + a_h + a_h = q3 4 5 3

La procedura di

do ie = 1, nt loop sugli elementi finiti (ef) (13.80)
    do i = 1, 3 loop sulle righe delle matrici (ef)
        ir = el(ie, i)
        eQ(ir) = Q(ir) + Q(i)
        do j = 1, 3 loop sulle colonne delle matrici (ef)
            ic = el(ie, j)
            eA(ir, ic) = A(ir, ic) + A(i, j)
            eB(ir, ic) = B(ir, ic) + B(i, j)
        enddo
    enddo
enddo

13.4. INTEGRAZIONE NEL TEMPO DI SISTEMI DI EQUAZIONI DIFFERENZIALI ORDINARIE 233
Sia dato un sistema di equazioni differenziali ordinarie, rappresentato in forma matriciale come segue: (13.81)
Ah + B ḣ = Q

1. Differenze in avanti