Metodo Ai Volumi Finiti Per I Fenomeni Di Diffusione

Considriamo un volume $\Omega\subset \mathbb{R}^N$ e osserviamo che data la densità di particelle in un dato istante temporale che si trovano in un punto dello spazio, i.e $\rho=\rho(t,\vec{x})$, possiamo mettere in relazione la variazione di quest'ultima con il flusso delle particelle $\vec{J}(t,\vec{x})$:

$$ \partial_t \int_\Omega \rho(t,\vec{x}) d\vec{x} = -\int_{\partial\Omega} \vec{J}(t,\vec{x})\cdot \vec{n}$$

quest'ultima equazione integrale infatti ci dice che la virazione di massa all'interno del volume $\Omega$ devo corrispondere a un flusso di particelle al di fuori del volume ( all'interno). Supponiamo ora di lavoare con funzioni lisce, allora possiamo manipolare l'equazione precedente per ottenere la seguente formula, utilizzando il teorema di Gauss:

$$\int_\Omega \partial_t \rho(t,x) + \nabla\cdot \vec{J}(t,\vec{x}) d\vec{x}= 0$$

siccome le funzioni sono lisce e la precedente equazione vale su ogni dominio $\Omega\subset \mathbb{R}$ abbiamo che l' integranda è nulla, di conseguenza otteniamo l'equazione alle derivate parziali nota con il nome di equzione di continuità:

$$\partial_t \rho + \nabla \cdot \vec{J} = 0$$

Studiamo ora dei metodi che ci permettono di risolvere in generale equazioni scritte in questa forma, i metodi ai volumi finiti (FVM) Partiamo dall'equaione di continuità in una dimensione e consideriamo come dominio l'intervallo chiuso $[a,b]$. Allora dividiamo l'intervallo in dei volumi finiti, ovvero dei sotto intervalli $V_i=(x_{i-1},x_i)$ dove:

$$x_0 = a\qquad x_n = b \qquad x_i = i\cdot h_x$$

a questo punto siamo pronti a definire la quantità $\overline{\rho}_i(t) = \frac{1}{h_x}\int_{V_i} \rho(t,x) dx$, che prende il nome di densità media di volume sull' i-esimo volume finito. A questo punto integriamo l'equazione di continuità per ottenere la seguente euqazione:

$$\int^{t_{j+1}}_{t_j} \partial_t \rho dt = -\int_{t_j}^{t_j{+1}}\partial_x J dt \Rightarrow \rho(t_{j+1},\vec{x})=\rho(t_j,\vec{x}) - \int_{t_j}^{t_{j+1}} \partial_x J dt$$

dove $t_j = j*h_t$. Integriamo ora in spazio l'equazione precedente e dividiamo per $h_x$ in modo tale da ottenere un equazione in termini delle densità medie di volume:

$$\overline{\rho}_i (t_{j+1}) = \overline{\rho}(t_j) - \int_{x_{i-1}}^{x_i} dx\frac{1}{h_x}\int_{t_j}^{t_{j+1}} \partial_x J dt$$

Siccome stiamo lavorando con funzioni lisce possiamo invertire l'ordine di integrazione (per il teorema di Fubini-Tonelli):

$$\overline{\rho}_i(t_{j+1}) = \overline{\rho}_i(t_{j}) - \int_{t_j}^{t_{j+1}} dt \frac{J_{i}(t) - J_{i-1}(t)}{h_x}$$$$\overline{\rho}_i(t_{j+1}) = \overline{\rho}_i(t_{j}) - \frac{1}{h_x}\int_{t_j}^{t_{j+1}}J_{i}(t) - J_{i-1}(t)$$

portando al limite $h_t\to 0$ (e utilizzando il teorema della media integrale che vedremo a breve) otteniamo la seguente discrettizzazione in spazio dell' equazione di continuità:

$$\partial_t \overline{\rho}_i = \frac{1}{h_x} \Big(J_i(t)-J_{i-1}(t)\Big)$$

Questa equazione è uno schema numerico per l'approsimazione di continuità, in una dimensione questa approssimazione corrisponde a un approssimazione del primo ordine utilizzando il metodo delle differenze finite, invece in più dimensioni i due metodi si distinguono. Osserviamo ora che abbiamo trovato uno schema numerico per approssimare la soluzione dell'equazione di continuità, poichè per un $h_x$ sufficentemente piccolo $\overline{\rho}_{i+1}(t)\approx \rho_{i+\frac{1}{2}}(t)$.

Teorema (Media Integrale) Data una funzione $f:[a,b]\to \mathbb{R}$ continua, esiste $c\in (a,b)$ tale che:

$$\frac{1}{b-a}\int_a^bf(x) = f(c)$$

Corollario Fissato un punto $x_0\in (a,b)$ abbiamo che:

$$\rho(t,x_0) = \frac{1}{h_x}\int_{x_0-\frac{h_x}{2}}^{x_0 + \frac{h_x}{2}} \rho(t,x) dx = \overline{\rho}(t,x_0)$$

Esempio (Equazione Del Calore) Consideriamo ora come legge fisica la legge di Fourier, che ci dice che le particelle si muvono contro gradiente:

$$ \vec{J}(t,\vec{x}) = -\nabla \rho (t, \vec{x})=-\partial_x \rho(t,x)$$

dunque dobbiamo approssimare la quantità, $J_i(t)-J_{i-1}(t)$, per fare questo utilizziamo le differenze finite del primo ordine in avanti e all'indietro:

$$J_i(t) = k\bigg[\frac{\rho_{i+1}-\rho_{i}}{h_x}\bigg] \qquad J_{i-1}(t) = k\bigg[\frac{\rho_{i}-\rho_{i-1}}{h_x}\bigg]$$

$$J_i(t) - J_{i-1}(t) = \frac{k}{h_x}\bigg(\rho_{i+1}-2 \rho_i + \rho_{i-1}\bigg)$$

dunque l'equazione che rappresenta il metodo ai volumi finiti diventa:

$$ \partial_t \overline{\rho} = \frac{1}{h_x^2}\bigg(\rho_{i+1}-2\rho_i+\rho_{i-1}\bigg)$$

in virtu del precedente corollario possiamo sostuire $\rho_i$ con $\overline{\rho}_{i}$, per ottenere:

$$ \partial_t \overline{\rho} = \frac{1}{h_x^2}\bigg(\overline{\rho}_{i+1}-2\overline{\rho}_i+\overline{\rho}_{i-1}\bigg)$$

integrando questa semi-discretizzazione sapaziale con il metodo di Eulero otteniamo il seguente schema numerico:

$$ \overline{\rho}(t_{j+1}) = \overline{\rho}(t_{j})+\frac{h_t}{h_x^2}\bigg(\overline{\rho}_{i+1}(t_j)-2\overline{\rho}_i(t_j)+\overline{\rho}_{i-1}(t_j)\bigg)$$

per chi è famigliare con il metodo delle differenze finite ques'ultimo schema numerico corrisponde a un approssimazione alle differenze finite del secondo ordine.

In [2]:
%matplotlib inline
from matplotlib import animation, rc
from IPython.display import HTML
In [3]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0
    
#Inizializziamo l' asse x
a= 0; b = 4; T=0.3;
hx = 0.02; ht = 0.0001
print("Curant Number: {}".format(ht/hx**2))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.5)*hx)
plt.plot(t,Rho)
Rho[0]=0.0; Rho[-1]=0.0;
for i in range(1,M+1):
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]+(ht/(hx**2))*(Rho[i+1]-2*Rho[i]+Rho[i-1]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
plt.plot(t,RhoN)
plt.title("Heat at t:"+str(M*ht))

Curant Number: 0.25
Out[3]:
Text(0.5, 1.0, 'Heat at t:0.2999')
In [4]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0

#Inizializziamo l' asse x
a= 0; b = 4; T=3;
hx = 0.05; ht = 0.001
print("Curant Number: {}".format(ht/hx**2))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.05)*hx)
Rho[0]=0.0; Rho[-1]=0.0; #Impongo condizioni al bordo
#Costruiamo il plot
fig, ax = plt.subplots()

ax.set_xlim(( 0, 4))
ax.set_ylim((0, 1.1))

line, = ax.plot([], [], lw=2)
plt.title("Heat Equation")
def init():
    line.set_data([], [])
    return (line,)
def animate(j):
    global Rho;
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]+(ht/(hx**2))*(Rho[i+1]-2*Rho[i]+Rho[i-1]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
    line.set_data(t, Rho)
    return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=1000, interval=5, blit=True);
HTML(anim.to_html5_video())

Curant Number: 0.3999999999999999
Out[4]:

Osservazione (Numero Di Courant)

In [7]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0

#Inizializziamo l' asse x
a= 0; b = 4; T=3;
hx = 0.05; ht = 0.00252
print("Curant Number: {}".format(ht/hx**2))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.05)*hx)
Rho[0]=0.0; Rho[-1]=0.0; #Impongo condizioni al bordo
#Costruiamo il plot
fig, ax = plt.subplots()

ax.set_xlim(( 0, 4))
ax.set_ylim((0, 1.1))

line, = ax.plot([], [], lw=2)
plt.title("Heat Equation")
def init():
    line.set_data([], [])
    return (line,)
def animate(j):
    global Rho;
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]+(ht/(hx**2))*(Rho[i+1]-2*Rho[i]+Rho[i-1]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
    line.set_data(t, Rho)
    return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=1000, interval=5, blit=True);
HTML(anim.to_html5_video())

Curant Number: 1.0079999999999998

/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:36: RuntimeWarning: overflow encountered in double_scalars
/usr/local/lib/python3.7/site-packages/ipykernel_launcher.py:36: RuntimeWarning: invalid value encountered in double_scalars
Out[7]:

Osserviamo che perchè il metodo numerico sia stabile dobbiamo richiedere che il numero di Courant rispetti la seguente disuguaglianza:

$$ \frac{1}{k}\frac{h_t}{h_x^2} \leq 1$$

Esempio (Equazione Dei Mezzi Porosi) Prendiamo ora come legge fisica la legge di Darcy, che ci dice che la velocità delle particelle nei mezzi porosi è data dal equazione:

$$ \vec{v} = -\frac{k}{\mu} \nabla \rho = -\frac{k}{\mu}\partial_x \rho$$

siccome $\vec{J}=\rho \vec{v}$ possiamo scrivere l'equazione:

$$J = \frac{k}{\mu} \rho \partial_x \rho$$

sostituendo quest'ultimo risultato all'interno dell' equazione di continuità otteniamo:

$$\partial_t \rho = -\nabla\cdot (\rho v) = - \partial_x \Bigg( \rho \Big( -\frac{k}{\mu} \partial_x \rho \Big) \Bigg)$$

Dunque anche in questo caso dobbiamo trovare un modo per approssimare il flusso e faremo queto utilizzando le differenze finite in avanti e all' indietro:

$$J_i(t) = \frac{k}{\mu} \rho_i(t) \bigg( \frac{\rho_{i+1}(t) -\rho_i(t)}{h_t}\bigg)\qquad J_{i-1}(t) = \frac{k}{\mu} \rho_i(t) \bigg( \frac{\rho_{i}(t) -\rho_{i-1}(t)}{h_t}\bigg)$$$$J_i(t)-J_{i-1}(t) = -\frac{k}{\mu}\rho_i(t)\bigg(\frac{\rho_i(t)-\rho_{i-1}(t)}{h_t}\bigg)$$

questo mi porta alla seguente semi discretizzazione numerica per il metodo l'equazione dei mezzi porosi:

