Difusión en 1D

Igual que con la convección, podemos intentar discretizar la ecuación:

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$$

o podemos remontarnos a la deducción de la misma

Deducción

Se suponen cambios de un campo $u$ sólo en la dirección $x$

El cambio en la cantidad total $A \Delta x \, u_i$ será:

$$\frac{d }{d t} (A \Delta x \, u_i ) = \Phi_{i-1/2} - \Phi_{i+1/2}$$

Flujos, convección

Antes los flujos por las caras venían dados por:

$$\Phi_{i-1/2} = A c \, u_{i-1/2}$$

$$\Phi_{i+1/2} = A c \, u_{i+1/2}$$

(Recordemos que había que aproximar $u_{i-1/2}$ y $u_{i+1/2}$, que son desconocidas, y vimos con python el "centered diff", el "upwind" y el "downwind")

Flujos, difusión de Fick

Ahora los flujos vienen de una ley física, la ley de Fick:

$$\Phi = - D \frac{\partial u}{\partial x} .$$

Lo que expresa es que la difusión tiende a "allanar" los campos, creando un flujo desde zonas con valores altos a otras con valores bajos.

(En realidad, Fick es para concentración; para temperatura se llama "ley de Fourier".)

Flujos, difusión

Así pues, parece natural

$$\Phi_{i-1/2} \approx - A D \frac{ u_i - u_{i-1/2}}{\Delta x}$$

$$\Phi_{i+1/2} \approx - A D \frac{ u_{i+1} - u_i }{\Delta x}$$

Los cuales, enchufados en

$$\frac{d }{d t} (A \Delta x \, u_i ) = \Phi_{i-1/2} - \Phi_{i+1/2} \ldots ...$$

Ecuación de difusión

... LLevan a:

$$\frac{d u_i }{d t} = D \frac{u_{i-1} - 2 u_i + u_{i+1} }{\Delta x^2}$$

Lo cual, por supuesto, lleva en el continuo a la ecuación tan conocida:

$$\frac{\partial u}{\partial t} = D \frac{\partial^2 u}{\partial x^2}$$

Algoritmo

Ya sólo queda discretizar el tiempo. Por ejemplo, explícitamente:

$$u_i^{n+1} = u_i^n + D \frac{\Delta t}{\Delta x^2 } (u^n_{i-1} - 2 u^n_i + u^n_{i+1} )$$

O, implícitamente:

$$u_i^{n+1} = u_i^n + D \frac{\Delta t}{\Delta x^2 } (u^{n+1}_{i-1} - 2 u^{n+1}_i + u^{n+1}_{i+1} )$$

Algoritmo final

... O incluso métodos mixtos, como el famoso de Crank-Nicolson $$\bar{u}_i = (u^{n+1}_i+u^n_i)/2$$

$$u_i^{n+1} = u_i^n - D \frac{\Delta t}{\Delta x^2 } (\bar{u}_{i-1} - 2 \bar{u}_i + \bar{u}_{i+1} )$$

Además, siempre aparece un nuevo "número de Courant difusivo" (en realidad, "número difusivo"):

$$\mathrm{Co} := D \frac{\Delta t}{\Delta x^2 }$$

Hala, a implementar todo esto

Referencias

12 Steps to Navier-Stokes