Proyecto final: Método de Elementos Finitos.

Aldo Sayeg Pasos Trejo

Introducción a los Medios Continuos. Prof. Roberto #Romero Arias

Posgrado en Ciencias Matemáticas. Universidad Nacional Autónoma de México

Introducción

A pesar de que los métodos de diferencias finitas nos permiten encontrar soluciones numéricas con una buena aproximación para la gran mayoría de los problemas, a la hora de plantear un problema que requiera de una geometría más exótica (es decir, cualquier cosa que no sea un cilindro, una esfera, un rectángulo o una región que no se pueda transformar continuamente en las anteriores), el problema se vuelve muy complicado de resolver.

La importancia de problemas con geometría exótica es muy grande y, por lo general, surgen frecuentemente en problemas enfocados en sistemas reales, normalmente estudiados por distintas ramas de las ingenierías. Para resolver dichos problemas, una opción es el método de elementos finitos (“Finite Element Method”, FEM), utilizado muy frecuentemente en aplicaciones que requieren mucha precisión.

A cambio de poder tratar geometrías más ricas, FEM requiere de un trasfondo matemático más fuerte pues requiere, en primer lugar, teselar nuestro dominio en varias subregiones poligonales para que podamos replantear la ecuación diferencial ordinaria o parcial a resolver como un problema de cálculo de variaciones, cuya solución estará dada por un conjunto de ecuaciones algebraicas sobre los vértices de la teselación.

Con la descripción anterior, inmeadiatamente surgen dos dudas sobre los FEM

1. ¿Cómo debo teselar mi dominio para que sea óptima mi solución?
2. ¿Las ecuaciones algebráicas que surgan, se podrán resolver de manera sencilla?

En el presente trabajo, utilizaremos los métodos de elemento finito para resolver un problema sencillo de dinámica de fluidos: un tunel de viento en el que un fluido incompresile se encuentra con un obstáculo rígido e inamovible.

Marco Teórico

Triangulaciones

Antes de proceder a lo que propiamente se denomina elemento finito, procedemos a dar explicaciones sobre una manera muy sencilla de teselar un dominio finito de dos dimensiones: mediante triángulos. Supondremos que $\Gamma \subset \mathbb{R}^2$ es un conjunto cerrad y acotado, cuya frontera $\partial \Gamma$ es una curva con a lo más un número finito de discontinuidades.

Una triangulación es un par ordenado $G,\xi$ tal que $G$ es una gráfica $3-$regular, $\xi : V(G) \to \Gamma$ es un encaje y $\partial \Gamma$ está contenido en todos la unión de todos aristas de $G$ en el plano. En términos más coloquiales, una triangulación es una teselación de una región $\Gamma$ mediante triángulos. En la siguiente figura mostramos un ejemplo de una región triangulada

Dado un número finito de puntos dentro de una región, podemos distinguir a la triangulación de Delauney como la triangulación tal que la circunferencia circunscrita un triangulo no contiene a ningún otro punto más que a los vértices del triángulo. Un ejemplo se muestra en la siguiente figura:

Las triangulaciones de Delauney cumplen con propiedades geométricas importantes. En primer lugar, la gráfica dual corresponde a al diagrama de Voronoi de los vértices de los triángulos. El diagrama de Voronoi tesela al plano en polígonos tales que el vértice de la triangulación más cercano a un punto corresponde al polígono de ese punto. La triangulación de Delauney también maximiza el ángulo mínimo de todos los triángulos de la triangulación.

Existen diversos algoritmos para generar una triangulación de Delauney dado un conjunto de puntos, o también para triangular una región dada un área máxima o alguna otra restricción sobre los triángulos.

Método de elemento finito

Los métodos de elemento finito son un conjunto de métodos empleados principalmente para resolver ecuaciones diferenciales parciales en dominios cuya geometría los vuelve demasiado complicados para ser tratados con el método de elementos finitos. En lugar de discretizar la ecuación diferencial parcial, la función solución a la ecuación se asume como una función que se interpola a partir de los valores en un conjunto discreto de puntos del dominio utilizando una base dada.

En nuestro caso, en cual los puntos del dominio están definidos por una triangulación, la función se aproxima como una función lineal (un plano) en cada triangulo del dominio. Es claro que, dados los valores que toma en tres vértices no colineales de un triángulo, la ecuación de un plano está completamente determinada. La deducción del método escapa a los alcances del presente trabajo y referimos al lector a los capítulos 3 y 4 de la referencia [1] para ver dichos procedimientos.

Túnel de Viento

Las ecuaciones de Navier-Stokes para un fluido incompresible son las siguientes
$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \mathbf{\nabla})\mathbf{u} + \nabla p - \nu \nabla^2 \mathbf{u} = f \\ \hspace{1cm}\\ \\ \mathbf{\nabla} \cdot \mathbf{u} = 0$
Para el caso de un tunel de viento rectangular, con un obstáculo definido por la región $\Gamma_{obs}$, con entrada $\Gamma_{in}$, salida $\Gamma_{out}$ y paredes $\Gamma_{W}$, las condiciones iniciales y a la frontera toman la siguiente forma
$\mathbf{u} \mid_{\Gamma_W} = \mathbf{u} \mid_{\Gamma_{obs}} = \mathbf{0} \\ \hspace{1cm}\\ \mathbf{u} \mid_{\Gamma_{in}} = (U(y) ,0) \\ \hspace{1cm}\\ p \mid_{\Gamma_{out}} = 0$
Por lo general, se asume un perfil parabólico para $U(y)$ de la siguiente forma
$U(y) = \frac{4U_{max} y (L_y - y)}{L_y^2}$
Con $L_y$ la longitud en la dirección $y$ del canal.

Solución numérica

Una manera simple de resolver las ecuaciones de Navier-Stokes mostradas es el método de paso fraccionario primero aproximar la derivada temporal con una diferencia finita para obtener la siguiente ecuación
$\frac{\mathbf{u}_{l+1} - \mathbf{u}_{l}}{k} + (\mathbf{u}_l \cdot \mathbf{\nabla})\mathbf{u}_l + \nabla p_l - \nu \nabla^2 \mathbf{u}_l = f_l$
Con $l$ el índice sobre los tiempos discretos separados a un intervalo $k$. Si ahora añadimos y restamos un $0 = \frac{\mathbf{u}_{*} - \mathbf{u}_{*}}{k}$, obtenemos la siguiente ecuación:
$\frac{\mathbf{u}_{l+1} -\mathbf{u}_{*} + \mathbf{u}_{*} - \mathbf{u}_{l}}{k} + (\mathbf{u}_l \cdot \mathbf{\nabla})\mathbf{u}_l + \nabla p_l - \nu \nabla^2 \mathbf{u}_l = f_l$
Podemos separar la ecuación en dos y obtener que es equivalente:
$\frac{ \mathbf{u}_{*} - \mathbf{u}_{l}}{k} = - (\mathbf{u}_l \cdot \mathbf{\nabla})\mathbf{u}_l + \nu \nabla^2 \mathbf{u}_l + f_l$
Y
$\frac{\mathbf{u}_{l+1} -\mathbf{u}_{*} }{k} = -\nabla p_l$
Al separar de esta forma las ecuaciones de Navier-Stokes podemos trabajar independientemente con la velocidad y la presión. Si volvemos a aplicar el operador $\nabla$ a la ecuación de la presión, entonces se cumple que
$\nabla \cdot \frac{\mathbf{u}_{l+1} -\mathbf{u}_{*} }{k} = -\nabla^2 \cdot p_l$
Suponiendo el flujo $\mathbf{u}_{l+1}$ ya es incompresible, entonces obtenemos
$\nabla \cdot \frac{\mathbf{u}_{*} }{k} = \nabla^2 p_l$
Esta es una ecuación de tipo Poisson para la presión. Así, ya tenemos un posible esquema de solución para las ecuaciones de Navier-Stokes, que se puede plantear de la siguiente manera:

1. Calcular $\mathbf{u}_{*}$ a partir de $\mathbf{u}_{*} = \mathbf{u}_{l}+ k \left( - (\mathbf{u}_l \cdot \mathbf{\nabla})\mathbf{u}_l + \nu \nabla^2 \mathbf{u}_l + f_l \right)$
2. Calcular $p_l$ mediante la ecuación $\nabla \cdot \frac{\mathbf{u}_{*} }{k} = \nabla^2 p_l$
3. Calcular $\mathbf{u}_{l+1}$ mediante $\mathbf{u}_{l+1} = \mathbf{u}_{*} -k\nabla p_l$
   Las condiciones a la frontera se involucran en la solución de cada una de estas ecuaciones, claramente.

En el caso particular de dos dimensiones, con $\mathbf{u} = (u,v)$, para una aproximación con elementos finitos, podemos representar la solución para cada paso del método anterior de la siguiente forma:

1. Calcular $u_{*} =M^{-1}[(C_l + \nu A)u_{l} - b_x]$ y $v_{*} =M^{-1}[(C_l + \nu A)v_{l} - b_y]$
2. Calcular $p_l$ de la ecuación $kAp_l = -(B_x u_* + B_y v_*)$
3. Calcular $u_{l+1} = u_* - k M^{-1}B_xp_l$ y $v_{l+1} = v_* - k M^{-1}B_yp_l$

