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;Volúmenes Finitos y la ecuación de transporte
Métodos Numéricos Avanzados. 2022-2023
Volúmenes Finitos
Leyes de conservación
\[ \newcommand{\diff}{\,\mathrm{d}} \DeclareMathOperator{\diver}{div} \newcommand{\flux}{{F}} \]
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. \]
Ecuación de Transporte
Ecuación de transporte
El caso más sencillo es \(\flux = a > 0\)
\[\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.
\(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=$(n*Δt)",
linewidth=2,thickness_scaling = 1.5,label="")
end