Test 2: Caso 1D - Couette Generalizado

En el flujo de Couette generalizado, el fluido se encuentra confinado entre placas planas paralelas infinitas separadas por una distancia constante $L$ y se encuentra sometido a un gradiente de presión en dirección paralela a las placas. Adicionalmente, las velocidades de las placas son de $U_0$ y $U_1$ (en dirección x) respectivamente. De acuerdo con esto, se tienen las siguientes consideraciones:

• El flujo se presenta únicamente en la dirección $x$ ($v=w=0$).
• El gradiente de presión al que está sometido el fluido es constante ($\frac{\partial p}{\partial x} = \frac{\Delta p}{\ell}$).
• El flujo no depende de la coordenada $z$ ($\frac{\partial}{\partial z}=0$).
• El fluido se encuentra en estado estacionario ($\frac{\partial}{\partial t}=0$).
• No hay fuerzas externas afectando el flujo ($f_x=f_y=f_z=0$).

De aquí, las ecuaciones de continuidad y de Navier-Stokes se reducen a:

\begin{align} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho \vec{u})=0 \quad &\to\quad \frac{\partial u}{\partial x}=0 \\ \rho\left(\frac{\partial \vec{u}}{\partial t}+\vec{u}\cdot\nabla\vec{u}\right)=\rho\vec{f}-\nabla p+\mu\nabla^2\vec{u}\quad &\to\quad 0=-\frac{\Delta p}{\ell}+\mu\frac{\partial^2u}{\partial y^2} \end{align}

Se puede ver entonces que la ecuación resultante es un caso particular del problema de difusión estable:

$$k\frac{d^2 f}{d x^2}+s(x)=0$$

Donde $k = \mu$, y $s(x)=-\frac{\Delta p}{\ell}$, y las condiciones de frontera corresponden a valores conocidos (tipo Dirichlet):

$$u(0)=U_0\quad u(L)=U_1$$

Ahora, para aplicar el método de elementos espectrales, es necesario realizar los siguientes pasos:

• Introducción de la información del problema
• Definición de los elementos a usar (cantidad $N_E$ y distribución, así como sus órdenes $N_P(i)$)
  • Selección de los nodos de interpolación a usar
  • Generación de la matriz de conectividad
  • Definición de la función de entrada
• Generación de la matriz del problema
  • Generación de la matriz de difusión y el vector b
• Solución del sistema resultante
• Gráfica de la solución obtenida

Cada uno de los pasos se tratará en detalle a continuación

Introducción de la información del problema

Se introducen los parámetros del problema. En este caso se tienen los parámetros $\mu$, $\nabla p$, $L$ y $U_0$.

Definición de los elementos a usar

Se define primero el número de elementos y su distribución, así como los órdenes de cada uno de los mismos. Por simplicidad, inicialmente se toman elementos uniformemente espaciados.

Selección de los nodos de interpolación a usar

Conociendo el número de elementos y su distribución, es posible generar los nodos correspondientes a los extremos de cada elemento. Sin embargo, como no se trabaja con elementos lineales (cada elemento tiene un orden), es necesario generar los nodos intermedios para cada elemento. Para ello, es necesario definir la familia de polinomios sobre la cual van a estar basados estos nodos. En este caso se utilizan los nodos de interpolación de Lobatto.

Generación de la matriz de conectividad

Una vez se definen los nodos, se utiliza la matriz de conectividad para almacenar de manera ordenada la numeración local y global de los mismos

Definición de la función de entrada

Especifica $s(x)$ para cada uno de los nodos

Generación de la matriz del problema

Para la generación de la matriz del problema es necesario definir cada una de las matrices de difusión y de masa por cada elemento:

$$\Theta_{ij}=\int_{-1}^1 \psi_i(\xi)\psi_j(\xi)d\xi \quad \Psi_{ij}=\int_{-1}^1 \frac{d\psi_i}{d\xi}\frac{d\psi_j}{d\xi}d\xi$$

En este caso se aplica la cuadratura de Lobatto para la construcción de dichas matrices, con lo que se obtienen los siguientes resultados:

$$\Theta=\sum_{p=1}^{N_p+1} m_{ip}m_{jp}w_p \quad \Psi_{ij}=\sum_{p=1}^{N_p+1} d_{ip}d_{jp}w_p$$

Generación de la matriz de difusión y el vector b

Se recorre sobre los elementos para acoplar las matrices de cada elemento en un sistema global, aprovechando la matriz de conectividad

Solución del sistema resultante

Teniendo ya construidos los componentes del sistema a resolver, basta aplicar un esquema de solución razonable para el tipo de matriz obtenida. En este caso se utiliza el solucionador por defecto del módulo linalg. Dado que los valores para el primer y último nodo son conocidos, se omiten los valores correspondientes.

Gráfica de la solución obtenida

Una vez resuelto el sistema, es posible visualizar el resultado obtenido por medio de una gráfica. En este caso se conoce que la solución analítica corresponde a la linea recta entre las dos condiciones extremas, y se incluye también su gráfica.