Debemos aclarar cada término de la ecuación. Supongamos que tenemos una triangulación del dominio a resolver cn $n_p$ puntos y $n_t$ triángulos. Así, $u_i,v_i,p_i$ son vectores columna de dimensión $n_p$ tales que en cada entrada tienen el valor de la velocidad en x, velocidad en y o la presión en dicho punto.

Así, $M$ es la matriz de masas una matriz de $n_p \times n_p$ que contiene la información sobre el peso de cada vértice de las triangulaciones. $A$ es la matriz de rigidez, una matriz de $n_p \times n_p$ análoga a la matriz de masas pero midiendo la variación en las derivadas direccionales de cada punto. $C_l$ es la matriz de convección asociada al campo vectorial $\mathbf{u}_l$. $B_x$ y $B_y$ son matrices de $n_p \times n_p$ representando las derivadas en la dirección $x$ y $y$, respectivamente. Y, por último, $b_x$ y $b_y$ son vectores columna de $n_p$ que contienen la proyección de las fuerzas de cuerpo en cada punto de la triangulación.

Debemos mencionar que cada una de las matrices mencionadas está definida en términos de integrales y proyecciones sobre las funciones base que utilizamos para expresar nuestra función solución. En nuestro caso, las funciones base son las funciones $\phi_i$ funciones lineales en cada triangulo de la triangulación, definidas por tomar los siguientes valores en un punto $x_j$ de la triangulación.
$\phi_i (x_j) = \begin{cases} 1 & \text{si } i =j\\ 0 & \text{si } i \neq j\\ \end{cases}$
El uso de dichas funciones como base y la aproximación de la integral de una función dentro de un triángulo $K$ con vértices $N_1,N_2,N_3$mediante la siguiente fórmula
$\int_{K} f(x,y) \: dx \; dy \approx \frac{|K|}{3} \sum_{i=1}^{3} f(N_i)$
Permite que las expresiones para las matrices mencionadas sean sencillas y se puedan obtener calculando primero las aportaciones de cada triángulo de la triangulación del dominio.

Resultados

Para nuestro programa, utilizamos la geometría del problema DFG: un canal rectangular de dimensión $[0,2.2]\times[0,0.41]$ que contiene un obstáculo, con un valor de $U_{max}=0.3$ para el perfil parabólico. Hay dos posibilidades para el obstáculo que se estudiaron en el presente trabajo:

1. Un círculo de diámetro $0.1$ centrado en $(0.2,0.2)$
2. Un ala de joukowski obtenida de aplicar dicha transformación conforme a un círculo de radio $1.25$, centrado en $(-0.2,0.2)$.

El ala fue escalada a $10 \%$ de su tamaño original y trasladada para que pudiera entrar en el dominio y para tener un ancho sobre el eje y similar a $0.1$. Notemos que el número de Reynolds del sistema está dado por
$Re = \frac{UL}{\nu}= \frac{0.1*0.2}{\nu}$
El paso de tiempo se fijó en $k =0.01$, mientras que, para cualquier obstáculo, la triangulación se construyó tomando $20$ puntos uniformemente separados sobre la superficie del obstáculo y una restricción del area máxima de un triángulo en $0.001$. Las triangulaciones de ambos obstáculos se puden observar en las figuras siguientes:

Las soluciones para $\nu = 0.001$ ($Re = 20$) se pueden observar en las siguientes animaciones. Para el círculo:

Para el ala de Joukowsky:

El comportamiento se vuelve todavía más intesante si tomamos $\nu = 0.00005$ ($Re = 400$). Las siguientes animaciones muestran los resultados para el obstáculo circular:

Para el ala de Joukowsky tenemos:

Es claro que en estos dos casos logramos obesrvar vórtices de Von Karman conforme avanza el tiempo. En general, el flujo se vuelve inestable debido a la baja viscosidad cinética.

Conclusiones

El método de elemento finito nos permitió realizar explorar tanto un problema ya trabajeble con elementos finitos, el canal de viento con un obstáculo circular, como un problema ques difícil de tratar directamente con diferencias finitas: un canal de viento con un obstáculo en forma de ala de Joukowsky.

A su vez, el problema nos permitio modelar el flujo de manera más compleja. Aunque incompresible, el flujo se obtiene resolviendo de manera directa las ecuaciones de Navier-Stokes y, además, exhibe fenómenos no lineales como son los vórtices de Von Karman para un número de Reynolds alto.

Concluimos que el método de elemento finito resulta muy útil para explorar problemas con geometrías no simples.

Referencias

1. Larson, M; Bengzon, F. The Finite Element Method: Theory, Implementations and Applications. Springer. Primera Edición. 2010.