$$\partial_t \overline{\rho}_i(t) = -\frac{k}{\mu}\overline{\rho}_i(t)\bigg(\frac{\overline{\rho}_i(t)-\overline{\rho}_{i-1}(t)}{h_t}\bigg)$$
In [6]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0

#Inizializziamo l' asse x
a= 0; b = 4; T=3;
hx = 0.05; ht = 0.001
print("Curant Number: {}".format(ht/hx**2))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.05)*hx)
Rho[0]=0.0; Rho[-1]=0.0; #Impongo condizioni al bordo
#Costruiamo il plot
fig, ax = plt.subplots()

ax.set_xlim(( 0, 4))
ax.set_ylim((0, 1.1))

line, = ax.plot([], [], lw=2)
plt.title("Porous Media Equation")
def init():
    line.set_data([], [])
    return (line,)
def animate(j):
    global Rho;
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]+(ht/(hx**2))*Rho[i]*(Rho[i+1]-2*Rho[i]+Rho[i-1]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
    line.set_data(t, Rho)
    return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=1000, interval=5, blit=True);
HTML(anim.to_html5_video())

Curant Number: 0.3999999999999999
Out[6]:

Esercizio (Equazione Del Trasporto e Metodi UpWind) Consideriamo ora come legge fisica per il flusso, la legg del trasporto:

$$\vec{J}(t,\vec{x}) = v\rho(t,x)$$

l'equazione che si ottiene dall 'equazione di continuità in questo caso è:

$$\partial_t \rho = v\partial_x \rho$$

si osservi cosa succede allo schema numerico quando abbiamo una velocità postiva e una velocità negativa.

In [4]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0

#Inizializziamo l' asse x
a= 0; b = 4; T=3;
hx = 0.05; ht = 0.001
print("Curant Number: {}".format(ht/hx))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.05)*hx)
Rho[0]=0.0; Rho[-1]=0.0; #Impongo condizioni al bordo
#Costruiamo il plot
fig, ax = plt.subplots()

ax.set_xlim(( 0, 4))
ax.set_ylim((0, 1.1))

line, = ax.plot([], [], lw=2)
plt.title("Transport Equation")
def init():
    line.set_data([], [])
    return (line,)
def animate(j):
    global Rho;
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]-(ht/(hx))*(Rho[i+1]-Rho[i]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
    line.set_data(t, Rho)
    return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=1000, interval=5, blit=True);
HTML(anim.to_html5_video())

Curant Number: 0.02
Out[4]:
In [5]:
import numpy as np
import matplotlib.pyplot as plt
def g(x):
    if  1 <= x <= 2:
        return 1
    else:
        return 0

#Inizializziamo l' asse x
a= 0; b = 4; T=3;
hx = 0.05; ht = 0.001
print("Curant Number: {}".format(ht/hx))
N = int((b-a)/hx);
M = int(T/ht);
t = np.zeros(N+1);
Rho = np.zeros(N+1)
for i in range(0,N+1):
    t[i] = (i+0.5)*hx
    Rho[i] = g((i+.05)*hx)
Rho[0]=0.0; Rho[-1]=0.0; #Impongo condizioni al bordo
#Costruiamo il plot
fig, ax = plt.subplots()

ax.set_xlim(( 0, 4))
ax.set_ylim((0, 1.1))

line, = ax.plot([], [], lw=2)
plt.title("Transport Equation")
def init():
    line.set_data([], [])
    return (line,)
def animate(j):
    global Rho;
    RhoN = np.zeros(N+1)
    for i in range(1,N+1-1):
        RhoN[i] = Rho[i]+(ht/(hx))*(Rho[i+1]-Rho[i]);
    RhoN[0]=0.0; RhoN[-1]=0.0;
    Rho = RhoN
    line.set_data(t, Rho)
    return (line,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=1000, interval=5, blit=True);
HTML(anim.to_html5_video())

Curant Number: 0.02
Out[5]: