Método de elementos finitos

Un método de elementos finitos se caracteriza por una formulación variacional , una estrategia de discretización, uno o más algoritmos de solución y procedimientos de postprocesamiento.

Ejemplos de la formulación variacional son el método de Galerkin , el método de Galerkin discontinuo, los métodos mixtos, etc.

Una estrategia de discretización se define como un conjunto de procedimientos claramente definidos que abarcan (a) la creación de mallas de elementos finitos, (b) la definición de funciones base en elementos de referencia (también llamadas funciones de forma) y (c) el mapeo de elementos de referencia a los elementos de la malla. Ejemplos de estrategias de discretización son la versión h, la versión p , la versión hp , el método de elementos finitos x , el análisis isogeométrico , etc. Cada estrategia de discretización tiene ciertas ventajas y desventajas. Un criterio razonable para seleccionar una estrategia de discretización es lograr un rendimiento casi óptimo para el conjunto más amplio de modelos matemáticos dentro de una clase de modelos particular.

Los diversos algoritmos de solución numérica se pueden clasificar en dos grandes categorías: solucionadores directos e iterativos. Estos algoritmos están diseñados para aprovechar la escasez de matrices, que depende de la formulación variacional y la estrategia de discretización elegidas.

Los procedimientos de postprocesamiento están diseñados para extraer los datos de interés de una solución de elementos finitos. Para cumplir con los requisitos de verificación de la solución, los postprocesadores deben proporcionar una estimación de error a posteriori en términos de las cantidades de interés. Cuando los errores de aproximación superan lo considerado aceptable, la discretización debe modificarse, ya sea mediante un proceso adaptativo automatizado o por acción del analista. Algunos postprocesadores muy eficientes permiten la superconvergencia .

Los dos problemas siguientes demuestran el método de elementos finitos.

P1 es un problema unidimensional $${\displaystyle {\text{ P1 }}:{\begin{cases}u''(x)=f(x){\text{ en }}(0,1),\\u(0)=u(1)=0,\end{cases}}}$$ dónde${\displaystyle f}$se da,${\displaystyle u}$es una función desconocida de${\displaystyle x}$, y${\displaystyle u''}$es la segunda derivada de${\displaystyle u}$con respecto a${\displaystyle x}$.

P2 es un problema bidimensional ( problema de Dirichlet ). $${\displaystyle {\text{P2 }}:{\begin{cases}u_{xx}(x,y)+u_{yy}(x,y)=f(x,y)&{\text{ en }}\Omega ,\\u=0&{\text{ en }}\partial \Omega ,\end{cases}}}$$

dónde${\displaystyle \Omega }$es una región abierta conectada en el${\displaystyle (x,y)}$plano cuyo límite${\displaystyle \partial \Omega }$es agradable (por ejemplo, una variedad suave o un polígono ), y${\displaystyle u_{xx}}$y${\displaystyle u_{yy}}$denotamos las segundas derivadas con respecto a${\displaystyle x}$y${\displaystyle y}$, respectivamente.

El problema P1 se puede resolver directamente calculando antiderivadas . Sin embargo, este método para resolver el problema de valores en la frontera (PVF) solo funciona cuando hay una dimensión espacial. No se generaliza a problemas de dimensiones superiores o problemas como${\displaystyle u+V''=f}$Por este motivo, desarrollaremos el método de elementos finitos para P1 y describiremos su generalización a P2.

Nuestra explicación se desarrollará en dos pasos, que reflejan dos pasos esenciales que se deben seguir para resolver un problema de valores en la frontera (PVF) utilizando el método de elementos finitos (MEF).

• En el primer paso, se reformula el problema de contorno original en su forma débil. Generalmente, este paso requiere pocos o ningún cálculo. La transformación se realiza a mano sobre papel.
• El segundo paso es la discretización, donde la forma débil se discretiza en un espacio de dimensión finita.

Después de este segundo paso, tenemos fórmulas concretas para un problema lineal grande pero de dimensión finita cuya solución resolverá aproximadamente el problema de contorno original. Este problema de dimensión finita se implementa luego en una computadora .

El primer paso es convertir P1 y P2 en sus formulaciones débiles equivalentes .

Tomamos el intervalo${\displaystyle (0,1)}$, elegir${\displaystyle n}$valores de${\displaystyle x}$con${\displaystyle 0=x_{0}<x_{1}<\cdots <x_{n}<x_{n+1}=1}$y definimos${\displaystyle V}$por: $${\displaystyle V=\{v:[0,1]\to \mathbb {R} \;:v{\text{ is continuous, }}v|_{[x_{k},x_{k+1}]}{\text{ is linear for }}k=0,\dots ,n{\text{, and }}v(0)=v(1)=0\}}$$

