Recordemos que Navier-Stokes no es mas que Euler donde hemos añadido viscosidad \mu, o equivalentemente pero al revés, Euler no es mas que Navier-Stokes donde hemos eliminado la viscosidad (\mu = 0). La viscosidad, de hecho, solo modifica las ecuaciones del momento.

¿Y como quedan estas en el caso de fluidos incompresibles (asumimos, como hasta ahora, que \rho = 1) cuando se encuentran en estado estacionario? Pues son las siguientes (expandidas, que es como las necesitamos al tratarlas numéricamente). Para el momento, nos queda:

\frac{\partial}{\partial x} (v_x^2 + p) + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z = \mu \Delta v_x

\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} (v_y^2 + p) + \frac{\partial}{\partial z} v_y v_z = \mu \Delta v_y

\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} (v_z^2+p) = \mu \Delta v_z

Y la de continuidad queda:

\frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z = 0

Una forma de proceder sería la indicada en Solución numérica a las ecuaciones de Euler incompresible n-D (estado estacionario).

Sin embargo, el proceder con Navier-Stokes era para tener operadores elípticos. Podríamos pensar en resolver numéricamente el sistema

\begin{cases} \Delta v_x = \frac{1}{\mu} [ \frac{\partial}{\partial x} (v_x^2 + p) + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z ] \\ \Delta v_y = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} (v_y^2 + p) + \frac{\partial}{\partial z} v_y v_z] \\ \Delta v_z = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} (v_z^2+p) ] \\ 0 = \frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z \end{cases}

pero existen una serie de problemas.

Para empezar, la última ecuación sigue siendo con derivadas primeras. Además, las tres ecuaciones que aparentemente nos deberían llevar a encontrar v=(v_x, v_y, v_z) presentan, en sus respectivos términos fuente, la también desconocida p, de manera que, tal y como esta escrito, no podemos abordar el problema numéricamente (¿qué valores ponemos en p). Si en la ecuación de continuidad apareciera p, podríamos escribir ésta en función del campo de velocidades y sustituir en las ecuaciones del momento, pero tampoco es el caso, pues en esta última solo se explicita una condición que debe de cumplir v. Por tanto, la p también tiene que aparecer en la parte izquierda de las ecuaciones

\begin{cases} \Delta v_x - \frac{\partial}{\partial x} p  = \frac{1}{\mu} [ \frac{\partial}{\partial x} v_x^2 + \frac{\partial}{\partial y} v_x v_y + \frac{\partial}{\partial z} v_x v_z ] \\ \Delta v_y - \frac{\partial}{\partial y} p= \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_y + \frac{\partial}{\partial y} v_y^2 + \frac{\partial}{\partial z} v_y v_z] \\ \Delta v_z - \frac{\partial}{\partial z} p = \frac{1}{\mu} [\frac{\partial}{\partial x} v_x v_z + \frac{\partial}{\partial y} v_y v_z + \frac{\partial}{\partial z} v_z^2 ] \\ 0 = \frac{\partial}{\partial x} v_x + \frac{\partial}{\partial y} v_y + \frac{\partial}{\partial z} v_z \end{cases}

La discretización del término de la izquierda de la primera ecuación quedaría:

\frac{v_{x_{i-1,j,k}}-2v_{x_{i,j,k}}+v_{x_{i+1,j,k}}}{h_x^2} + \frac{v_{x_{i,j-1,k}}-2v_{x_{i,j,k}}+v_{x_{i,j+1,k}}}{h_y^2}

+ \frac{v_{x_{i,j,k-1}}-2v_{x_{i,j,k}}+v_{x_{i,j,k+1}}}{h_z^2} - \frac{p_{i+1,j,k}-p_{i-1,j,k}}{2h_x}

Puesto que en el esquema iterativo solo puede quedar una incógnita en la parte izquierda de la iteración ¿Estamos consiguiendo algo? Parece que no… Ya que, como hemos dicho, una de las ecuaciones solo incluye velocidades, tendríamos que reconvertir alguna de las del Laplaciano en una iteración para la p. Sin embargo, junto con el orden de la discretización para incluir el punto actual,  el problema que tenemos es que la presión aparece con primera derivadas…

Anuncios