Test 1: Caso 1D - Couette

En el flujo de Couette, el fluido se encuentra confinado entre placas planas paralelas infinitas separadas por una distancia constante $L$. Adicionalmente, la velocidad de las placas es 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 flujo no depende de la coordenada $z$ ($\frac{\partial}{\partial z}=0$).
El fluido se encuentra en estado estacionario ($\frac{\partial}{\partial t}=0$).
La presión a la que está sometido el fluido es constante ($\nabla p = 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{\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 = 1$, y $s(x)=0$, 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
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



In [1]:

    
# Prepara el ambiente de trabajo con las librerías a utilizar

import numpy as np               # Ofrece funciones para operaciones básicas con arreglos
import scipy.linalg as lin       # Ofrece funciones para operaciones de álgebra lineal
import matplotlib.pyplot as plt  # Permite incluir un ambiente de visualización
import timeit                    # Permite la medición de tiempos para los distintos algoritmos
from sem1D import *               # Agrupa las funciones externas en un archivo

Introducción de la información del problema

Se introducen los parámetros del problema. En este caso se tienen únicamente los parámetros $L$ y $U_0$.



In [2]:

    
L = 2.0
U0 = 15.0
U1 = 25.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.



In [3]:

    
Ne = 5
Np = np.array((6,4,5,3,6))

# Genera los puntos correspondientes a los nodos extremos de cada elemento
xe = np.linspace(0,L,Ne+1)

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.

def genNodes(Ne,Np,xe): #================================================= # Genera los nodos de interpolacion segun el # numero de elementos y sus ordenes. Esta es la # numeracion local de los nodos: (l,i) # # Ne = Np.size # max(Np) < 6 #================================================= if(Ne != Np.size): raise NameError('La cantidad de elementos no coincide con el orden') xint = np.zeros((Ne,max(Np)+1)) for l in range(Ne): mid = 0.5 * (xe[l+1]+xe[l]) dis = 0.5 * (xe[l+1]-xe[l]) xint[l,0] = xe[l] vals = lobatto(Np[l]-1) xint[l,1:Np[l]] = mid + vals*dis xint[l,Np[l]] = xe[l+1] return xint



In [4]:

    
xint = genNodes(Ne,Np,xe)

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

def genConn(Ne,Np,xint): #================================================= # Genera los vectores de numeracion global y la # matriz de conectividad # # max(Np) < 6 #================================================= xglob = np.zeros(np.sum(Np)+1) C = np.zeros((Ne,max(Np)+1),int) # Contador Global cont = 0 for l in range(Ne): C[l,0:Np[l]+1] = np.arange(cont,cont+Np[l]+1) xglob[cont:cont+Np[l]+1] = xint[l,0:Np[l]+1] cont+=Np[l] return xglob, C



In [5]:

    
(xglob, C) = genConn(Ne,Np,xint)



In [6]:

    
# Permite visualizar todos los nodos a utilizar

plt.plot(xe,np.zeros(xe.shape),'bo',xglob,np.zeros(xglob.shape),'r.')
plt.xlabel('Posicion (coord y)')
plt.show()

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

def matDiff(Ne,Np,xe,xglob,C,U0,U1,mu=1,s=lambda x: 0): #================================================= # Genera la matriz de difusion/rigidez y el vector # del lado derecho del sistema sujeto a condicio- # nes de Dirichlet en ambos extremos # # max(Np) < 6 #================================================= Ng = xglob.size Mat_dif = np.zeros((Ng,Ng)) vec_b = np.zeros(Ng) s = np.vectorize(s) for l in range(Ne): N = Np[l] h = xe[l+1]-xe[l] elm_mm = 0.5*h*emm(N) elm_dm = 2.0*edm(N)/h Mat_dif[C[l,0]:C[l,N]+1,C[l,0]:C[l,N]+1] += elm_dm vec_b[C[l,0:N+1]] += np.dot(elm_mm,s(xglob[C[l,0:N+1]]))/mu # Implementa las condiciones de frontera de Dirichlet # Extremo derecho # los valores de N, h y elm_dm ya corresponden a los del ultimo nodo vec_b[C[Ne-1,0:N+1]] -= elm_dm[0:N+1,N]*U1 # Extremo izquierdo N = Np[0] h = xe[1]-xe[0] elm_dm = 2.0*edm(Np[0])/h vec_b[C[0,0:N+1]] -= elm_dm[0:N+1,0]*U0 return Mat_dif, vec_b



In [7]:

    
# Permite visualizar la manera en que se acoplan los bloques

(Mat_dif, vec_b) = matDiff(Ne,Np,xe,xglob,C,U0,U1)

plt.spy(Mat_dif)
plt.show()

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.



In [8]:

    
sol = np.zeros(xglob.size-2)
A = Mat_dif[1:-1,1:-1]
sol = lin.solve(A,vec_b[1:-1])
solU = np.concatenate([[U0],sol,[U1]])

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.



In [9]:

    
plt.plot(xe,np.linspace(U0,U1,Ne+1),'-.or',ms=8)
plt.hold(True)
plt.plot(xglob,solU,':*b')
plt.xlabel('Posicion (coord y)')
plt.ylabel('Velocidad (U)')
plt.show()



In [ ]: