Another FEM: Rectangular element of 4 nodes

Resumen

Tanto en el tradicional MEF como en la metodología que se propone en este artículo, las Ecuaciones de Compatibilidad se consiguen a partir de las mismas funciones de interpolación, por lo que la matriz Elemental que relaciona las deformaciones en los nodos del elemento con los desplazamientos de dichos nodos es idéntica en ambos procedimientos. Pero, mientras que en aquél se aplica el PTV para la determinación de las Fuerzas Nodales Equivalentes que se han de utilizar para la consecución de las Ecuaciones de Equilibrio de todos los nodos de la estructura, en éste dichas fuerzas se consiguen bajo la hipótesis de que las tensiones en los lados del elemento varían linealmente a lo largo de dichos lados. Además de esto, puesto que el sistema de ecuaciones que en este procedimiento se utiliza viene dado en todas las incógnitas del problema en forma explícita, es posible imponer a cualquiera de dichas incógnitas todo tipo de restricción, por lo que, aparte de las condiciones de apoyo imprescindibles para evitar el movimiento como sólido rígido, se imponen también condiciones de equilibrio a todos y cada uno de los elementos que discretizan la estructura, así como las condiciones de tensión en el perímetro del dominio; condiciones que, al no ser posible contemplar con el tradicional MEF, éste se ve obligado a sumir unos errores originados por un mayor alejamiento entre la modelización a adoptar y el hecho físico real. Los ejemplos prácticos que se estudian, utilizando como base el elemento rectangular de 4 nodos, muestran unos resultados un tanto mejorados respecto a los que se obtienen con el tradicional MEF.

Palabras clave: Elemento, matriz, compatibilidad, equilibrio, constitutivas, fuerzas nodales

Abstract

In both of the traditional FEM and the proposed methodology in this article, the compatibility equations are obtained from the same interpolation functions, so the elementary matrix that relates deformations in the nodes of the element with the displacements of that nodes are the same in both procedures. But, whereas in that one the PVW is applied to establish the Equivalent Nodal Forces to be used for the achievement of the equilibrium equations of all the nodes of the structure, in which these forces are obtained from the hypothesis that the stresses on all four sides of the elemental rectangle vary linearly along these sides and are thus to replace said stresses as statically equivalent, concentrated forces. Since the system of equations used in this procedure is given in all the unknowns of the problem explicitly, and it is possible to impose any kind of restriction on any of these unknowns, so that, in addition to the conditions of essential support to avoid movement as a solid rigid, conditions of equilibrium are imposed on every one of the elements that discretize the structure, as well as the conditions of tension in the perimeter of the domain; conditions that, since it is not possible to be considered with the traditional FEM, this one is forced to sink some errors originated by a greater distance between the modeling to be adopted and the physical real fact. The practical examples studied, using as a base the rectangular element of 4 nodes, show somewhat improved results with respect to those obtained with the traditional FEM.

Keywords: Matrix, element, compatibility, balance, constitutive, nodal forces.

1. Antecedentes

El cálculo de estructuras por el Método de los Elementos Finitos tradicional (MEF) está fundamentalmente sustentado en, por una parte, la elección de unas determinadas funciones (funciones de interpolación) y, por otra, en la aplicación del Principio de los Trabajos Virtuales (PTV). Pero, comoquiera que dichas funciones no pueden reproducir de forma exacta el comportamiento mecánico de la estructura en estudio y todas las magnitudes involucradas en el proceso de cálculo dependen de forma directa de tales funciones, los resultados que se obtendrán (desplazamientos, tensiones, etc.) sólo podrán ser más o menos aproximados; aproximación que, además de las causas ya apuntadas, se verá afectada por la integración numérica utilizada habitualmente.

Nota: en adelante se hará una distinción entre nodos y nudos, considerando nodos a los vértices de los elementos (que vendrán signados por i, j, k… para cada elemento e), los cuales estarán sometidos a tensiones y deformaciones, mientras que los nudos (signados por un subíndice numérico en la forma n₁, n_2, n₃...) son los puntos donde inciden uno o más nodos y sólo pueden sufrir desplazamientos.

2. Objetivos

Para el análisis de toda estructura formada por elementos discretos conectados entre sí en unos puntos que llamaremos nudos, es imprescindible el uso de cuatro tipos de ecuaciones básicas: a) Ecuaciones de Equilibrio de fuerzas que inciden en cada nudo, b) Ecuaciones de Compatibilidad que obliguen a que se cumplan unas determinadas relaciones entre los desplazamientos de los nudos y las deformaciones que sufren los elementos, c) Ecuaciones Constitutivas del material, que relacionan las deformaciones con las tensiones y d) Ecuaciones que eviten el movimiento de la estructura como sólido rígido, esto es, las Condiciones Cinemáticas de sustentación. El tratamiento de este último tipo de ecuaciones por el Método de los Multiplicadores de Lagrange permite obtener, junto con los tres tipos anteriores, un sistema en todas las incógnitas involucradas en el análisis de la estructura. Habitualmente se suelen reducir las dimensiones de tal sistema mediante la aplicación del PTV, obteniendo otro menor en las incógnitas de desplazamientos y de Multiplicadores de Lagrange (reacciones con sentido opuesto); sin embargo, en el procedimiento que en este trabajo se presenta, se mantienen de forma explícita todas las ecuaciones mencionadas, lo que permitirá introducir cualquier tipo de restricción a cualquiera de las incógnitas involucradas, por lo que puede obligarse a que se cumpla el equilibrio en cada elemento aislado sin más que añadir al sistema anterior un quinto tipo de ecuaciones: las Condiciones de Restricción Estática interna y, finalmente, un último grupo de ecuaciones que describan el estado tensional del perímetro del dominio: las Condiciones de Restricción estática externa o perimetral.

Es, pues, objetivo de este artículo la descripción de un nuevo procedimiento aplicado al cálculo por elementos finitos de estructuras en el espacio 2D, tomando como base el elemento rectangular de 4 nodos y utilizando el sistema de ecuaciones antes descrito (primera diferencia respecto al método clásico, excepto las ecuaciones de compatibilidad, que son idénticas), sistema que contempla el equilibrio de cada elemento aislado (segunda diferencia fundamental) así como el estado tensional del perímetro del dominio (tercera diferencia) y sin la utilización de integración alguna (cuarta diferencia), lo que supone una formulación totalmente distinta a la usada hasta ahora.

3. Herramientas del método

El procedimiento que aquí se propone es un método matricial, por lo tanto las unidades operativas que maneja son matrices Elementales agrupadas en seis tipos, en directa correspondencia con las ecuaciones mencionadas en el Apartado 2.: a) matrices de Fuerzas Nodales, b) matrices de Compatibilidad, c) matrices Constitutivas, d) matrices de Restricciones Cinemáticas (restricciones en aquellos aparatos de apoyo que sustentan a la estructura), e) matrices de Restricciones Estáticas internas (condiciones de equilibrio de cada uno de los elementos finitos con los que se ha discretizado ésta) y f) matrices de Restricciones Estáticas externas (estado tensional del perímetro del dominio).

3.1. Matrices Elementales de Fuerzas Nodales

Son aquellas matrices que permiten determinar las dos componentes de fuerza concentrada (X, Y), en cada uno de los nodos del elemento (cuatro para el elemento que se ha elegido en esta presentación) en función de las tensiones existentes en dichos puntos. Coherentemente con las función de interpolación adoptada, las tensiones variarán linealmente a lo largo de cada uno de los lados del rectángulo, de dimensiones 2a×2b y espesor t, y por tanto la distribución será entones como se muestra en las Figuras 1 (a) y (b), donde se ha tomado para dichas tensiones el criterio convencional de signos y se ha supuesto que el elemento no contiene cargas volumétricas ni de superficie.

Figura 1. Elemento rectangular de 4 nodos aislado. (a) Distribución de tensiones normales. b) Tensiones tangenciales. c) Fuerzas Nodales Equivalentes

Signando por i, j, k y m a los nodos del elemento, las componentes de Fuerza Nodal Equivalente en el nodo i (Figura 1 (c)) podrán obtenerse en la forma ^[1]§7.3.2.

${\begin{array}{c}X_{i}=-t{\frac {2b}{6}}(2{\sigma }_{x}^{i}+{\sigma }_{x}^{m})-t{\frac {2a}{6}}(2{\tau }^{i}+{\tau }^{j})=-{\frac {t}{3}}(2b{\mbox{ }}{\sigma }_{x}^{i}+b{\mbox{ }}{\sigma }_{x}^{m}+2a{\mbox{ }}{\tau }^{i}+a{\mbox{ }}{\tau }^{j})\\Y_{i}=-t{\frac {2a}{6}}(2{\sigma }_{y}^{i}+{\sigma }_{y}^{j})-t{\frac {2b}{6}}(2{\tau }^{i}+{\tau }^{m})=-{\frac {t}{3}}(2a{\mbox{ }}{\sigma }_{y}^{i}+a{\mbox{ }}{\sigma }_{y}^{j}+2b{\mbox{ }}{\tau }^{i}+b{\mbox{ }}{\tau }^{m})\end{array}}$

(1)

donde X_i e Y_i son, respectivamente, la componente horizontal y vertical de la fuerza concentrada en el nodo i, fuerza que sustituye a las tensiones distribuidas a lo largo de los lados i-j e i-m del elemento. Llamando σ^e al vector columna que contiene las componentes de tensión existentes en los cuatro vértices del elemento (en el orden i, j, k, m), éste vendrá expresado como

${\mathbf {\mbox{σ}} ^{\mbox{e}}}={\left[{\begin{array}{ccc|ccc|ccc|ccc}{\sigma }_{x}^{i}&{\sigma }_{y}^{i}&{\tau }^{i}&{\sigma }_{x}^{j}&{\sigma }_{y}^{j}&{\tau }^{j}&{\sigma }_{x}^{k}&{\sigma }_{y}^{k}&{\tau }^{k}&{\sigma }_{x}^{m}&{\sigma }_{y}^{m}&{\tau }^{m}\end{array}}\right]}^{T}$

(2)

por lo que las componentes dadas en la ec.(1) se podrán escribir matricialmente en la forma

$\left[{\begin{array}{c}X_{i}\\Y_{i}\end{array}}\right]={\frac {t}{3}}\left[{\begin{array}{ccc|ccc|ccc|ccc}-2b&0&-2a&0&0&-a&0&0&0&-b&0&0\\0&-2a&-2b&0&-a&0&0&0&0&0&0&-b\end{array}}\right]\left[{\mbox{σ}}^{\mbox{e}}\right]$

(3)

y signando por ${\mbox{F}}_{\mbox{i}}^{\mbox{e}},$${\mbox{F}}_{\mbox{j}}^{\mbox{e}},$${\mbox{F}}_{\mbox{k}}^{\mbox{e}}$ y ${\mbox{F}}_{\mbox{m}}^{\mbox{e}}$ a cada uno de los vectores que contienen a las dos componentes de fuerza nodal de los nudos i, j, k y m, respectivamente, la anterior vendrá dada para el nodo i por

$\mathbf {{\mbox{F}}_{\mbox{i}}^{\mbox{e}}} =\mathbf {{\mbox{I}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}}$

(4)

donde se ha llamado

$\mathbf {{\mbox{I}}^{\mbox{e}}} ={\frac {t}{3}}\left[{\begin{array}{ccc|ccc|ccc|ccc}-2b&0&-2a&0&0&-a&0&0&0&-b&0&0\\0&-2a&-2b&0&-a&0&0&0&0&0&0&-b\end{array}}\right]$

(5)

a la matriz de Fuerza Nodal asociada al nodo i, que, como se ha visto, permite calcular las componentes de dicha fuerza en función de las dimensiones de los lados del rectángulo y de las tensiones a las que están sometidos los cuatro nodos del elemento.

Procediendo de la misma forma para los nodos j, k y m, las expresiones para las correspondientes matrices de Fuerza Nodal para estos nodos serán

$\mathbf {{\mbox{J}}^{\mbox{e}}} ={\frac {t}{3}}\left[{\begin{array}{ccc|ccc|ccc|ccc}0&0&-a&2b&0&-2a&b&0&0&0&0&0\\0&-a&0&0&-2a&2b&0&0&b&0&0&0\end{array}}\right]$

(6)

$\mathbf {{\mbox{K}}^{\mbox{e}}} ={\frac {t}{3}}\left[{\begin{array}{ccc|ccc|ccc|ccc}0&0&0&b&0&0&2b&0&2a&0&0&a\\0&0&0&0&0&b&0&2a&2b&a&0&0\end{array}}\right]$

(7)

$\mathbf {{\mbox{M}}^{\mbox{e}}} ={\frac {t}{3}}\left[{\begin{array}{ccc|ccc|ccc|ccc}-b&0&0&0&0&0&0&0&a&-2b&0&2a\\0&0&-b&0&0&0&0&a&0&0&2a&-2b\end{array}}\right]$

(8)

Las expresiones dadas en las ecs.(5)-(8) son las unidades operativas que llamaremos matrices Elementales de Fuerzas Nodales, con cuyo pertinente ensamblaje se obtendrá una matriz H, matriz de Equilibrio de la Estructura, originándose de esta forma un primer conjunto de ecuaciones lineales simultáneas en las incógnitas de tensiones de todos los nodos de la estructura, que obligará a que se cumpla el equilibrio estático en cada uno de sus nudos, siendo el vector de términos independientes de este subsistema las acciones que actúan sobre todos ellos, incluidos los que pertenezcan a los aparatos de apoyo.

3.2. Matrices Elementales de Compatibilidad

Los valores de las deformaciones y de los desplazamientos no pueden ser independientes entre sí, por lo que necesariamente ha de existir un conjunto de ecuaciones que los relacionen. En este trabajo se adoptan las mismas relaciones que habitualmente se han tomado desde la aparición del MEF, las cuales son consecuencia de la adopción de unas determinadas funciones de interpolación con las que se obtienen otras funciones (funciones de forma) de las que se deriva finalmente la conocida expresión

$\mathbf {{\mbox{ε}}_{\mbox{(x,y)}}^{\mbox{e}}} =\mathbf {{\mbox{B}}_{\mbox{(x,y)}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{d}}_{\mbox{(x,y)}}^{\mbox{e}}}$

(9)

donde $\mathbf {{\mbox{ε}}_{\mbox{(x,y)}}^{\mbox{e}}}$ y $\mathbf {{\mbox{d}}_{\mbox{(x,y)}}^{\mbox{e}}}$ son, respectivamente, los vectores columna de deformación y de desplazamiento de un punto genérico (x,y) del elemento, mientras que la matriz que los relaciona, Matriz de Compatibilidad para el elemento rectangular que nos ocupa, viene dada por la conocida expresión

$\mathbf {{\mbox{B}}_{\mbox{(x,y)}}^{\mbox{e}}} ={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc|cc|cc|cc}y-b&0&b-y&0&b+y&0&-b-y&0\\0&x-a&0&-a-x&0&a+x&0&a-x\\x-a&y-b&-a-x&b-y&a+x&b+y&a-x&-b-y\end{array}}\right]$

(10)

siendo el origen de coordenadas el centro de gravedad del elemento, de acuerdo con la Figura 1. Particularizando la expresión (10) para cada uno de los cuatro nodos del elemento, la (9) desarrollada tendrá la forma

$\left[{\begin{array}{c}{\mbox{ε}}_{x}^{i}\\{\mbox{ε}}_{y}^{i}\\{\gamma }^{i}\\{\mbox{ε}}_{x}^{j}\\{\mbox{ε}}_{y}^{j}\\{\gamma }^{j}\\{\mbox{ε}}_{x}^{k}\\{\mbox{ε}}_{y}^{k}\\{\gamma }^{k}\\{\mbox{ε}}_{x}^{m}\\{\mbox{ε}}_{y}^{m}\\{\gamma }^{m}\end{array}}\right]={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc|cc|cc|cc}-b&0&b&0&0&0&0&0\\0&-a&0&0&{\mbox{ }}{\mbox{ }}0{\mbox{ }}{\mbox{ }}&{\mbox{ }}{\mbox{ }}0{\mbox{ }}{\mbox{ }}&0&a\\-a&-b&0&b&0&0&a&0\\\hline -b&0&b&0&0&0&0&0\\0&0&0&-a&0&a&0&0\\0&-b&-a&b&a&0&0&0\\\hline 0&0&0&0&b&0&-b&0\\0&0&0&-a&0&a&0&0\\0&0&-a&0&a&b&0&-b\\\hline 0&0&0&0&b&0&-b&0\\0&-a&0&0&0&0&0&a\\-a&0&0&0&0&b&a&-b\end{array}}\right]\,\,\,\,\left[{\begin{array}{c}u_{i}\\v_{i}\\\hline u_{j}\\v_{j}\\\hline u_{k}\\v_{k}\\\hline u_{m}\\v_{m}\end{array}}\right]{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\acute {o}}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{ε}}^{\mbox{e}}} {\mbox{=}}\mathbf {{\mbox{B}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{d}}^{\mbox{e}}}$

(11)

que permite el cálculo de las deformaciones ε^e en los cuatro nodos en función de los desplazamientos d^e de ellos. Particionando la matriz B^e anterior se obtienen las expresiones

$\mathbf {{\mbox{i}}^{\mbox{e}}} ={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc}-b&0\\0&-a\\-a&-b\\\hline -b&0\\0&0\\0&-b\\\hline 0&0\\0&0\\0&0\\\hline 0&0\\0&-a\\-a&0\end{array}}\right];$${\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{j}}^{\mbox{e}}} ={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc}b&0\\0&0\\0&b\\\hline b&0\\0&-a\\-a&b\\\hline 0&0\\0&-a\\-a&0\\\hline 0&0\\0&0\\0&0\end{array}}\right];$${\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{k}}^{\mbox{e}}} ={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc}{\mbox{ }}{\mbox{ }}0{\mbox{ }}{\mbox{ }}&{\mbox{ }}{\mbox{ }}0{\mbox{ }}{\mbox{ }}\\0&0\\0&0\\\hline 0&0\\0&a\\a&0\\\hline b&0\\0&a\\a&b\\\hline b&0\\0&0\\0&b\end{array}}\right];$${\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{m}}^{\mbox{e}}} ={\frac {1}{2a{\mbox{ }}b}}\left[{\begin{array}{cc}0&0\\0&a\\a&0\\\hline 0&0\\0&0\\0&0\\\hline -b&0\\0&0\\0&-b\\\hline -b&0\\0&a\\a&-b\end{array}}\right];$

(12)

que llamaremos matrices Elementales de Compatibilidad asociadas a cada uno de los nodos del elemento. La primera de ellas determina las deformaciones que se producirán en los cuatro nodos como consecuencia de los desplazamientos del nodo i, la segunda las deformaciones originadas en dichos nodos debidas a los desplazamientos del j, etc. De forma semejante a las matrices de Fuerzas Nodales deducidas en el apartado anterior, las dadas en la ec.(12) no son más que las unidades operativas con cuyo ensamblaje se obtendrá otra matriz B, matriz de Compatibilidad de la Estructura, originándose un segundo conjunto de ecuaciones lineales simultáneas y homogéneas en las incógnitas de desplazamientos d de los nudos y de las deformaciones ε de los nodos que obligarán al cumplimiento de la compatibilidad entre éstas y aquéllos en todo el sistema estructural.

3.3. Matriz Elemental Constitutiva

Como es bien sabido, la relación existente entre las tensiones y deformaciones en un nodo genérico n viene dada por la ley de Hooke generalizada

$\mathbf {{\mbox{σ}}_{\mbox{n}}} =\mathbf {{\mbox{D}}_{\mbox{n}}} {\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}}$

(13)

donde $\mathbf {{\mbox{σ}}_{\mbox{n}}}$ y $\mathbf {{\mbox{ε}}_{\mbox{n}}}$ son, respectivamente, los vectores columna de tensión y de deformación en dicho nodo y, ciñéndonos exclusivamente al estado de Tensión Plana, la Matriz Elemental Constitutiva viene dada por

$\mathbf {{\mbox{D}}_{\mbox{n}}} =E\left[{\begin{array}{ccc}f_{1}&f_{2}&0\\f_{2}&f_{1}&0\\0&0&f_{3}\end{array}}\right];{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}siendo{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}f_{1}=$${\frac {1}{1-{\nu }^{2}}};{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}f_{2}=$$\nu {\mbox{ }}f_{1};{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}f_{3}=$${\frac {1}{2(1+\nu )}};$

(14)

en donde E es el Módulo de Elasticidad y ν la relación de Poisson del material. Aunque puede operarse con estas expresiones, la introducción del valor numérico de E en la ecuación (13) suele acarrear algún problema de condicionamiento por ser su orden de magnitud muy superior a los del resto de valores involucrados en el proceso de cálculo, por lo que efectuando el cambio de variable $\mathbf {{\mbox{ε}}{\mbox{'}}_{\mbox{n}}} =E{\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}}$ la ec.(13) vendrá dada en la forma

$\mathbf {{\mbox{σ}}_{\mbox{n}}} =\mathbf {{\mbox{K}}_{\mbox{n}}} {\mbox{ }}E{\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}} {\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}es{\mbox{ }}{\mbox{ }}decir{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{σ}}_{\mbox{n}}} =$$\mathbf {{\mbox{K}}_{\mbox{n}}} {\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}^{'}}$

(15)

siendo

$\mathbf {{\mbox{K}}_{\mbox{n}}} =\left[{\begin{array}{ccc}f_{1}&f_{2}&0\\f_{2}&f_{1}&0\\0&0&f_{3}\end{array}}\right]$

(16)

quedando de esta forma eliminado dicho módulo. De acuerdo con esto, y teniendo en cuenta que el elemento que nos ocupa tiene cuatro nodos, llamaremos matriz Elemental Constitutiva K^e a la matriz diagonal por bloques de dimensiones $4\times 4$

$\mathbf {K} ^{\mathbf {e} }=\left[{\begin{array}{cccc}\mathbf {K} _{\mathbf {n} }&\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &\mathbf {K} _{\mathbf {n} }&\mathbf {0} &\mathbf {0} \\\mathbf {0} &\mathbf {0} &\mathbf {K} _{\mathbf {n} }&\mathbf {0} \\\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {K} _{\mathbf {n} }\end{array}}\right]$

(17)

unidad operativa con cuyo ensamblaje se obtendrá una tercera matriz K, matriz Constitutiva de la Estructura, y, por tanto, un tercer sistema de ecuaciones simultáneas lineales, asimismo homogéneas, en las incógnitas de tensiones σ y deformaciones ε' de todos los nodos de la estructura en estudio. Para introducir en la ec.(11) dicho cambio de variable, bastará con multiplicar m.a.m. esta ecuación por E, con lo que la ecuación de compatibilidad para el elemento quedará en la forma

$\mathbf {{\mbox{ε}}{\mbox{'}}^{\mbox{e}}} {\mbox{=}}\mathbf {{\mbox{B}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{d}}{\mbox{'}}^{\mbox{e}}}$

(18)

donde $\mathbf {{\mbox{ε}}{\mbox{'}}^{\mbox{e}}} =E{\mbox{ }}\mathbf {{\mbox{ε}}^{\mbox{e}}}$ y $\mathbf {{\mbox{d}}{\mbox{'}}^{\mbox{e}}} =E{\mbox{ }}\mathbf {{\mbox{d}}^{\mbox{e}}} .$

3.4. Matrices Elementales de Restricciones Cinemáticas

Con los tres grupos de ecuaciones indicadas anteriormente, se puede formar un sistema de ecuaciones simultáneas y lineales para toda la estructura en las incógnitas de desplazamientos $\mathbf {\mbox{d'}} ,$ tensiones σ y deformaciones $\mathbf {\mbox{ε'}} ,$ cuya matriz de coeficientes será cuadrada, pero singular, puesto que aún no se han introducido las condiciones de sustentación que eviten el movimiento de la estructura como sólido rígido. Estas condiciones-ecuaciones cinemáticas, que por brevedad supondremos puntuales y homogéneas, aunque pueden ser de cualquier otro tipo, tendrán aquí un tratamiento basado en los Multiplicadores de Lagrange, con lo cual no sólo se evitará la reordenación y eliminación de filas y columnas para proceder a la resolución del sistema resultante (como es habitual en el tradicional MEF), sino que de esta forma se obtendrá directamente el valor de las reacciones en los aparatos de apoyo, pues como se sabe, dichos multiplicadores no son más que éstas con el signo contrario. Aunque las unidades matriciales operativas necesarias para la determinación de tales ecuaciones de restricción pueden ser de varios tipos (piénsese, p.e., en condiciones cinemáticas no concordantes, condiciones multipunto, etc.) en este trabajo nos ceñiremos exclusivamente a las más habituales de desplazamientos nulos en los aparatos de apoyo. Dichas condiciones vendrán dadas en formato matricial por la expresión general (para un nudo genérico n)

$\mathbf {{\mbox{R}}_{\mbox{n}}{\mbox{ }}{\mbox{d}}_{\mbox{n}}^{'}={\mbox{0}}}$

(19)

siendo $\mathbf {{\mbox{d}}_{\mbox{n}}^{'}}$ el vector de desplazamientos (multiplicados por E) del nudo restringido n, y R_n una de las tres unidades matriciales operativas siguientes, según el tipo de aparato de apoyo en contacto con él

$\mathbf {{\mbox{R}}_{\mbox{n}}} =\left[{\begin{array}{cc}1&0\\0&1\end{array}}\right];{\mbox{ }}\mathbf {{\mbox{R}}_{\mbox{n}}} =\left[{\begin{array}{cc}0&1\end{array}}\right];{\mbox{ }}\mathbf {{\mbox{R}}_{\mbox{n}}} =\left[{\begin{array}{cc}1&0\end{array}}\right]$

(20)

donde el subíndice designará el nudo restringido, correspondiendo la primera de ellas a un apoyo fijo (nulidad de desplazamiento horizontal u_n y vertical v_n), la segunda a un apoyo móvil con plano de rodamiento horizontal (nulidad de desplazamiento vertical v_n) y la tercera a un apoyo móvil con rodamiento vertical (desplazamiento horizontal u_n nulo). El ensamblaje de las matrices Elementales de Restricción Cinemática dadas en (20) proporcionará una matriz R, Matriz de Restricciones Cinemáticas, que genera un primer grupo de restricciones en las incógnitas $\mathbf {\mbox{d'}}$ que, junto con las tres anteriores y aplicando el método de los Multiplicadores de Lagrange, forman el sistema matricial

$\left[{\begin{array}{ccc|c}\mathbf {0} &\mathbf {H} &\mathbf {0} &\mathbf {R} ^{\mathbf {T} }\\\mathbf {B} &\mathbf {0} &\mathbf {-I} &\mathbf {0} \\\mathbf {0} &\mathbf {-I} &\mathbf {K} &\mathbf {0} \\\hline \mathbf {R} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {\mbox{σ}} \\\mathbf {\mbox{ε'}} \\\hline {\mbox{ }}{\mbox{ }}{\mathbf {\mbox{λ}} }_{\mathbf {R} }\end{array}}\right]=\left[{\begin{array}{c}\mathbf {v} _{\mathbf {H} }\\\mathbf {0} \\\mathbf {0} \\\hline \mathbf {0} \end{array}}\right]$

(21)

que ahora ya es compatible y determinado, siendo respectivamente $\mathbf {\mbox{d'}}$ , σ, $\mathbf {\mbox{ε'}}$ y λ_R los vectores columna de desplazamientos de nudos (multiplicados por E), tensiones en los nodos, deformaciones en éstos (multiplicadas por E) y Multiplicadores de Lagrange asociados a las restricciones de los apoyos (reacciones con sentido contrario), mientras que v_H es el vector asociado a la matriz H, esto es, el vector de acciones que actúan en todos los nudos de la estructura, incluidos los apoyos. Como puede observarse, la matriz de coeficientes del sistema (21) no puede ser simétrica, como así ocurre en el Método de la Rigidez para estructuras de barras, debido al incumplimiento de la relación de Contragradiencia; incumplimiento que es consecuencia, por una parte, del carácter aproximado de las funciones de interpolación con las que se ha generado la matriz de Compatibilidad B y, por otra, a la naturaleza también aproximada de la matriz de Equilibrio H, la cual se ha obtenido a partir de la ya comentada hipótesis de variación lineal de las tensiones a lo largo de los lados de los elementos que modelizan la estructura.

3.5. Matrices Elementales de Restricciones Estáticas Internas

El sistema (21) tiene, pues, solución, mas en general ésta será inadmisible, pues no está asegurado aún que se cumpla el equilibrio estático de cada uno de los elementos que discretizan a la estructura. Pero, como ya se señaló anteriormente, el procedimiento que aquí se presenta permite exigir que se cumpla tal condición, es decir, obligar a que las fuerzas nodales que actúan en cada elemento estén en estricto equilibrio; por lo tanto, se ha de introducir para todo elemento aislado e, por una parte la ecuación matricial de equilibrio de Fuerzas Nodales Equivalentes:

$\mathbf {{\mbox{F}}_{\mbox{i}}^{\mbox{e}}} +\mathbf {{\mbox{F}}_{\mbox{j}}^{\mbox{e}}} +$$\mathbf {{\mbox{F}}_{\mbox{k}}^{\mbox{e}}} +\mathbf {{\mbox{F}}_{\mbox{m}}^{\mbox{e}}} =$$(\mathbf {{\mbox{I}}^{\mbox{e}}} +\mathbf {{\mbox{J}}^{\mbox{e}}} +\mathbf {{\mbox{K}}^{\mbox{e}}} +$$\mathbf {{\mbox{M}}^{\mbox{e}}} ){\mathbf {\mbox{σ}} ^{\mbox{e}}}={\mbox{0}}$

(22)

y por otra el equilibrio de momentos de dichas fuerzas. Tomando, pues, momentos respecto al centro del rectángulo (o cualquier otro punto) se tiene que cumplir la ecuación matricial, de acuerdo con la Figura 1 (c)

$\left[{\begin{array}{cc}b&-a\end{array}}\right]\mathbf {{\mbox{F}}_{\mbox{i}}^{\mbox{e}}} +\left[{\begin{array}{cc}b&a\end{array}}\right]\mathbf {{\mbox{F}}_{\mbox{j}}^{\mbox{e}}} +\left[{\begin{array}{cc}-b&a\end{array}}\right]\mathbf {{\mbox{F}}_{\mbox{k}}^{\mbox{e}}} +\left[{\begin{array}{cc}-b&-a\end{array}}\right]\mathbf {{\mbox{F}}_{\mbox{m}}^{\mbox{e}}} ={\mbox{0}}$

(23)

Sustituyendo en las ecs.(22) y (23) las expresiones dadas en (5)-(8), operando y simplificando, se obtienen finalmente las tres ecuaciones homogéneas de equilibrio para cada elemento e aislado

$\left[{\begin{array}{ccc|ccc|ccc|ccc}-b&0&-a&b&0&-a&b&0&a&-b&0&a\\0&-a&-b&0&-a&b&0&a&b&0&a&-b\\-b^{2}&a^{2}&0&b^{2}&-a^{2}&0&-b^{2}&a^{2}&0&b^{2}&-a^{2}&0\end{array}}\right]\mathbf {{\mbox{σ}}^{\mbox{e}}} =\left[{\begin{array}{c}0\\0\\0\end{array}}\right]$

(24)

y signando por T^e a la matriz de coeficientes, ésta vendrá dada, pues, por

$\mathbf {{\mbox{T}}^{\mbox{e}}} =\left[{\begin{array}{ccc|ccc|ccc|ccc}-b&0&-a&b&0&-a&b&0&a&-b&0&a\\0&-a&-b&0&-a&b&0&a&b&0&a&-b\\-b^{2}&a^{2}&0&b^{2}&-a^{2}&0&-b^{2}&a^{2}&0&b^{2}&-a^{2}&0\end{array}}\right]$

(25)

que no es más que la unidad matricial operativa que llamaremos matriz Elemental de Restricciones Estáticas Internas, con cuyo ensamblaje se obtendrá otra matriz T, matriz de Restricciones Estáticas Internas para toda la estructura y, por tanto, un segundo grupo de restricciones en las incógnitas de tensión σ que obligará a que todos y cada uno de los elementos que la conforman estén en equilibrio estático, extremo éste al que ha de renunciar el tradicional MEF. Añadiendo, pues, esta última matriz al sistema (21) se obtiene

$\left[{\begin{array}{ccc|cc}\mathbf {0} &\mathbf {H} &\mathbf {0} &{\mbox{ }}{\mbox{ }}\mathbf {R} ^{\mathbf {T} }&\mathbf {0} \\\mathbf {B} &\mathbf {0} &\mathbf {-I} &\mathbf {0} &\mathbf {T} ^{\mathbf {T} }\\\mathbf {0} &\mathbf {-I} &\mathbf {K} &\mathbf {0} &\mathbf {0} \\\hline \mathbf {R} &\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &\mathbf {T} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {\mbox{σ}} \\\mathbf {\mbox{ε'}} \\\mathbf {\mbox{λ}} _{\mathbf {R} }\\\mathbf {\mbox{λ}} _{\mathbf {T} }\end{array}}\right]=\left[{\begin{array}{c}\mathbf {v} _{\mathbf {H} }\\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \end{array}}\right]$

(26)

sistema compatible y determinado, aunque, como se ha dicho, no simétrico, debido a que H y B no son recíprocamente transpuestas una de la otra, y en donde λ_T es el vector que contiene los valores -con signo contrario- de la deformación que ha de sufrir cada elemento para que se cumplan las condiciones de equilibrio dadas en la ec.(24). Obsérvese cómo, a consecuencia de las condiciones introducidas, las relaciones de compatibilidad (18) -segundo bloque de (21)- se ven modificadas, tal y como expresa el segundo bloque de (26), lo que significa un nuevo apartamiento de este procedimiento respecto al tradicional.

3.6. Matrices Elementales de Restricciones Estáticas Externas o Perimetrales

La resolución del sistema (26) devolverá unos valores de las incógnitas que, si bien mejoran considerablemente a los que se obtienen por el MEF tradicional -lo que puede comprobarse resolviendo este sistema para los ejemplos que figuran al final de este trabajo- no deben darse como definitivos puesto que, de forma similar a las anteriores restricciones, aún falta por introducir un tercer y último grupo de ellas que describan el estado tensional del contorno del dominio, estado que puede ser de cualquier tipo pero que por brevedad aquí supondremos nulo, como suele ser habitual. La restricción estática de tensión normal nula en los lados i-j ( ${\sigma }_{y}^{i}={\sigma }_{y}^{j}=0$ ) y k-m ( ${\sigma }_{y}^{k}={\sigma }_{y}^{m}=0$ ) de un elemento e perteneciente al perímetro vendrá dada, respectivamente, por las ecuaciones homogéneas

$\mathbf {{\mbox{P}}_{\mbox{i-j}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =$$\mathbf {0} {\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}y{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{P}}_{\mbox{k-m}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =$$\mathbf {0}$

(27)

mientras que la condición de tensión normal nula en lo lados i-m ( ${\sigma }_{x}^{i}={\sigma }_{x}^{m}=0$ ) y j-k ( ${\sigma }_{x}^{j}={\sigma }_{x}^{k}=0$ ) será, de forma similar,

$\mathbf {{\mbox{P}}_{\mbox{i-m}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =$$\mathbf {0} {\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}y{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {{\mbox{P}}_{\mbox{j-k}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =$$\mathbf {0}$

(28)

donde cada una de las matrices de (27) y (28) vienen dadas por

${\begin{array}{c}\mathbf {{\mbox{P}}_{\mbox{i-j}}^{\mbox{e}}} =\left[{\begin{array}{ccc|ccc|ccc|ccc}0&1&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&0&0&0&0\end{array}}\right]\\\mathbf {{\mbox{P}}_{\mbox{k-m}}^{\mbox{e}}} =\left[{\begin{array}{ccc|ccc|ccc|ccc}0&0&0&0&0&0&0&1&0&0&0&0\\0&0&0&0&0&0&0&0&0&0&1&0\end{array}}\right]\\\mathbf {{\mbox{P}}_{\mbox{i-m}}^{\mbox{e}}} =\left[{\begin{array}{ccc|ccc|ccc|ccc}1&0&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&1&0&0\end{array}}\right]\\\mathbf {{\mbox{P}}_{\mbox{j-k}}^{\mbox{e}}} =\left[{\begin{array}{ccc|ccc|ccc|ccc}0&0&0&1&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&1&0&0&0&0&0\end{array}}\right]\end{array}}$

(29)

que no son más que las unidades operativas con cuyo ensamblaje se obtendrá una matriz P_, matriz de Restricciones Estáticas Perimetrales y, por tanto, un tercer y último grupo de restricciones en las incógnitas de tensión σ que describirá el estado tensional del contorno del dominio en estudio -extremo éste al que también ha de renunciar el MEF tradicional- con lo que el sistema contará ya con la información suficiente para que su resolución pueda devolver unos resultados aceptablemente satisfactorios, sobre todo si la discretización es lo suficientemente tupida.

4. Ensamblaje de las matrices anteriores y del vector de términos independientes

Como ya se indicó, el método que aquí se propone trata de conseguir mediante una exclusiva manipulación matricial el sistema de ecuaciones de equilibro, compatibilidad y constitutivas, así como unas determinadas condiciones de restricción estática y cinemática, cuya resolución bajo el Método de los Multiplicadores de Lagrange devolverá directamente los valores de todas las incógnitas que describen el comportamiento mecánico de la estructural en estudio. Dicho sistema se conseguirá a partir de los correspondientes ensamblajes de las matrices elementales antes deducidas y de los vectores de términos independientes, cada uno de tales ensamblajes se detalla en los siguientes apartados.

4.1. Ensamblaje de la Matriz de Equilibrio (H)

Una vez signados los nodos de cada uno de los elementos por i, j, k y m, de acuerdo con la Fig.1 y partiendo de un arreglo de tantas filas como nudos n y tantas columnas como elementos e tenga la estructura, el ensamblaje de las matrices elementales dadas en (5)-(8) se podrá abordar por filas o por columnas. En el primer caso se ubicará en la fila f asociada al nudo f las submatrices I^e, J^e, K^e o M^e, correspondientes a los elementos e cuyos vértices incidan en dicho nudo, colocando cada una de ellas en la columna asociada al elemento respectivo. De forma recíproca, se si procede al ensamblaje por columnas, en la columna c, asociada al elemento c, se ubicarán las cuatro submatrices (5)-(8) de dicho elemento en la correspondiente fila (nudo) donde inciden cada uno de los cuatro nodos del mismo. Los ejemplos que más adelante se presentan aclararán mejor dicho ensamblaje.

4.2. Ensamblaje de la Matriz de Compatibilidad (B)

En las estructuras de barras, tanto de nudos rígidos como articulados, es sabido que se cumple la Relación de Contragradiencia, por lo que en estos casos la Matriz de Compatibilidad vendrá dada directamente por H^T, sin necesidad, pues, de proceder a un nuevo ensamblaje para la obtención de B (o, recíprocamente, ensamblar B y obtener H por transposición de B). Pero dado el carácter sólo aproximado de las funciones de interpolación que sustentan, como ya se ha dicho, al procedimiento que nos ocupa, no es posible determinar la Matriz de Compatibilidad a partir de la ya calculada matriz H. Será necesario, pues, proceder al pertinente ensamblaje de las matrices (12), ensamblaje recíproco del anterior. No obstante, la terminología utilizada en este trabajo (signando con mayúsculas a las matrices elementales de Fuerzas Nodales, expresiones (5)-(8), y con letras minúsculas a las correspondientes matrices de compatibilidad ) permite evitar dicho ensamblaje, pues bastará con transponer la matriz de equilibrio H ya ensamblada (y antes de que tome sus valores numéricos) y, a continuación, sencillamente cambiar la mayúsculas por las minúsculas para obtener la matriz B, como puede observarse en el ejemplo 1.

4.3. Ensamblaje de la Matriz Constitutiva (K)

El ensamblaje de esta matriz es inmediato, pues no es más que un arreglo diagonal por bloques, de dimensiones e×e, donde e es el número de elementos, ubicando la expresión (17) en dicha diagonal.

4.4. Ensamblaje de la Matriz de Restricciones Cinemáticas (R)

Partiendo de un arreglo de tantas filas como aparatos de apoyo tenga la estructura y tantas columnas como nudos totales n tenga ésta, se ubicará en cada una de las filas una de las matrices (20) (según sea el tipo de apoyo con el que se modele la sustentación) y en la columna correspondiente al nudo restringido.

4.5. Ensamblaje de la Matriz de Restricciones Estáticas Internas (T)

De la misma forma que la Matriz Constitutiva, esta otra también se obtiene a partir del mismo arreglo diagonal por bloques, en cuya diagonal principal se ha de colocar ahora la matriz elemental dada en (25).

4.6. Ensamblaje de la Matriz de Restricciones Estáticas Perimetrales (P)

Partiendo de un arreglo con tantas filas como lados a restringir y tantas columnas como elementos totales tenga la estructura, se ubicarán en cada fila una de las matrices dadas en (29), según sea el lado a restringir, y en la columna correspondiente al elemento al que pertenece dicho lado.

4.7. Ensamblaje del vector de términos independientes (v)

Según todo lo anterior, el vector de términos independientes será un vector cuyas componentes por bloques son todas nulas, excepto el término asociado a las ecuaciones de equilibrio, que contendrá las componentes de acción v_H que actúan en los nudos de la estructura.

5. Sistema Completo

Ordenando apropiadamente todas las matrices y vectores anteriores, se llega finalmente al Sistema Completo siguiente

${\begin{array}{ccc}{\mbox{ }}&{\begin{array}{ccccccc}{\mbox{ }}&{\underset {\begin{array}{c}\downarrow \\(2n)\end{array}}{\mathbf {d'} }}&{\underset {\begin{array}{c}\downarrow \\(12e)\end{array}}{{\mathbf {\mbox{σ}} }_{}}}&{\underset {\begin{array}{c}\downarrow \\(12e)\end{array}}{\mathbf {\mbox{ε}} \mathbf {'} _{}}}&{\underset {\begin{array}{c}\downarrow \\(?)\end{array}}{{\mathbf {\mbox{λ}} }_{\mathbf {R} }}}&{\underset {\begin{array}{c}\downarrow \\(3e)\end{array}}{{\mathbf {\mbox{λ}} }_{\mathbf {T} }}}&{\underset {\begin{array}{c}\downarrow \\(?)\end{array}}{{\mathbf {\mbox{λ}} }_{\mathbf {P} }}}&{\mbox{ }}\end{array}}&{\mbox{ }}\\{\begin{array}{cc}{\mbox{Ecuaciones de equilibrio (2n)}}&\rightarrow \\\qquad \qquad \quad {\mbox{''}}\qquad {\mbox{de compatibilidad (12e)}}&\rightarrow \\\quad \qquad {\mbox{''}}\qquad {\mbox{constitutivas (12e)}}&\rightarrow \\{\mbox{Restricciones cinemáticas (?)}}&\rightarrow \\\quad \qquad {\mbox{''}}\qquad {\mbox{estáticas internas (3e)}}&\rightarrow \\\quad \qquad {\mbox{''}}\qquad {\mbox{estáticas externas (?)}}&\rightarrow \end{array}}&{\left[{\begin{array}{ccc|ccc}\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {H} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &{\mbox{ }}{\mbox{ }}\mathbf {R} ^{\mathbf {T} }&\mathbf {0} &\mathbf {0} \\\mathbf {B} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {-I} &\mathbf {0} &\mathbf {T} ^{\mathbf {T} }&\mathbf {P} ^{\mathbf {T} }\\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {-I} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {K} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\hline \mathbf {R} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {T} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {P} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]}&{\!\!\!\!\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {\mbox{σ}} \\\mathbf {\mbox{ε'}} \\\hline \mathbf {\mbox{λ}} _{\mathbf {R} }\\\mathbf {\mbox{λ}} _{\mathbf {T} }\\\mathbf {\mbox{λ}} _{\mathbf {P} }\end{array}}\right]=\left[{\begin{array}{c}\mathbf {v} _{\mathbf {H} }\\\mathbf {0} \\\mathbf {0} \\\hline \mathbf {0} \\\mathbf {0} \\\mathbf {0} \end{array}}\right]{\acute {o}}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {M} {\mbox{ }}\mathbf {s=v} }\end{array}}$

(30)

siendo s el vector solución y en donde, para facilitar su lectura, se han colocado a su izquierda la denominación de cada uno de los bloques de ecuaciones con sus respectivas dimensiones y en la parte superior las incógnitas y dimensiones asociadas a cada columna. De nuevo las relaciones de compatibilidad (segundo bloque del anterior sistema) sufren una última modificación respecto al sistema (26) debido a las condiciones estáticas añadidas, las cuales originan un último vector incógnita λ_P que contiene los valores -con signo contrario- de la nueva deformación que ha de sufrir cada elemento para que se cumplan las condiciones estáticas perimetrales dadas en (27) y (28).

6. Sistema Reducido

En la publicación “Las tecnologías cuánticas” (Eduardo Arroyo, de National Geographic de 2017) se preconiza que “El desarrollo de la computación cuántica anuncia una auténtica revolución en el campo de la tecnología”, por lo que, por compleja que sea la estructura, la capacidad de los futuros ordenadores puede ser presumiblemente suficiente para la resolución del Sistema Completo . No obstante, éste puede reducirse si se desea eliminado las deformaciones, ya que éstas no suelen ser magnitudes fundamentales en el análisis estructural (eliminación que en absoluto ha de significar la renuncia a su cálculo, pues éste puede hacerse a posteriori). En efecto, despejando el vector de deformaciones $\mathbf {\mbox{ε'}}$ del segundo bloque de (30) y sustituyéndolo en el tercero, se obtiene la ecuación matricial

$\mathbf {K} {\mbox{ }}\mathbf {B} {\mbox{ }}\mathbf {d'-I} {\mbox{ }}\mathbf {\mbox{σ+K}} {\mbox{ }}\mathbf {T} ^{\mathbf {T} }{\mathbf {\mbox{λ}} }_{\mathbf {T} }\mathbf {+K} {\mbox{ }}\mathbf {P} ^{\mathbf {T} }{\mathbf {\mbox{λ}} }_{\mathbf {P} }\mathbf {=0}$

(31)

por lo que el Sistema Completo quedará modificado y en la forma siguiente

$\left[{\begin{array}{cc|ccc}\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {H} &{\mbox{ }}{\mbox{ }}\mathbf {R} ^{\mathbf {T} }&\mathbf {0} &{\mbox{0}}\\\mathbf {K} {\mbox{ }}\mathbf {B} &{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {-I} &\mathbf {0} &{\mbox{ }}\mathbf {K} {\mbox{ }}\mathbf {T} ^{\mathbf {T} }&{\mbox{ }}\mathbf {K} {\mbox{ }}\mathbf {P} ^{\mathbf {T} }\\\hline \mathbf {R} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {T} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {P} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {\mbox{σ}} \\\mathbf {\mbox{λ}} _{\mathbf {R} }\\\mathbf {\mbox{λ}} _{\mathbf {T} }\\\mathbf {\mbox{λ}} _{\mathbf {P} }\end{array}}\right]=\left[{\begin{array}{c}\mathbf {v} _{\mathbf {H} }\\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \end{array}}\right]$

(32)

que es el Sistema Reducido, habiéndose eliminado, pues, 12e incógnitas y otras tantas ecuaciones, lo cual supone una gran disminución de las dimensiones del sistema, puesto que el número de elementos e generalmente será en la práctica muy elevado. El anterior sistema puede adoptar una forma más simple pre-multiplicando el segundo bloque por la inversa de K, la cual es una matriz simétrica y extraordinariamente simple, observando la (16) y (17), con lo que se obtiene finalmente

$\left[{\begin{array}{cc|ccc}\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {H} &{\mbox{ }}{\mbox{ }}\mathbf {R} ^{\mathbf {T} }&\mathbf {0} &{\mbox{0}}\\\mathbf {B} &{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {-K} ^{-1}&\mathbf {0} &{\mbox{ }}\mathbf {T} ^{\mathbf {T} }&{\mbox{ }}\mathbf {P} ^{\mathbf {T} }\\\hline \mathbf {R} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {0} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {T} &\mathbf {0} &\mathbf {0} &\mathbf {0} \\\mathbf {0} &{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}\mathbf {P} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {\mbox{σ}} \\\mathbf {\mbox{λ}} _{\mathbf {R} }\\\mathbf {\mbox{λ}} _{\mathbf {T} }\\\mathbf {\mbox{λ}} _{\mathbf {P} }\end{array}}\right]=\left[{\begin{array}{c}\mathbf {v} _{\mathbf {H} }\\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \end{array}}\right]$

(33)

De forma semejante se pueden eliminar también las tensiones, pues despejando este vector en (31)

$\mathbf {{\mbox{σ=K}}{\mbox{ }}{\mbox{B}}{\mbox{ }}{\mbox{d'+K}}{\mbox{ }}{\mbox{T}}^{\mbox{T}}{\mbox{λ}}_{\mbox{T}}{\mbox{+K}}{\mbox{ }}{\mbox{P}}^{\mbox{T}}{\mbox{λ}}_{\mbox{P}}}$

(34)

y sustituyéndolo en (33) se obtiene el sistema

$\left[{\begin{array}{cccc}\mathbf {HKB} &\mathbf {HKT} ^{\mathbf {T} }&\mathbf {HKP} ^{\mathbf {T} }&\mathbf {R} ^{\mathbf {T} }\\\mathbf {TKB} &\mathbf {TKT} ^{\mathbf {T} }&\mathbf {TKP} ^{\mathbf {T} }&\mathbf {0} \\\mathbf {PKB} &\mathbf {PKT} ^{\mathbf {T} }&\mathbf {PKP} ^{\mathbf {T} }&\mathbf {0} \\\mathbf {R} &\mathbf {0} &\mathbf {0} &\mathbf {0} \end{array}}\right]\left[{\begin{array}{c}\mathbf {d'} \\\mathbf {{\mbox{λ}}_{\mbox{T}}} \\\mathbf {{\mbox{λ}}_{\mbox{P}}} \\\mathbf {{\mbox{λ}}_{\mbox{R}}} \end{array}}\right]=\left[{\begin{array}{c}\mathbf {{\mbox{v}}_{\mbox{H}}} \\\mathbf {0} \\\mathbf {0} \\\mathbf {0} \end{array}}\right]$

(35)

en las incógnitas exclusivas de desplazamiento $\mathbf {d'}$ y de Multiplicadores de Lagrange. Pero como quiera que los multiplicadores asociados a las restricciones estáticas no tendrán, en general, interés práctico alguno -serán las deformaciones, con sentido contrario, que hay que imponer a los elementos para que se cumplan tales condiciones- también convendrá eliminarlas, lo que podrá hacerse con ayuda de unas determinadas matrices de Transformación utilizando el procedimiento que se expone en la Ref. [3], Cap.8 y que aquí no se incluyen por motivos de espacio y por no ser objeto primordial de esta investigación. Asimismo será fácil la eliminación de $\mathbf {{\mbox{λ}}_{\mbox{R}}} ;$ e incluso de aquellos elementos de $\mathbf {d'}$ que se sabe que son nulos, cuales son los desplazamientos en los aparatos de apoyo, consiguiéndose entonces un sistema de idénticas dimensiones al correspondiente del Método Directo de la Rigidez, pero con la salvedad de que se obtendrán unos resultados muy mejorados respecto a los que se obtendrían por el MEF tradicional, como se demostrará con los ejemplos que a continuación se detallan, pues como ya se ha indicado, el MEF tradicional parte de un sistema con mucha menos información que el método que aquí se propone. No obstante, aunque el sistema ha de constar de 2n+3e incógnitas, más las asociadas a los condiciones de contorno (cinemáticas y estáticas) hay que tener en cuenta que en la práctica será necesario modelizar el dominio en estudio con un elevado número de elementos para conseguir una solución aceptable y, por lo tanto, el número de condiciones de contorno será muy pequeño frente a 2n+3e, por lo que no merecerá la pena reducir el sistema mediante tales eliminaciones. Por otra parte, si bien el sistema (35)-y los derivados de éste si se eliminan algunos o todos los Multiplicadores de Lagrange- tiene muchas menos incógnitas que los anteriores, será obligado proceder al cálculo posterior de deformaciones y tensiones, post-proceso que es imprescindible en el tradicional MEF, pero innecesario en el método que aquí se propone utilizando cualquiera de los sistemas (30), (32) o (33), los cuales devuelven de forma directa los valores de todas la incógnitas.

7. Ejemplo práctico 1: voladizo en flexión pura y en flexión simple

La estructura que se muestra en la Figura 2 es el conocido “Patch test” aplicado a una viga en voladizo y analizada en la Ref. [2]^§4.14.1 con los valores E=1500, ν = 0.25 y t = 1 para dos tipos de cargas: a) Flexión pura, producida por dos fuerzas horizontales iguales y contrarias aplicadas en los vértices de la cara derecha y b) Flexión simple por aplicación de sendas cargas verticales de igual magnitud y dirección aplicadas también en dichos puntos, tal y como se muestra en la figura. Aunque algunos autores prefieren suponer uno de los apoyos como móvil, parece ser lo más coherente para modelizar el empotramiento considerar fijos a ambos y en base a esto se analizará la estructura con la metodología que aquí se propone.

Figura 2. Viga en voladizo Dimensiones, numeración de nudos, de elementos y cargas en el extremo libre

Se pretende calcular la flecha en el punto A y la tensión normal σ_x en B utilizando una discretización de 5 elementos rectangulares de 4 nodos (por lo que los valores que se obtendrán por el método tradicional diferirán de los que figuran en la referencia, donde los elementos tienen una geometría irregular). A continuación se expone el desarrollo completo del cálculo con el programa Mathematica 10.0 para las cargas a), evitando por comodidad de escritura subíndices y superíndices en los Input a introducir, habiéndose sombreado éstos para facilitar su lectura.

7.1 Introducción de datos generales

De acuerdo con la Figura 2 y la expresión :

7.2 Ensamblaje de las matrices de Equilibrio H y de Compatibilidad B

Puesto que todos los elementos son idénticos, llamando I⁰, J⁰, etc. a las matrices dadas en (5)-(8), y de los apartados §4.1. y §4.2.:

y cambiando las mayúsculas por minúsculas en la matriz transpuesta que se obtiene con la anterior orden, el esamblaje de la matriz B viene dado de forma inmediata por:

donde es de señalar cómo estas matrices presentan una gran simplicidad y unas densidades muy inferiores a las que se obtienen al ensamblar la Matriz de Rigidez en el método tradicional.

7.3 Matrices elementales de Fuerzas Nodales y de Compatibilidad

De las (5)-(8) y (12):

7.4 Valores numéricos de las Matrices de Equilibrio H y de Compatibilidad B de la Estructura

Una vez leídas las anterores matrices, la expresión numérica de las matrices de Equilibrio y Compatibilidad se obtendrán mediante las órdenes:

7.5 Valor numérico de la Matriz Constitutiva K de la Estructura

De la (16), (17) y del §4.3:

7.6 Valor numérico de la Matriz de Restricciones Cinemáticas R

De la primera de (20) y del §4.4.:

7.7 Valor numérico de la Matriz de Restricciones Estáticas Internas T

De la (25) y del §4.5.:

7.8 Valor numérico de la Matriz de Restricciones Estáticas Perimetrales P

De (29) y del §4.6., y teniendo en cuenta que los lados a restringir son los bordes superior e inferior de la barra, esto es, 10 lados:

7.9 Vector de términos independientes v

Llamando v_H, v_B, v_K, v_R, v_T y v_P a los vectores asociados a las respectivas matrices H, B, K, R, T y P e I_d a la matriz identidad, se habrá de escribir, de acuerdo con §4.7.:

7.10 Sistema Completo y cálculo del vector solución s

Aunque bastaría con tomar el sistema (32), el (33) o el (35) -que tendría unas dimensiones de 2n más multiplicadores de Lagrange- aquí se resolverá el (30), de mayores dimensiones, para que el lector pueda extraer de dicha resolución los valores de todas las incógnitas involucradas en el cálculo de la estructura en estudio, lo que se consigue mediante el Input

con lo que se obtiene el vector solución s, del que extraeremos aquéllos valores que se requerían en el enunciado y que se muestran en la Tabla 1 en la que se incluyen tambien los resultados para las cargas b) -repitiendo todos los cálculos anteriores mediante la modificación del input del §7.9- así como los valores teóricamente exactos y los errores en %,

Tabla 1. Resultados del Ejemplo 1.
Método	Caso a): flexión pura				Caso b): flexión simple
	vA		σxB		vA		σxB
	Valor	Error %	Valor	Error %	Valor	Error %	Valor	Error %
MEF (Integ. analítica)	68.1818	31.82	-2181.82	27.27	70	31.79	-2945.45	27.27
Propuesto	100	0	-3000	0	101.5	0.011	-4050	0
Valor exacto	100	-	-3000	-	102.625	-	-4050	-

Tabla en la que es de resaltar cómo, para la metodología propuesta, todos los valores son exactos, excepto el correspondiente a la flecha para la flexión simple, caso que se estudiará con más detalle en el apartado siguiente.

Para la flexión pura, los desplazamientos de todos los nudos y las tensiones en los nodos de los cinco elementos se exponen, respectivamente, en las Tablas 2 y 3.

Tabla 2. Desplazamientos de los nudos de la estructura para la flexión pura.
u1= u7	v1= v7	u2= - u8	v2= v8	u3= - u9	v3= v9	u4= - u10	v4= v10	u5= - u11	v5= v11	u6= - u12	v6= v6
0	0	4	4	8	16	12	36	16	64	20	100

Tabla 3. Tensiones en los nodos de todos los elementos para la flexión pura.
Elemento	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{i}}}$	${\textstyle {\mbox{τ}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{j}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{j}}}$	${\textstyle {\mbox{τ}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{k}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{k}}}$	${\textstyle {\mbox{τ}}^{\mbox{k}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{m}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{m}}}$	${\textstyle {\mbox{τ}}^{\mbox{m}}}$
1, 2, 3, 4 y 5	3000	0	0	3000	0	0	3000	0	0	3000	0	0

7.11 Estudio de la convergencia

Puesto que, como se ha señalado, los resultados de todas las incógnitas son exactos para la flexión pura, aún para una modelización tan grosera, no cabe analizar su convergencia, por lo que se repetirán los cálculos sólo para la flexión simple con modelizaciones de 10, 15 y 20 elementos, lo que puede comprobarse a la vista de la Tabla 4 en la que se repiten los valores para 5 elementos y en la que claramente se observa cómo el valor del desplazamiento vertical del punto A tiende paulatinamente a su valor exacto.

Tabla 4. Análisis de convergencia del Ejemplo1 en flexión simple.
Nº de elementos	Flecha en A	Tensión en B
5	101.5	-4050
10	102.25	-4050
15	102.389	-4050
20	102.437	-4050
Valor exacto	102.625	-4050

Por motivos de espacio sólo se detallan a continuación los resultados correspondientes a la modelización de 10 elementos, figurando en la Tabla 5 los desplazamientos de los nudos de la cara inferior, cuyos desplazamientos v son idénticos a los correspondientes de la cara superior, mientras que los u serán de signo contrario.

Tabla 5. Desplazamientos de los nudos de la estructura para la flexión simple.
u1	v1	u2	v2	u3	v3	u4	v4	u5	v5	u6	v6	u7	v7	u8	v8	u9	v9	u10	v10	u11	v11
0	0	2.85	1.625	5.4	6.05	7.65	12.825	9.6	21.7	11.25	32.375	12.6	44.55	13.65	57.925	14.4	72.2	14.85	87.075	15	102.25

En la Tabla 6 se muestran los valores de las tensiones -que normalmente son las incógnitas de más interés para el calculista- en los nodos de los diez elementos.

Tabla 6. Tensiones en los nodos de los elementos para flexión simple y 10 elementos.
Elemento	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{i}}}$	${\textstyle {\mbox{τ}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{j}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{j}}}$	${\textstyle {\mbox{τ}}^{\mbox{i}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{k}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{k}}}$	${\textstyle {\mbox{τ}}^{\mbox{k}}}$	${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{m}}}$	${\textstyle {\mbox{σ}}_{\mbox{y}}^{\mbox{m}}}$	${\textstyle {\mbox{τ}}^{\mbox{m}}}$
1	4275	0	150	4275	0	150	-4275	0	150	-4275	0	150
2	3825	0	150	3825	0	150	-3825	0	150	-3825	0	150
3	3375	0	150	3375	0	150	-3375	0	150	-3375	0	150
4	2925	0	150	2925	0	150	-2925	0	150	-2925	0	150
5	2475	0	150	2475	0	150	-2475	0	150	-2475	0	150
6	2025	0	150	2025	0	150	-2025	0	150	-2025	0	150
7	1575	0	150	1575	0	150	-1575	0	150	-1575	0	150
8	1125	0	150	1125	0	150	-1125	0	150	-1125	0	150
9	675	0	150	675	0	150	-675	0	150	-675	0	150
10	225	0	150	225	0	150	-225	0	150	-225	0	150

A la vista de dicha tabla se han de hacer las siguientes puntualizaciones:

- Siguiendo el habitual criterio de adoptar para la tensión en un nudo la media de las tensiones de los nodos que él inciden, puede comprobarse cómo las tensiones en cada nudo coinciden con la solución exacta, excepto en los nudos extremos. Así pues, la tensión en la dirección horizontal en el nudo B se obtiene a partir de la correspondiente tensión en el nodo k del elemento 1 y la del nodo m del elemento 2, en negritas en la tabla. Con estos valores la tensión en este punto es de -4050, que es el valor exacto.

- Todas las tensiones en la dirección vertical son nulas, lo que se ha conseguido al imponer las condiciones estáticas perimetrales introducidas en el §7.8, extremo éste al que ha de renunciar el procedimiento tradicional del MEF y asumir necesariamente algún grado de imprecisión en los resultados. Concretamente, para este ejemplo puede comprobarse que dichas tensiones tienen un valor de 545.455.

- La tensión tangencial es constante e igual a 150 en cualquier sección transversal del voladizo, distribución que se aparta de la real, pues ésta ha de ser parabólica y con valor nulo en los nodos; no obstante se cumple que la resultante de las tensiones en ambas distribuciones son idénticas, siendo su valor 300, que no es más que el cortante originado por las cargas aplicadas. Con el MEF habitual el valor de esta tensión en los nodos vale 818.82, por lo que su resultante sería muy superior al valor exacto del cortante.

8. Ejemplo práctico 2: voladizo en flexión simple

La estructura que se muestra en la Figura 3 está analizada las Refs. 4^§6.11, 5^§9.8 y 6^§4.5 con los datos de la Tabla 7 para cada una de ellas y en la que se pretende calcular la flecha en el punto C y la tensión normal σ_x en los puntos A y B utilizando una discretización de sólo 5 elementos rectangulares de 4 nodos

Figura 3. Viga en voladizo Dimensiones, numeración de nudos, de elementos y cargas en el extremo libre

Tabla 7. Datos para cada Ref.
Autor	h	E	𝜈	t	P
Ref. 4	2	104	0.3	1	10
Ref. 5*	2	1500	0.25	1	150
Ref. 6	1	2×108	0.2	0.1	450

^*Aunque estos datos no aparecen de forma explícita en la referencia, los resultados que figuran en ella concuerdan con los valores de la Tabla.

Procediendo de idéntica forma que en el ejemplo anterior, se obtienen los resultados que figuran en la Tabla 8 en la que, como puede observarse, los resultados dados por el MEF tradicional son altamente inexactos debido, entre otras cosas, a la grosera modelización adoptada, sobre todo para los datos de la Ref. [6] en la que los errores se disparan debido a un valor mayor de la relació a/b del elemento rectangular elegido, relación que no afecta negativamente a los valores obtenidos por el método propuesto, sino más bien al contrario, como puede observarse comparando los valores de la última fila de las Tablas 8 y 9. Para solventar este problema, Babuška y Rheinboldt en 1970, y más recientemente en las referencias [5], [6] y [8], se proponen procedimientos alternativos que, aunque con un mayor número de elementos, consiguen unos resultados que pueden considerarse totalmente aceptables. Concretamente, en la Ref. [6] figuran unos resultados muy similares a los que se obtienen con el procedimiento que en este trabajo se propone, aunque en éste el número de elementos es de 5, mientras que en aquél son necesarios al menos 1192. Por último, se resumen en la Tabla 4 los errores relativos en % para cada uno de los tres ejercicios propuestos.

Tabla 8. Resultados del Ejemplo 2.
Método	Para los datos de la Ref. [4]			Para los datos de la Ref. [5]			Para los datos de la Ref. [6]
	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B
MEF (Int. analítica)	0.693	200	111	70.14	2945.45	1636.36	0.067	1.90×105	1.04×105
Propuesto	1.016	270	150	101.5	4050	2250	0.179	4.86×105	2.70×105
Valor exacto (Timoshenko)	1.0275	300	150	102.625	4500	2250	0.18(Bernoulli)	5.40×105	2.70×105

Tabla 9. Errores relativos en % del Ejemplo 2.
Método	Ejercicio de la Ref. [3]			Ejercicio de la Ref. [4]			Ejercicio de la Ref. [5]
	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B
MEF (Int. analítica)	33	33	26	32	34	27	63	65	61
Propuesto	1.4	10	0	1.4	10	0	0.56	10	0

En cuanto a los valores de tensiones normales en la dirección vertical y las tensiones tangenciales, cabe hacer las mismas puntualizaciones que se hicieron en el ejercicio anterior.

8.1 Análisis de la convergencia

Para mostrar cómo a medida que la discretización se hace más tupida los resultados convergen hacia la solución teórica exacta, se ha calculado la misma estructura con modelados de 10, 15 y 20 elementos y con los datos de cada una de las referencias. Los resultados se muestran en la Tabla 10, en la que se repiten los obtenidos anteriormente para una discretización de 5 elementos y en la que se puede observar cómo todos los valores se acercan a su valor exacto al aumentar el número de elementos, siendo de resaltar cómo, según ya se comentó en el ejemplo anterior, la tensión en el punto B es exacta aún para sólo 5 elementos, pero no así en el punto A, donde la tensión se va haciendo más precisa a medida que aumenta la densidad de la discretización.

Tabla 10. Análisis de convergencia del Ejemplo 2.
Nº de elementos	Para los datos de la Ref.3			Para los datos de la Ref.4			Para los datos de la Ref.5
	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B	Flecha en C	Tensión σxen A	Tensión σxen B
5	1.016	270	150	101.5	4050	2250	0.179	4.86×105	2.70×105
10	1.0235	285	150	102.25	4275	2250	0.18063	5.13×105	2.70×105
15	1.02483	290	150	102.389	4350	2250	0.18088	5.22×105	2.70×105
20	1.02537	292	150	102.438	4387.5	2250	0.18097	5.265×105	2.70×105
Valor exacto (Timoshenko)	1.0275	300	150	102.625	4500	2250	0.1811	5.4×105	2.70×105

9. Ejemplo práctico 3

La estructura de la Figura 4, extraída de la Ref. [6],^§4.5 es una viga simplemente apoyada sometida al peso propio, con un peso total de 4000 unidades de carga y los mismos valores de las restantes características de la estructura anterior. Se pretende calcular la flecha y tensión en el punto central A de la cara inferior de la viga para una modelización de 5 elementos.

Figura 4. Viga simplemente apoyada sometida a peso propio. Dimensiones, numeración y de elementos.

Siguiendo los mismos pasos expuestos en el §7 y extrayendo del vector solución s los valores que se desean calcular se confecciona la Tabla 11 donde se comparan los valores obtenidos por el método tradicional y el método propuesto

Tabla 11. Resultados del Ejemplo 3 para 5 elementos.
Procedimiento	Flecha en A	Tensión σxen A
MEF (Integ. analítica)	-0.01089	1.013×105
Propuesto	-0.02746	2.64×105
Valor exacto	0.03125	3×105

9.1 Análisis de la convergencia

Finalmente se calcula la misma estructura para modelizaciones de 10, 15 y 20 elementos con objeto de analizar la convergencia de resultados, obteniendo los valores de la Tabla 12 en la que se repiten los obtenidos para 5 elementos y ya expuestos en la tabla anterior y en la que nuevamente se observa la tendencia de los valores hacia la solución exacta al aumentar el número de elementos.

Tabla 12. Análisis de convergencia del Ejemplo 3.
Nº de elementos	Flecha en A	Tensión en A
5	-0.02746	2.64×105
10	-0.031531	2.94×105
15	-0.031824	2.95×105
20	-0.032056	2.98×105
Valor exacto	-0.032153	3×105

10. Ejemplo práctico 4

Como último ejemplo, un tanto diferente a los anteriores, se analiza la estructura de la Figura 5, extraída de la Ref. [7],^§6.4.1 modelizándose como una viga de gran canto simplemente apoyada y sometida a carga triangular en su cara superior. Debido a la simetría geométrica y de acciones, se utilizará la mitad izquierda para la determinación de los desplazamientos de los nudos y distribución de tensiones, para lo cual bastará con seguir los pasos del §7 con los cambios pertinentes de los datos y de la modelización adoptada. Puesto que en dicha referencia se detallan profusamente todos los pasos -excepto el cálculo de las matrices elementales mediante la integración- que el método requiere, pueden contrastarse éstos con los de esta metodología con objeto de comparar el grado de complejidad de que cada uno de los procedimientos.

Figura 5. Viga de gran canto bajo carga triangular.

Adoptando el mismo criterio que el de la referencia, se analiza la estructura para modelizaciones de 4, 9, 16 y 36 elementos, de los que se detalla el primero de ellos, de acuerdo con la Figura 5 (b). Como ya se comentó en la exposición teórica, §3.6, se puede resolver el sistema (26), esto es, sin introducir condición estática perimetral alguna, con lo que se obtienen unos resultados ya mejorados respecto al método tradicional (p.e. la flecha en el nudo 3 sería de -0.9163×10^-4 frente a -0.8754×10^-4 que figura en el original); pero puede comprobarse cómo a medida que se van introduciendo dichas condiciones los resultados van tendiendo a la solución exacta, pues cada condición añadida representa una nueva información para el sistema a resolver, por lo que éste devolverá unos resultados cada vez más precisos. De acuerdo con esto, observando la Figura 5 (b), las tensiones normales en la dirección horizontal de todos los nodos del lateral izquierdo han de ser nulas (el método tradicional asigna una tensión en el nudo 1 de valor 385.185, es decir, muy lejos del valor real) y al introducir dichas condiciones en el sistema la flecha en el nudo 3 pasa a valer -1.0582×10^-4, esto es, más correcta que la anterior. Si además se añaden las condiciones de tensiones tangenciales nulas en los nodos del lado derecho debido a la simetría de la estructura (el método tradicional asigna una tensión tangencial en el nudo 9 de valor -329.167, cuando debería ser nulo), dicha flecha pasa a valer -1.3515×10^-4, valor más preciso que los dos anteriores. Así, pues, las matrices elementales de restricciones estáticas perimetrales a introducir en el §7.8 han de contemplar que para los elementos 1 y 3: ${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{i}}}$ = ${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{m}}}$ = 0; y para el 2 y el 4: ${\textstyle {\mbox{τ}}^{\mbox{j}}}$ = ${\textstyle {\mbox{τ}}^{\mbox{k}}}$ = 0. Con objeto de comparar resultados, en la Tabla 13 se muestran los valores de los desplazamientos de los nudos con el método tradicional y con el de la propuesta, para una modelización de 4 rectángulos, en la que se ha resaltado en negritas el valor de la flecha en el nudo 3, punto medio del lado inferior de la estructura original y cuyo valor exacto es 1.86×10^-4.

Tabla 13. Desplazamientos de los nudos para modelización de 4 elementos del ejemplo 4.
Procedimiento	u1	v1	u2	v2	u3	v3	u4	v4	u5	v5	u6	v6	u7	v7	u8	v8	u9	v9
MEF (Integ. numér.)	-0.3976	0	-0.2692	-0.7152	0	-0.8754	0.0457	-0.3287	0.0191	-0.7037	0	-0.9028	0.3062	-0.4259	0.2309	-0.7757	0	-0.9951	×10-4
Propuesto	-0.9456	0	-0.3149	-0.980	0	-1.3515	-0.0330	-0.4102	0.0237	-0.9152	0	-1.4126	0.5429	-0.4926	0.2767	-0.9945	0	-1.5059

Como se observa en la tabla, así como en todas las anteriores, los desplazamientos del método propuesto son siempre mayores y más correctos que los del tradicional, de lo que se infiere que en esta metodología los elementos adoptan un comportamiento menos rígido que en aquélla. Así, p.e., el error en la flecha del nudo 3 (en negritas en la Tabla) en el procedimiento tradicional es del 52.93%, mientras que para el método propuesto es del 27.34, es decir se reduce prácticamente a la mitad de aquél. Las tensiones normales en la dirección horizontal en los nudos del eje de simetría valen ${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{3}}=982.469}$ , ${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{6}}=-50.26}$ y ${\textstyle {\mbox{σ}}_{\mbox{x}}^{\mbox{9}}=-826,}$ por lo que la distribución de tensiones en dicha sección sería el que se muestra en la Figura 6, donde también se incluye la gráfica correspondiente al método tradicional. Comparando ambas gráficas se puede observar que la correspondiente al método propuesto se acerca algo más a la distribución real de las tensiones, de acuerdo con los valores de los respectivos momentos que figuran en la Tabla 14.

Figura 6. Distribución de tensiones normales en la sección central.

Si se supone, como en la referencia, despreciable la tensión del nudo central 6 frente a las de los otros dos, el momento flector en dicha sección sería de 120.56, siendo de 99.99 el que se obtiene por el MEF clásico y de 133.33 el exacto, con errores, pues, del 9.6% y 25%, respectivamente.

10.1 Análisis de la convergencia

Por último se confecciona la Tabla 14 en la que figuran las flechas del nudo 3 y el momento en la sección central de la estructura para modelizaciones de 3, 9, 16 y 36 elementos, así como los errores para cada discretización. Es de señalar que el método propuesto proporciona valores mucho más precisos que el tradicional para una misma discretización, sobre todo para las tensiones, parámetros de primordial importancia en el análisis por elementos finitos; véase, p.e., cómo los valores con 16 elementos con el método propuesto superan en exactitud a los obtenidos con 36 por el MEF.

Tabla 14. Flecha en el nudo 3 y momento en la sección central
Elementos	Procedimiento	v3(×10-4)	Error %	Momento	Error %
4	MEF	-0.8754	52.94	99.99	25.01
	Propuesto	-1.3515	27.34	120.56	9.58
9	MEF	-1.0794	41.97	116.68	12.49
	Propuesto	-1.3670	26.50	126.16	5.38
16	MEF	-1.2082	35.04	123.48	7.39
	Propuesto	-1.4155	23.90	129.14	3.14
36	MEF	-1.3758	26.03	128.95	3.29
	Propuesto	-1.5256	17.98	131.28	1.54

11. Dimensiones de los sistemas

A primera vista se puede tener la impresión de que la metodología propuesta adolece del inconveniente de exigir un sistema de ecuaciones de dimensiones muy superiores a los del método tradicional, y ello efectivamente es así, pues si se considera un dominio cualquiera discretizado con n nudos y e elementos y utilizando, p.e., el sistema reducido (33) cuya resolución devuelve de forma directa las tensiones, este sistema tendría unas dimensiones de 2n+12e más el número de Multiplicadores de Lagrange consecuencia de las restricciones, mientras que el MEF tradicional requeriría sólo de 2n ecuaciones, menos el número de restricciones de apoyo, es decir, muchas menos que aquél. Pero esta comparación solamente sería válida si se obtuviesen soluciones igualmente aceptables con ambos procedimientos, lo cual no es así, como ampliamente se ha demostrado en todos los ejemplos numéricos expuestos. Por lo tanto es obligado plantearse tal comparación en los términos de qué dimensiones han de tener los sistemas de cada procedimiento para que su resolución devuelva unos resultados equivalentes y, además, aceptables; y desde esta óptica, la afirmación de que el método propuesto adolece de tal inconveniente no será, en general, cierta. De los ejemplos antes expuestos fijémonos p.e., en los ejercicios 2 y 3 y extraídos de la Ref. [6], ya que en ésta figuran las gráficas que relacionan los resultados en función del número de g.d.l:

- Para el ejemplo 2 (§8): adoptando el sistema (33), la flecha en el punto C con la metodología propuesta y una modelización de 5 elementos vale 0.1793 (valor que puede considerarse aceptable, ya que el exacto es 0.18); la matriz de dicho sistema tiene unas dimensiones de 123×123, y en cuanto a la tensión en el punto B este sistema devuelve el valor exacto, que es 2.70×10⁵. Para conseguir un valor semejante e igualmente aceptable para la flecha por el MEF (0.1788) se requiere una modelización de 1192 elementos, que origina un sistema de dimensiones 2044×2044, esto es, 16.6 veces mayor que aquél; el valor de la tensión en B con este sistema es 2.5798×10⁵, menos preciso que el que se ha obtenido anteriormente. El sistema de 123 ecuaciones devuelve, pues, mejores resultados que el de 2044, sobre todo para las tensiones, magnitudes prioritarias en todo análisis estrutural.

- Para el ejemplo 3 (§9): el valor de la flecha para una modelización de 10 elementos por el método propuesto (-0.031531, según la Tabla 12) es equivalente al del MEF con 1192 elementos (-0.0315, según la referencia), lo que suponen sendos sistemas de ecuaciones de dimensiones 233×233 para el primero y de 2599×2599 para el segundo, 11 veces mayor. En cuanto a las tensiones en el punto A los resultados son de 2.94×10⁵ y 2.877×10⁵ para los respectivos procedimientos, siendo el valor exacto 3×10⁵.

Es necesario recordar, además, que la resolución del sistema con la metodología propuesta devuelve de forma directa los valores de las tensiones en todos los nodos y las reacciones en los apoyos (amén de los Multiplicadores de Lagrange, magnitudes que, en general, no tendrán utilidad desde el punto de vista práctico) mientras que el sistema tradicional solamente proporciona los desplazamientos de los nudos y, por tanto, será obligado un post-proceso para determinar las tensiones, magnitudes imprescindibles en todo análisis estructural.

12. Consideraciones últimas

Como ya se indicó, esta investigación tiene como objetivo fundamental la exposición de una metodología para el análisis y cálculo por elementos finitos totalmente distinta -y por ello novedosa- respecto al método tradicional universalmente conocido. Para comprobar si el método es válido para todo tipo de elemento finito elegido en la discretización del dominio en estudio sería imprescindible un detallado estudio para cada uno de dichos tipos de elementos, pero como quiera que éstos pueden adoptar infinidad de geometrías y un número arbitrario de nodos en cada una de ellas, sería imposible un análisis general que las incluya a todas, y, por ello, se ha de elegir una cualquiera que sirva de modelo para esta exposición. De dichos elementos, el Triángulo de Turner (CST) es el más simple, pero sus especiales características respecto a los demás -concretamente la hipótesis de tensión constante en todos sus puntos- lo hacen totalmente inadecuado para una descripción completa del funcionamiento de la metodología que aquí se propone. En efecto, como ya se ha visto a lo largo de este artículo, una de las novedades más significativas de este procedimiento (amén de la de evitar toda integración) es la inclusión de dos tipos de condiciones que para una modelización más cercana el hecho físico obligatoriamente han de cumplir las magnitudes estáticas asociadas a los nodos, condiciones a las que ha de renunciar el MEF habitual: 1º las tensiones han de ser tales que cada uno de los elementos que discretizan al dominio ha de estar en equilibrio, y 2º que en aquellos nodos del perímetro en los que no existan cargas aplicadas tales tensiones han de ser nulas. Pero ninguna de las dos condiciones son aplicables a dicho elemento, pues imponer la primera sería introducir una redundancia en el sistema (el equilibrio del elemento CST está intrínsecamente asegurado debido a dicha hipótesis), mientras que la introducción de la segunda sería una contradicción al tratar de imponer a los nodos de un elemento situados en el contorno unas tensiones diferentes a las tensiones de éste en los demás puntos de su geometría, lo que es contrario a dicha hipótesis de tensión constante. Por ello, para discretizar el dominio, se ha elegido como modelo para este estudio el rectángulo de 4 nodos, el más simple de todos los restantes. A lo largo de esta exposición se ha demostrado que el comportamiento de este elemento bajo esta nueva metodología mejora considerablemente los resultados respecto a la metodología tradicional con el mismo elemento, de lo que no puede en absoluto inferirse que para otros tipos de elementos se hayan de repetir las mismas conclusiones que para éste, aunque en principio no hay por qué descartarlas, debido a la ya reiterada circunstancia de que, sea cual fuere el elemento elegido, el sistema de ecuaciones de la propuesta contendrá unas condiciones que describirán mejor el hecho físico real, condiciones a las cuales no tiene acceso el habitual MEF. Si, no obstante, se comprobara que el comportamiento de otros elementos no fuese satisfactorio para este nueva metodología, no por ello se ha de renunciar al elemento aquí utilizado y que se ha mostrado altamente eficaz, ya que éste siempre se podría usar al menos para los dominios que lo permitiesen, que serían aquéllos en los que todos los ángulos de los vértices de su perímetro fuesen rectos (y si además tuviesen huecos internos, los lados de éstos fuesen paralelos a dicho perímetro). Hasta ahora sólo se ha estudiado también el elemento triangular de tres nodos (CST), el de seis nodos (LST) y el rectángulo de ocho nodos (Q8), con unas conclusiones similares a las del elemento rectangular aquí estudiado (y que, naturalmente, no pueden incluirse aquí por motivos de espacio) por lo que para una geometría irregular del contorno del dominio podría utilizarse el elemento rectangular de cuatro nodos para discretizar la mayor parte de él y en las regiones de su perímetro -que normalmente requerirán un número de elementos muy inferior a la del resto- se podría utilizar cualquiera de los dos elementos triangulares antes indicados, salvo que se compruebe que el elemento cuadrilátero irregular necesario para modelar también tales regiones se muestre igual de adecuado que el rectangular, lo que requerirá un estudio totalmente independiente del aquí expuesto y que, en su momento, habrá de ser objeto del correspondiente escrutinio.

13. Conclusiones

No debe extrañar que la metodología aquí expuesta proporcione unos resultados más cercanos al hecho físico que el del MEF tradicional, no solamente porque evita la integración numérica, origen de algunas inexactitudes, sino fundamentalmente porque permite describir la estructura real con un mayor grado de fiabilidad, ya que el hecho de partir de sistemas de ecuaciones en los que todas las incógnitas vienen dadas de forma explícita, posibilita la introducción de todo tipo de condiciones a cualesquiera de ellas, sean del tipo que fueren, lo que le hace más versátil y general para su aplicación a todo tipo de análisis por elementos finitos.

Agradecimientos

Agradezco la generosa colaboración del delineante J. Antonio Muñoz por los dibujos y la de mi colega J. Lázaro Bailón por la redacción del ABSTRACT en inglés.

Referencias

[1] Carlos A. Felippa. “Introduction to Finite Element Methods”. Department of Aerospace Engineering Sciences and Center for Aerospace Structures University of Colorado, Boulder, Colorado 80309-0429, USA, (2004).

[2] J. Tomás Celigüeta Lizarra. “Método de los Elementos Finitos para Análisis Estructural”. Campus Tecnológico de la Universidad de Navarra, San Sebastián, (2011).

[3] F. Martín Chica. “Propuesta de un método matricial bajo restricciones cinemáticas” Tesis Doctoral, Universidad de Granada. Abril de 2009.

[4] Robert D. Cook, David S. Malkus, Michael E. Plesha, Robert J. Witt. “Concepts and application of Finite Element Analisis”, 4ª Ed. University of Wisconsin-Madison. JOHN WILEY & SONS, INC. (2002).

[5] O.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu “The Finite Element Method: Its Basis & Fundamentals” 7ª Ed. ELSEVIER (2013).

[6] Eugenio Oñate. “Structural Analysis with the Finite Element Method. Linear Statics” Volume 1. Basis and Solids. CINE (2009).

[7] Walter Wunderlich y Walter D. Pilkey. “Mechanics of Structures. Variational and Computational Methods”. CRS PREESS, London (2003).

[8] Juan José Ródenas y otros. “Sobre la necesidad de controlar el error de discretización de elementos finitos en optimización de forma estructural con algoritmos evolutivos”. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. (2011).