donde definimos${\displaystyle x_{0}=0}$y${\displaystyle x_{n+1}=1}$. Observe que las funciones en${\displaystyle V}$no son diferenciables según la definición elemental del cálculo. De hecho, si${\displaystyle v\in V}$entonces la derivada normalmente no está definida en ningún punto.${\displaystyle x=x_{k}}$,${\displaystyle k=1,\ldots ,n}$Sin embargo, la derivada existe en cualquier otro valor de${\displaystyle x}$y se puede utilizar esta derivada para la integración por partes .

Necesitamos${\displaystyle V}$ser un conjunto de funciones de${\displaystyle \Omega }$En la figura de la derecha, hemos ilustrado una triangulación de una región poligonal de 15 lados.${\displaystyle \Omega }$en el plano (abajo), y una función lineal a trozos (arriba, en color) de este polígono que es lineal en cada triángulo de la triangulación; el espacio${\displaystyle V}$Consistiría en funciones que son lineales en cada triángulo de la triangulación elegida.

Se espera que, a medida que la malla triangular subyacente se vuelva cada vez más fina, la solución del problema discreto (3) converja, en cierto sentido, a la solución del problema de contorno original P2. Para medir esta finura de la malla, la triangulación se indexa mediante un parámetro de valor real.${\displaystyle h>0}$que se considera muy pequeño. Este parámetro estará relacionado con el tamaño del triángulo más grande o promedio en la triangulación. A medida que refinamos la triangulación, el espacio de funciones lineales a trozos${\displaystyle V}$También debe cambiar con${\displaystyle h}$Por esta razón, a menudo se lee${\displaystyle V_{h}}$en lugar de${\displaystyle V}$en la literatura. Dado que no realizamos dicho análisis, no utilizaremos esta notación.

Interpolación de una función de Bessel

16 funciones base triangulares escaladas y desplazadas (colores) utilizadas para reconstruir una función de Bessel de orden cero J ₀ (negro)

La combinación lineal de funciones base (en amarillo) reproduce J ₀ (en negro) con la precisión deseada.

Para completar la discretización, debemos seleccionar una base de${\displaystyle V}$En el caso unidimensional, para cada punto de control${\displaystyle x_{k}}$Elegiremos la función lineal por partes.${\displaystyle v_{k}}$en${\displaystyle V}$cuyo valor es${\displaystyle 1}$en${\displaystyle x_{k}}$y cero en cada${\displaystyle x_{j},\;j\neq k}$, es decir, $${\displaystyle v_{k}(x)={\begin{cases}{x-x_{k-1} \over x_{k}\,-x_{k-1}}&{\text{ if }}x\in [x_{k-1},x_{k}],\\{x_{k+1}\,-x \over x_{k+1}\,-x_{k}}&{\text{ if }}x\in [x_{k},x_{k+1}],\\0&{\text{ otherwise}},\end{cases}}}$$

para${\displaystyle k=1,\dots ,n}$; esta base es una función de tienda desplazada y escalada . Para el caso bidimensional, elegimos de nuevo una función base.${\displaystyle v_{k}}$por vértice${\displaystyle x_{k}}$de la triangulación de la región planar${\displaystyle \Omega }$. La función${\displaystyle v_{k}}$es la función única de${\displaystyle V}$cuyo valor es${\displaystyle 1}$en${\displaystyle x_{k}}$y cero en cada${\displaystyle x_{j},\;j\neq k}$.

Según el autor, el término «elemento» en el método de elementos finitos se refiere a los triángulos del dominio, a la función base lineal por partes o a ambas. Por ejemplo, un autor interesado en dominios curvos podría reemplazar los triángulos por primitivas curvas y, por lo tanto, describir los elementos como curvilíneos. Por otro lado, algunos autores reemplazan «lineal por partes» por «cuadrático por partes» o incluso «polinomial por partes». En ese caso, podrían usar «elemento de orden superior» en lugar de «polinomial de grado superior». El método de elementos finitos no se limita a triángulos (tetraedros en 3D o símplexes de orden superior en espacios multidimensionales). También puede definirse en subdominios cuadriláteros (hexaedros, prismas o pirámides en 3D, etc.). Las formas de orden superior (elementos curvilíneos) pueden definirse con formas polinomiales e incluso no polinomiales (por ejemplo, elipses o círculos).

Ejemplos de métodos que utilizan funciones base polinómicas por partes de grado superior son el hp-FEM y el spectral FEM .

Las implementaciones más avanzadas (métodos de elementos finitos adaptativos) utilizan un método para evaluar la calidad de los resultados (basado en la teoría de estimación de errores) y modifican la malla durante la solución con el objetivo de lograr una solución aproximada dentro de ciertos límites de la solución exacta del problema continuo. La adaptabilidad de la malla puede utilizar diversas técnicas; las más populares son:

• nodos móviles (r-adaptabilidad)
• Elementos refinados (y no refinados) (h-adaptabilidad)
• cambio de orden de las funciones base (p-adaptabilidad)
• combinaciones de lo anterior ( hp-adaptabilidad ).

La principal ventaja de esta elección de base es que los productos internos $${\displaystyle \langle v_{j},v_{k}\rangle =\int _{0}^{1}v_{j}v_{k}\,dx}$$ y $${\displaystyle \phi (v_{j},v_{k})=\int _{0}^{1}v_{j}'v_{k}'\,dx}$$ será cero para casi todos${\displaystyle j,k}$. (La matriz que contiene${\displaystyle \langle v_{j},v_{k}\rangle }$en el${\displaystyle (j,k)}$La ubicación se conoce como la matriz de Gram . En el caso unidimensional, el soporte de${\displaystyle v_{k}}$es el intervalo${\displaystyle [x_{k-1},x_{k+1}]}$. Por lo tanto, los integrandos de${\displaystyle \langle v_{j},v_{k}\rangle }$y${\displaystyle \phi (v_{j},v_{k})}$son idénticamente cero siempre que${\displaystyle |j-k|>1}$.

De manera similar, en el caso planar, si${\displaystyle x_{j}}$y${\displaystyle x_{k}}$no comparten un borde de la triangulación, entonces las integrales $${\displaystyle \int _{\Omega }v_{j}v_{k}\,ds}$$ y $${\displaystyle \int _{\Omega }\nabla v_{j}\cdot \nabla v_{k}\,ds}$$ ambos son cero.

Si escribimos${\displaystyle u(x)=\sum _{k=1}^{n}u_{k}v_{k}(x)}$y${\displaystyle f(x)=\sum _{k=1}^{n}f_{k}v_{k}(x)}$entonces problema (3), tomando${\displaystyle v(x)=v_{j}(x)}$para${\displaystyle j=1,\dots ,n}$, se convierte en

Si denotamos por${\displaystyle \mathbf {u} }$y${\displaystyle \mathbf {f} }$los vectores columna${\displaystyle (u_{1},\dots ,u_{n})^{t}}$y${\displaystyle (f_{1},\dots ,f_{n})^{t}}$y si dejamos $${\displaystyle L=(L_{ij})}$$ y $${\displaystyle M=(M_{ij})}$$ sean matrices cuyas entradas sean $${\displaystyle L_{ij}=\phi (v_{i},v_{j})}$$ y $${\displaystyle M_{ij}=\int v_{i}v_{j}dx}$$ entonces podemos reformular (4) como

No es necesario asumir${\displaystyle f(x)=\sum _{k=1}^{n}f_{k}v_{k}(x)}$. Para una función general${\displaystyle f(x)}$, problema (3) con${\displaystyle v(x)=v_{j}(x)}$para${\displaystyle j=1,\dots ,n}$en realidad se vuelve más simple, ya que no hay matriz${\displaystyle M}$se utiliza,

dónde${\displaystyle \mathbf {b} =(b_{1},\dots ,b_{n})^{t}}$y${\displaystyle b_{j}=\int fv_{j}dx}$para${\displaystyle j=1,\dots ,n}$.

Como ya hemos comentado anteriormente, la mayoría de las entradas de${\displaystyle L}$y${\displaystyle M}$son cero porque las funciones base${\displaystyle v_{k}}$tienen poco soporte. Así que ahora tenemos que resolver un sistema lineal en la incógnita${\displaystyle \mathbf {u} }$donde la mayoría de las entradas de la matriz${\displaystyle L}$, que necesitamos invertir, son cero.

Estas matrices se conocen como matrices dispersas , y existen solucionadores eficientes para estos problemas (mucho más eficientes que invertir la matriz directamente). Además,${\displaystyle L}$Es simétrica y definida positiva, por lo que se prefiere una técnica como el método del gradiente conjugado . Para problemas que no son demasiado grandes, las descomposiciones LU dispersas y las descomposiciones de Cholesky siguen funcionando bien. Por ejemplo, el operador de barra invertida de MATLAB (que utiliza LU dispersa, Cholesky dispersa y otros métodos de factorización) puede ser suficiente para mallas con cien mil vértices.

La matriz${\displaystyle L}$se suele denominar matriz de rigidez , mientras que la matriz${\displaystyle M}$se la denomina matriz de masas .

En general, el método de elementos finitos se caracteriza por el siguiente proceso.

• Uno elige una cuadrícula para${\displaystyle \Omega }$En el tratamiento anterior, la cuadrícula estaba formada por triángulos, pero también se pueden utilizar cuadrados o polígonos curvilíneos.
• A continuación, se eligen las funciones base. En nuestra explicación utilizamos funciones base lineales a trozos, pero es común usar funciones base polinómicas a trozos.

Una consideración aparte es la suavidad de las funciones base. Para problemas de contorno elípticos de segundo orden , basta con una función base polinómica por partes que sea simplemente continua (es decir, las derivadas son discontinuas). Para ecuaciones diferenciales parciales de orden superior, se deben usar funciones base más suaves. Por ejemplo, para un problema de cuarto orden como${\displaystyle u_{xxxx}+u_{yyyy}=f}$, se pueden utilizar funciones base cuadráticas por partes que son${\displaystyle C^{1}}$.

Otra consideración es la relación del espacio de dimensión finita${\displaystyle V}$a su contraparte de dimensión infinita en los ejemplos anteriores${\displaystyle H_{0}^{1}}$. Un método de elemento conforme es aquel en el que el espacio${\displaystyle V}$es un subespacio del espacio de elementos para el problema continuo. El ejemplo anterior es un método de este tipo. Si no se cumple esta condición, obtenemos un método de elementos no conformes , un ejemplo del cual es el espacio de funciones lineales a trozos sobre la malla, que son continuas en cada punto medio de los bordes. Dado que estas funciones son generalmente discontinuas a lo largo de los bordes, este espacio de dimensión finita no es un subespacio del original.${\displaystyle H_{0}^{1}}$.

Normalmente, se dispone de un algoritmo para subdividir una malla dada. Si el método principal para aumentar la precisión es subdividir la malla, se dispone de un método h ( donde h suele ser el diámetro del elemento más grande de la malla). De esta manera, si se demuestra que el error con una cuadrícula${\displaystyle h}$está delimitado superiormente por${\displaystyle Ch^{p}}$, para algunos${\displaystyle C<\infty }$y${\displaystyle p>0}$, entonces se tiene un método de orden p . Bajo hipótesis específicas (por ejemplo, si el dominio es convexo), un polinomio por partes de orden${\displaystyle d}$El método tendrá un error de orden${\displaystyle p=d+1}$.

Si en lugar de disminuir h , se aumenta el grado de los polinomios utilizados en la función base, se obtiene un método p . Si se combinan estos dos tipos de refinamiento, se obtiene un método hp ( hp-FEM ). En el hp-FEM, los grados de los polinomios pueden variar de un elemento a otro. Los métodos de orden superior con un p uniforme grande se denominan métodos de elementos finitos espectrales ( SFEM ). Estos no deben confundirse con los métodos espectrales .

Para ecuaciones diferenciales parciales vectoriales, las funciones base pueden tomar valores en${\displaystyle \mathbb {R} ^{n}}$.

El método de elementos aplicados o AEM combina características tanto del método de elementos finitos (FEM) como del método de elementos discretos (DEM).

Yang y Lui introdujeron el Método de Elementos Finitos Aumentados, cuyo objetivo era modelar las discontinuidades débiles y fuertes sin necesidad de grados de libertad adicionales, como afirmó PuM.

El método de elementos finitos de corte se desarrolló en 2014. ^{[ 16 ]} El método consiste en "hacer que la discretización sea lo más independiente posible de la descripción geométrica y minimizar la complejidad de la generación de la malla, manteniendo al mismo tiempo la precisión y robustez de un método de elementos finitos estándar". ^{[ 17 ]}

El método generalizado de elementos finitos (GFEM) utiliza espacios locales formados por funciones, no necesariamente polinomios, que reflejan la información disponible sobre la solución desconocida y, por lo tanto, garantizan una buena aproximación local. Posteriormente, se utiliza una partición de la unidad para "unir" estos espacios y formar el subespacio de aproximación. La eficacia del GFEM se ha demostrado al aplicarlo a problemas con dominios de contornos complejos, problemas con microescalas y problemas con capas límite. ^{[ 18 ]}

El método de elementos finitos mixtos es un tipo de método de elementos finitos en el que se introducen variables independientes adicionales como variables nodales durante la discretización de un problema de ecuaciones diferenciales parciales.

El hp-FEM combina de forma adaptativa elementos con tamaño variable h y grado polinómico p para lograr tasas de convergencia exponenciales excepcionalmente rápidas. ^{[ 19 ]}

El hpk-FEM combina de forma adaptativa elementos con tamaño variable h , grado polinomial de las aproximaciones locales p y diferenciabilidad global de las aproximaciones locales ( k -1) para lograr las mejores tasas de convergencia.

El método de elementos finitos extendido (XFEM) es una técnica numérica basada en el método de elementos finitos generalizado (GFEM) y el método de partición de la unidad (PUM). Extiende el método clásico de elementos finitos al enriquecer el espacio de soluciones para ecuaciones diferenciales con funciones discontinuas. Los métodos de elementos finitos extendidos enriquecen el espacio de aproximación para reproducir de forma natural las características problemáticas asociadas al problema en cuestión: la discontinuidad, la singularidad, la capa límite, etc. Se ha demostrado que, para algunos problemas, esta incorporación de las características del problema en el espacio de aproximación puede mejorar significativamente las tasas de convergencia y la precisión. Además, el tratamiento de problemas con discontinuidades mediante XFEM elimina la necesidad de mallar y remallar las superficies de discontinuidad, lo que reduce los costos computacionales y los errores de proyección asociados a los métodos de elementos finitos convencionales, a costa de restringir las discontinuidades a los bordes de la malla.

Varios códigos de investigación implementan esta técnica en distintos grados:

1. GetFEM++
2. xfem++
3. openxfem++

XFEM también se ha implementado en códigos como Altair Radios, ASTER, Morfeo y Abaqus. Su adopción está creciendo entre otros programas comerciales de elementos finitos, con varios complementos e implementaciones básicas disponibles (ANSYS, SAMCEF, OOFELIE, etc.).

La introducción del método de elementos finitos de contorno escalado (SBFEM) provino de Song y Wolf (1997). ^{[ 20 ]} El SBFEM ha sido una de las contribuciones más provechosas en el área del análisis numérico de problemas de mecánica de fractura. Es un método semianalítico sin solución fundamental que combina las ventajas de las formulaciones y procedimientos de elementos finitos y la discretización de elementos de contorno. Sin embargo, a diferencia del método de elementos de contorno, no requiere una solución diferencial fundamental.

El método S-FEM (Métodos de Elementos Finitos Suavizados) es una clase particular de algoritmos de simulación numérica para la simulación de fenómenos físicos. Se desarrolló combinando métodos sin malla con el método de elementos finitos.

Los métodos de elementos espectrales combinan la flexibilidad geométrica de los elementos finitos con la precisión de los métodos espectrales. Estos métodos son la solución aproximada de ecuaciones parciales en forma débil, basadas en interpolantes lagrangianos de alto orden, y se utilizan únicamente con ciertas reglas de cuadratura. ^{[ 21 ]}

La iteración de Loubignac es un método iterativo en los métodos de elementos finitos.

El método de elementos finitos de plasticidad cristalina (CPFEM) es una herramienta numérica avanzada desarrollada por Franz Roters. Los metales pueden considerarse agregados cristalinos que presentan anisotropía bajo deformación, como la localización anómala de tensiones y deformaciones. El CPFEM, basado en el deslizamiento (velocidad de deformación por cizallamiento), permite calcular la dislocación, la orientación cristalina y otra información textural para considerar la anisotropía cristalina durante el proceso. Se ha aplicado en el estudio numérico de la deformación de materiales, la rugosidad superficial, las fracturas, etc.

El método de elementos virtuales (MEV), introducido por Beirão da Veiga et al. (2013) ^{[ 22 ]} como una extensión de los métodos de diferencias finitas miméticas (MFD), es una generalización del método de elementos finitos estándar para geometrías de elementos arbitrarias. Esto permite la admisión de polígonos generales (o poliedros en 3D) que son altamente irregulares y no convexos. El nombre " virtual" deriva del hecho de que no se requiere el conocimiento de la base de la función de forma local y, de hecho, nunca se calcula explícitamente.