Método de los elementos finitos de Galerkin

Como se indicó al principio de este capítulo, pocos tipos de problemas de mecánica de fluidos admiten ser expresados en forma variacional. Sin embargo, la formulación de Galerkin del método de los elementos finitos es en muchos casos equivalente al método de Ritz para resolver problemas variacionales. La mayoría de las aplicaciones del MEF en mecánica de fluidos emplean el método estándar de Galerkin (también denominado método de Bubnow-Galerkin, para distinguirlo del método de Petrov-Galerkin, en el que las funciones de ponderación no coinciden con las de aproximación), por lo que será el que principalmente se considere en lo que sigue.

Para empezar se va a hacer referencia a los dos aspectos más importantes que caracterizan el MEF de Galerkin frente al método de Galerkin de residuos ponderados. En primer lugar, los coeficientes $a_{j}$ que aparecen en la solución aproximada van a coincidir con los valores de las incógnitas en los nodos resultantes de la discretización espacial del dominio de cálculo. La solución aproximada de la ecuación (6.23) se expresa ahora directamente de la forma

$$\tilde{u}(x, y, z, t) = \sum_{j=1}^{M} u_{j}(t)\phi_{j}(x, y, z),$$ (6.45)

es decir, la solución se aproxima en cada punto del dominio mediante una interpolación de los valores de la solución en los nodos, $u_{j}$. Con ello resulta más fácil entender el MEF como un medio de discretizar las ecuaciones en derivadas parciales, lo cual resulta siempre más evidente en métodos de diferencias finitas.

La segunda característica relevante se refiere al tipo de funciones de aproximación que se utilizan, que suelen ser, casi exclusivamente, polinomios de bajo grado restringidos a determinados subdominios o elementos en los que, como se indicará más adelante, se divide el dominio de cálculo, que dan lugar a matrices con pocos elementos distintos de cero. Puede conseguirse además que estos elementos no nulos estén próximos a la diagonal principal, lo que facilitará la resolución del sistema de ecuaciones algebraicas a que da lugar la discretización.

El dominio $\Omega$ debe ser subdividido en subdominios no solapados, llamados elementos, de forma geométrica simple. Las integrales que aparecen en la formulación débil (ecuación (6.42)) pueden expresarse como suma de integrales sobre los elementos.

Una vez realizada la subdivisión del dominio físico en elementos, el problema de representar de forma aproximada la variación de una variable dependiente sobre todo el dominio resulta mucho más sencillo al referir dicha variación a los distintos elementos. Las funciones de aproximación $\phi_{j}$ suelen denominarse funciones de interpolación en la literatura matemática y funciones de forma en ingeniería (utilizando el símbolo $N_{j}$ en lugar de $\phi_{j}$ en este último caso).

En principio, y aunque también pueden emplearse funciones de interpolación cúbicas y de orden superior, en la práctica no es usual, en problemas de mecánica de fluidos, emplear interpolaciones de orden superior al cuadrático. En dinámica de fluidos computacional la eficiencia de los métodos empleados es particularmente importante, sobre todo en problemas tridimensionales con geometrías complicadas y modelos matemáticos complejos, por lo que el compromiso entre economía y precisión suele aconsejar descartar funciones de interpolación de orden elevado, que dan lugar a sistemas de ecuaciones con muchos más términos no nulos. Lógicamente, la interpolación cuadrática se hace progresivamente más precisa que la interpolación lineal a medida que la malla se hace más fina. En principio cabe esperar que, a medida que la malla se hace más fina, el error en la solución numérica se reduzca al mismo ritmo en que lo haga el error de interpolación. Para una cierta malla, generalmente el error en la solución será mayor que el de interpolación al existir un error adicional debido a que la solución en un nodo no coincide con la solución exacta. A modo de idea orientativa, puede decirse que el uso de funciones de forma lineales produce soluciones de precisión similar a la obtenida empleando métodos de diferencias finitas de segundo orden, mientras que el uso de funciones de forma cuadráticas proporciona la misma precisión que los métodos de diferencias finitas de tercer orden.

Una importante característica del MEF es la facilidad con la que permite tratar dominios de cálculo con mallas no uniformes y distorsionadas. Desde que comenzó a aplicarse el MEF se utilizaron muy frecuentemente elementos triangulares debido a que permiten representar formas geométricas complicadas. Sin embargo, también es posible representar dichas formas complejas utilizando elementos rectangulares distorsionados haciendo uso de la transformación isoparamétrica. Además, los elementos rectangulares tienen ventajas sobre los triangulares desde el punto de vista de la utilización de programas de generación de mallas, especialmente en problemas tridimensionales.

Obsérvese que al expresar, de acuerdo con la ecuación (6.45),

$$\tilde{u}(\vec{x},t) = \sum_{j=1}^{M} u_{j}(t)N_{j}(\vec{x})$$ (6.46)

