Matriz de rigidez

En el método de elementos finitos para la solución numérica de ecuaciones diferenciales parciales elípticas , la matriz de rigidez es una matriz que representa el sistema de ecuaciones lineales que debe resolverse para determinar una solución aproximada a la ecuación diferencial.

Para simplificar, primero consideraremos el problema de Poisson.

${\displaystyle -\nabla ^{2}u=f}$

en algún dominio Ω , sujeto a la condición de contorno u = 0 en el contorno de Ω . Para discretizar esta ecuación mediante el método de elementos finitos , se elige un conjunto de funciones base { φ ₁ , …, φ _n } definidas en Ω que también se anulan en el contorno. Luego se aproxima

${\displaystyle u\aproximadamente u^{h}=u_{1}\varphi _{1}+\cdots +u_{n}\varphi _{n}.}$

Los coeficientes u ₁ , u ₂ , …, u _n se determinan de manera que el error en la aproximación sea ortogonal a cada función base φ _i :

${\displaystyle \int _{\Omega }\varphi _{i}\cdot f\,dx=-\int _{\Omega }\varphi _{i}\nabla ^{2}u^{h}\,dx=-\sum _{j}\left(\int _{\Omega }\varphi _{i}\nabla ^{2}\varphi _{j}\,dx\right)\,u_{j}=\sum _{j}\left(\int _{\Omega }\nabla \varphi _{i}\cdot \nabla \varphi _{j}\,dx\right)u_{j}.}$

como consecuencia de las condiciones de contorno de Dirichlet homogéneas . La matriz de rigidez es la matriz cuadrada de n elementos A definida por

${\displaystyle \mathbf {A} _{ij}=\int _{\Omega }\nabla \varphi _{i}\cdot \nabla \varphi _{j}\,dx.}$

Al definir el vector de carga F con componentes${\textstyle \mathbf {F} _{i}=\int _{\Omega }\varphi _{i}f\,dx,}$Los coeficientes u _i se determinan mediante el sistema lineal Au = F. La matriz de rigidez es simétrica , es decir, A _ij = A _ji , por lo que todos sus autovalores son reales. Además, es una matriz estrictamente definida positiva , de modo que el sistema Au = F siempre tiene una solución única. (Para otros problemas, estas propiedades ventajosas se pierden).

Cabe destacar que la matriz de rigidez variará según la malla computacional utilizada para el dominio y el tipo de elemento finito empleado. Por ejemplo, la matriz de rigidez, al utilizar elementos finitos cuadráticos por partes, tendrá más grados de libertad que la obtenida con elementos lineales por partes.

La determinación de la matriz de rigidez para otras EDP sigue esencialmente el mismo procedimiento, pero puede complicarse por la elección de las condiciones de contorno. Como ejemplo más complejo, considérese la ecuación elíptica.

${\displaystyle -\sum _{k,l}{\frac {\partial }{\partial x_{k}}}\left(a^{kl}{\frac {\partial u}{\partial x_{l}}}\right)=f}$

dónde${\displaystyle \mathbf {A} (x)=a^{kl}(x)}$es una matriz definida positiva definida para cada punto x en el dominio. Imponemos la condición de contorno de Robin.

${\displaystyle -\sum _{k,l}\nu _{k}a^{kl}{\frac {\partial u}{\partial x_{l}}}=c(ug),}$

donde ν _k es la componente del vector normal unitario exterior ν en la dirección k . El sistema a resolver es

${\displaystyle \sum _{j}\left(\sum _{k,l}\int _{\Omega }a^{kl}{\frac {\partial \varphi _{i}}{\partial x_{k}}}{\frac {\partial \varphi _{j}}{\partial x_{l}}}dx+\int _{\partial \Omega }c\varphi _{i}\varphi _{j}\,ds\right)u_{j}=\int _{\Omega }\varphi _{i}f\,dx+\int _{\partial \Omega }c\varphi _{i}g\,ds,}$

como se puede demostrar utilizando una analogía de la identidad de Green . Los coeficientes u _i se siguen obteniendo resolviendo un sistema de ecuaciones lineales, pero la matriz que representa el sistema es notablemente diferente de la del problema de Poisson ordinario.

En general, a cada operador elíptico escalar L de orden 2 k , se le asocia una forma bilineal B en el espacio de Sobolev H ^k , de modo que la formulación débil de la ecuación Lu = f es

${\displaystyle B[u,v]=(f,v)}$

para todas las funciones v en H ^k . Entonces la matriz de rigidez para este problema es

${\displaystyle \mathbf {A} _{ij}=B[\varphi _{j},\varphi _{i}].}$

Para implementar el método de elementos finitos en una computadora, primero se debe elegir un conjunto de funciones base y luego calcular las integrales que definen la matriz de rigidez. Generalmente, el dominio Ω se discretiza mediante algún método de generación de malla , donde se divide en triángulos o cuadriláteros que no se superponen, denominados elementos. Las funciones base se eligen como polinomios de cierto orden dentro de cada elemento y continuas en sus límites. Las opciones más sencillas son las funciones lineales a trozos para elementos triangulares y las funciones bilineales a trozos para elementos rectangulares.

La matriz de rigidez del elemento A ^{[ k ]} para el elemento T _k es la matriz

${\displaystyle \mathbf {A} _{ij}^{[k]}=\int _{T_{k}}\nabla \varphi _{i}\cdot \nabla \varphi _{j}\,dx.}$

La matriz de rigidez del elemento es cero para la mayoría de los valores de i y j , para los cuales las funciones base correspondientes son cero dentro de T _k . La matriz de rigidez completa A es la suma de las matrices de rigidez de los elementos. En particular, para las funciones base que solo tienen soporte local, la matriz de rigidez es dispersa .

_{Para muchas elecciones estándar de funciones base, es decir ,} funciones base lineales a trozos en triángulos, existen fórmulas simples para las matrices de rigidez de los elementos. Por ejemplo, para elementos lineales a trozos, considérese un triángulo con vértices ₍ x₁ _, y₁ ) _, ( x₂ , y₂ ) , ( x₃ , y₃ ) _y defina la matriz de 2× ₃ .

${\displaystyle \mathbf {D} =\left[{\begin{matrix}x_{3}-x_{2}&x_{1}-x_{3}&x_{2}-x_{1}\\y_{3}-y_{2}&y_{1}-y_{3}&y_{2}-y_{1}\end{matrix}}\right].}$

Entonces la matriz de rigidez del elemento es

${\displaystyle \mathbf {A} ^{[k]}={\frac {\mathbf {D} ^{\mathsf {T}}\mathbf {D} }{4\operatorname {área} (T)}}.}$

Cuando la ecuación diferencial es más complicada, por ejemplo, al tener un coeficiente de difusión no homogéneo, la integral que define la matriz de rigidez del elemento se puede evaluar mediante cuadratura gaussiana .

El número de condición de la matriz de rigidez depende en gran medida de la calidad de la malla numérica. En particular, los triángulos con ángulos pequeños en la malla de elementos finitos generan valores propios elevados en la matriz de rigidez, lo que degrada la calidad de la solución.