Volúmenes Finitos y la ecuación de transporte

Métodos Numéricos Avanzados. 2023-2024

David Gómez-Castro

Universidad Complutense de Madrid

Volúmenes Finitos

Ecuación de continuidad

Trabajamos en dimensión \(d = 1\) y estudiamos leyes de conservación de la forma \[\begin{equation} \tag{C} \label{eq:conservacion} \frac{\partial}{\partial t} \rho = - \diver \flux \end{equation}\]

integrando en \([x_{i-\frac 1 2}, x_{i+\frac 1 2}]\) tomamos la semi-discretización espacial \[\begin{equation} \frac{\diff}{\diff t} \rho_i (t) = - \frac{ \flux_{i+\frac 1 2} (t) - \flux_{i-\frac 1 2} (t) }{x_{i+\frac 1 2} - x_{i-\frac 1 2}}, \qquad i = 1, \cdots, N. \end{equation}\]

En nuestro tipo de leyes de conservación, será frecuente que estudiemos EDPs cuyas soluciones exactas sean positivas: \(\rho \ge 0\).

Muchos de los razonamientos se pueden generalizar al caso de cambio de signo.

Conservación

Estos métodos son naturales para leyes de conservación porque \[\begin{align*} \frac{\diff}{\diff t} \sum_{i=1}^N \Big(x_{i+\frac 1 2} - x_{i-\frac 1 2}\Big)\rho_i (t) &= - \left( \flux_{\frac 3 2} (t) - \flux_{\frac 1 2} (t) \right) - \left( \flux_{\frac 5 2} (t) - \flux_{\frac 3 2} (t) \right) - \cdots - \left( \flux_{N-\frac 1 2} (t) - \flux_{N + \frac 1 2} (t) \right) \\ &= - F_{N+\frac 1 2} (t) + F_{\frac 1 2} (t) \end{align*}\]

Conservación de masa es equivalente a que \[ F_{N+\frac 1 2} (t) = F_{\frac 1 2} (t) \]

En algunas ocasiones se toman ambos flujos igual a 0 para que el problema esté bien determinado.

En otras ocasiones tomaremos condiciones periódicas en \(\rho\), es decir \[ \rho_0 = \rho_N, \qquad \rho_1 = \rho_{N+1}. \]

Ecuación de Transporte

Ecuación de transporte

El caso más sencillo es \(\flux(t,x) = a \rho (t,x) \ge 0\) con \(a\) constante

\[\frac{\partial \rho}{\partial t} + a\frac{\partial \rho}{\partial x} = 0\]

Soluciones exactas

Suelen llamarse características a curvas \((t,X(t))\) tales que la EDP se puede resolver aisladamente a lo largo de curva.

En este caso, buscamos curvas \(X(t)\) tales que \(\rho(t, X(t)) = \rho_0 (y)\), donde \(a(0) = 1\). \[ \partial_t \rho + \frac{\partial \rho}{\partial x} \frac{\diff X}{\diff t} = 0. \]

Luego \(\frac{\diff X}{\diff t}(t) = a\). Y \(X(t,y) = y + at\).

Finalmente resolviendo \(x = y + at\) tenemos \[ \rho(t,x) = \rho_0 (x - a t). \] Es decir, solución en el punto \((t,x)\) depende sólo del valor inicial en \((0,x - at)\).

Restricción a un intervalo cerrado

Si queremos trabajar un intervalo acotado, digamos \((0,1)\), no puede admitir condiciones de flujo nulo.

Puede, sin embargo, admitir condiciones de periódicas. Por ejemplo, para \(t \in (0,1)\) y \(a > 0\) \[ \rho(t,0) = \rho(t,1) = \rho_0(1 - at). \]

Imponer condiciones de flujo periódico consiste en decir \[ a \rho(t,1) = F(t,1) = F(t,0) = a \rho(t,0). \]

Desde el punto de vista exacto, es de nuevo imponer condiciones de periódicas en \(\rho\).

Discretización temporal

Hemos de elegir un método de discretizar \(\frac{\diff \rho_i}{\diff t}\).

Probaremos el método de Euler \[ \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = - \frac{ F_{i+\frac 1 2}^{n+1} - F_{i - \frac 1 2}^{n+1} }{\Delta x} \]

Para que sea explícito \(F_{i+\frac 1 2}^{n+1} = G_{i+\frac 1 2} ( \rho_i^n, \rho_{i+1}^n )\).

Es decir \(F_{i+\frac 1 2}^{n+1} = a\rho_i^n , a\rho_{i+1}^n\) u otro valor.

Experimentación numérica

Hagamos unos experimentos para ver qué elección de \(\rho_i\) funciona mejor.

using Plots

N  = 64
T  = 1
x  = range(-1.0,1.0,N+1); x = vec(x[1:end-1])
Δx = x[2] - x[1] 
ρ0 = min.(1,max.(-4*abs.(x).+2,0));
a  = 1;

\(F_{i+\frac 1 2}^{n+1} = a\rho_i^n\)

Tenemos el método \[ \rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i}^n - \rho_{i-1}^n ) \]

ρ = copy(ρ0); Δt = Δx/2
@gif for n=1:ceil(T/Δt)
    ρn = copy(ρ)
    
    ρ[1] = ρn[1] - Δt/Δx * a * (ρn[1] - ρn[end])
    for i = 2:length(ρ)
        ρ[i] = ρn[i] - a * Δt/Δx * ( ρn[i] - ρn[i-1] )
    end
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",
        linewidth=2,thickness_scaling = 1.5,label="") 
end

Ejercicio. \(F_{i+\frac 1 2}^{n+1} = a\rho_{i+1}^n\)

\[\rho^{n+1}_i = \rho^n_i - a\frac{\Delta t}{\Delta x} ( \rho_{i+1}^n - \rho_{i}^n )\]

Code

ρ = copy(ρ0); Δt = Δx/2
@gif for n=1:ceil(T/Δt)
    ρn = copy(ρ)
    
    ρ[1:end-1] = ρn[1:end-1] - Δt/Δx * a * ( ρn[2:end] - ρn[1:end-1] )
    ρ[end] = ρn[end] - a * Δt/Δx * ( ρn[1] - ρn[end] )
    
    plot(x,ρ,
        ylim=(0.0,1.2),
        title="t=$(round(n*Δt,digits=3))",
        linewidth=2,thickness_scaling = 1.5,label="") 
end