(en lo que sigue se sustituirá $\phi$ por $N$ para denotar las funciones de forma), se está asignando una función de forma a cada valor nodal o grado de libertad. Aunque estas funciones pueden ser de muy diversos tipos, ya se ha indicado que en los métodos estándar de elementos finitos son polinomios definidos localmente, de manera que son nulas fuera de los elementos a los que pertenece el nodo correspondiente. Es decir,

$$N^e_{j}(\vec{x})=0,\;\; \mbox{si $\vec{x}$ no está en $\Omega_e$} ,$$ (6.47)

siendo $N^e_{j}$ la función de forma local correspondiente al nodo $j$ del elemento $e$, que ocupa el subdominio $\Omega_e$. Dado que se ha tomado $\tilde{u}(\vec{x}_j)=u_j$, en cualquier punto $\vec{x}_i$ en el que está localizado un nodo $i$ del elemento $e$ se cumplirá

$$N^e_{j}(\vec{x}_i)=\delta_{ij}.$$ (6.48)

Finalmente, para que pueda representarse de forma exacta una función $\tilde{u}(\vec{x})=$constante sobre del elemento $e$, debe satisfacerse la condición

$$\sum_{j=1}^{n}N^e_{j}(\vec{x})=1,\;\; \mbox{para todo $\vec{x}\in \Omega_e$} ,$$ (6.49)

siendo $n$ el número de nodos en el elemento $e$. Aplicando en cada elemento el esquema de interpolación que se acaba de describir (con las propiedades de las funciones de forma expresadas en las ecuaciones (6.47), (6.48) y (6.49)),

$$\tilde{u} = \sum_{j=1}^{n} u_{j}N^e_{j},$$ (6.50)

se garantiza la continuidad del valor de la función aproximada en todo el dominio. Sumando la interpolación sobre todos los elementos, resulta

$$\tilde{u} = \sum_{e}\sum_{j,e} u_{j}N^e_{j} = \sum_{i=1}^M u_i N_i.$$ (6.51)

En el segundo miembro, el primer sumatorio se extiende sobre todos los elementos del dominio de cálculo, y el segundo, sobre todos los nodos del elemento $e$. Las funciones $N^e_j$, asociadas a cada elemento, se denominan funciones de forma locales o elementales, y las funciones $N_i$, funciones de forma globales. Las funciones de forma globales se obtienen ensamblando las contribuciones $N^e_i$ de todos los elementos a los que pertenece el nodo $i$.

Los elementos y funciones de forma asociadas que generalmente se utilizan pueden clasificarse en dos categorías, dependiendo del grado de continuidad que se requiera entre elementos y de los valores asociados a los nodos. Si los valores nodales se corresponden con los de las variables dependientes y el orden de las ecuaciones de conservación no es superior a dos, es suficiente una continuidad $C^{0}$ en los límites entre elementos. A los elementos con estas características se les suele denominar elementos lagrangianos. Cuando se consideran las derivadas primeras de las variables dependientes como grados de libertad adicionales, generalmente se impone una continuidad $C^1$ en los límites entre elementos. Los elementos con estas características se denominan elementos hermitianos. En los textos de Hirsch (1988) y Wendt (1992), por ejemplo, pueden encontrarse discusiones concisas sobre la elección de funciones de forma y elementos de distintos tipos.

Introduciendo la ecuación (6.45) en la ecuación (6.42), haciendo en ésta $W_i=N_i$, puede deducirse la ecuación discretizada siguiente para el nodo $i$:

$$-\sum_{j\in \Omega_i} u_j \int_{\Omega_i}\lambda\vec{\nabla}N_j\cdot\vec{\nabla}N_i \Omega + \int_{\Omega_i}N_i{Q} \Omega + \int_{\Gamma_{1}}N_i{q} \Gamma=0,$$ (6.52)

donde $\Omega_i$ es el subdominio formado por todos los elementos que contienen el nodo $i$, y el sumatorio sobre $j$ está extendido a todos los nodos del subdominio $\Omega_i$. La matriz formada por los elementos

$$K_{ji}=\int_{\Omega_i}\lambda\vec{\nabla}N_j\cdot\vec{\nabla}N_i \Omega$$ (6.53)

es la matriz de ``rigidez'', que sólo depende de la geometría de la malla y de los elementos utilizados si $\lambda$ es constante.

Obsérvese que la integración que debe hacerse en la formulación débil ha dado lugar a una suma de contribuciones sobre elementos a los que pertenece el nodo considerado. En lugar de esto, sin embargo, también es posible hacer en primer lugar la integración sobre todos los elementos separadamente, y a continuación construir la ecuación para cada nodo a partir de la suma de las contribuciones de los elementos a los que pertenece el nodo, siguiendo un procedimiento de ``ensamblado''. Este ha sido el procedimiento considerado, por ejemplo, en la Sección 6.6.2.