## 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.

## Resumen

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

## 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 n1, n2, n3...) son los puntos donde inciden uno o más nodos y sólo pueden sufrir desplazamientos.

## 2. Objetivos

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 Figs. 1a) y 1b), 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 (Fig. 1c) podrán obtenerse en la forma [1]§7.3.2.

 ${\displaystyle {\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 Xi e Yi 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

 ${\displaystyle {\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 (1) se podrán escribir matricialmente en la forma

 ${\displaystyle \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 ${\displaystyle {\mbox{F}}_{\mbox{i}}^{\mbox{e}},}$${\displaystyle {\mbox{F}}_{\mbox{j}}^{\mbox{e}},}$${\displaystyle {\mbox{F}}_{\mbox{k}}^{\mbox{e}}}$ y ${\displaystyle {\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

 ${\displaystyle \mathbf {{\mbox{F}}_{\mbox{i}}^{\mbox{e}}} =\mathbf {{\mbox{I}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} }$
(4)

 ${\displaystyle \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

 ${\displaystyle \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)
 ${\displaystyle \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)
 ${\displaystyle \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 (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

 ${\displaystyle \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 ${\displaystyle \mathbf {{\mbox{ε}}_{\mbox{(x,y)}}^{\mbox{e}}} }$ y ${\displaystyle \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

 ${\displaystyle \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 Fig.1. Particularizando la expresión (10) para cada uno de los cuatro nodos del elemento, la (9) desarrollada tendrá la forma

 ${\displaystyle \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 de de ellos. Particionando la matriz Be anterior se obtienen las expresiones

 ${\displaystyle \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];}$${\displaystyle {\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];}$${\displaystyle {\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];}$${\displaystyle {\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 (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

 ${\displaystyle \mathbf {{\mbox{σ}}_{\mbox{n}}} =\mathbf {{\mbox{D}}_{\mbox{n}}} {\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}} }$
(13)

donde ${\displaystyle \mathbf {{\mbox{σ}}_{\mbox{n}}} }$ y ${\displaystyle \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

 ${\displaystyle \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}=}$${\displaystyle {\frac {1}{1-{\nu }^{2}}};{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}f_{2}=}$${\displaystyle \nu {\mbox{ }}f_{1};{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}{\mbox{ }}f_{3}=}$${\displaystyle {\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 ${\displaystyle \mathbf {{\mbox{ε}}{\mbox{'}}_{\mbox{n}}} =E{\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}} }$ la (13) vendrá dada en la forma

 ${\displaystyle \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}}} =}$${\displaystyle \mathbf {{\mbox{K}}_{\mbox{n}}} {\mbox{ }}\mathbf {{\mbox{ε}}_{\mbox{n}}^{'}} }$
(15)

siendo

 ${\displaystyle \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 Ke a la matriz diagonal por bloques de dimensiones 4×4

 ${\displaystyle \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 (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

 ${\displaystyle \mathbf {{\mbox{ε}}{\mbox{'}}^{\mbox{e}}} {\mbox{=}}\mathbf {{\mbox{B}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{d}}{\mbox{'}}^{\mbox{e}}} }$
(18)

donde ${\displaystyle \mathbf {{\mbox{ε}}{\mbox{'}}^{\mbox{e}}} =E{\mbox{ }}\mathbf {{\mbox{ε}}^{\mbox{e}}} }$ y ${\displaystyle \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 ${\displaystyle \mathbf {\mbox{d'}} ,}$ tensiones σ y deformaciones ${\displaystyle \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)

 ${\displaystyle \mathbf {{\mbox{R}}_{\mbox{n}}{\mbox{ }}{\mbox{d}}_{\mbox{n}}^{'}={\mbox{0}}} }$
(19)

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

 ${\displaystyle \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 un y vertical vn), la segunda a un apoyo móvil con plano de rodamiento horizontal (nulidad de desplazamiento vertical vn) y la tercera a un apoyo móvil con rodamiento vertical (desplazamiento horizontal un 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 ${\displaystyle \mathbf {\mbox{d'}} }$ que, junto con las tres anteriores y aplicando el método de los Multiplicadores de Lagrange, forman el sistema matricial

 ${\displaystyle \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 ${\displaystyle \mathbf {\mbox{d'}} }$ , σ, ${\displaystyle \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 vH 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:

 ${\displaystyle \mathbf {{\mbox{F}}_{\mbox{i}}^{\mbox{e}}} +\mathbf {{\mbox{F}}_{\mbox{j}}^{\mbox{e}}} +}$${\displaystyle \mathbf {{\mbox{F}}_{\mbox{k}}^{\mbox{e}}} +\mathbf {{\mbox{F}}_{\mbox{m}}^{\mbox{e}}} =}$${\displaystyle (\mathbf {{\mbox{I}}^{\mbox{e}}} +\mathbf {{\mbox{J}}^{\mbox{e}}} +\mathbf {{\mbox{K}}^{\mbox{e}}} +}$${\displaystyle \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 Fig. 1c)

 ${\displaystyle \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 (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

 ${\displaystyle \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 Te a la matriz de coeficientes, ésta vendrá dada, pues, por

 ${\displaystyle \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

 ${\displaystyle \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 (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 ( ${\displaystyle {\sigma }_{y}^{i}={\sigma }_{y}^{j}=0}$ ) y k-m ( ${\displaystyle {\sigma }_{y}^{k}={\sigma }_{y}^{m}=0}$ ) de un elemento e perteneciente al perímetro vendrá dada, respectivamente, por las ecuaciones homogéneas

 ${\displaystyle \mathbf {{\mbox{P}}_{\mbox{i-j}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =}$${\displaystyle \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}}} =}$${\displaystyle \mathbf {0} }$
(27)

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

 ${\displaystyle \mathbf {{\mbox{P}}_{\mbox{i-m}}^{\mbox{e}}} {\mbox{ }}\mathbf {{\mbox{σ}}^{\mbox{e}}} =}$${\displaystyle \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}}} =}$${\displaystyle \mathbf {0} }$
(28)

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

 ${\displaystyle {\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 Ie, Je, Ke o Me, 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.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 vH 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

 ${\displaystyle {\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 ${\displaystyle \mathbf {\mbox{ε'}} }$ del segundo bloque de (30) y sustituyéndolo en el tercero, se obtiene la ecuación matricial

 ${\displaystyle \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

 ${\displaystyle \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

 ${\displaystyle \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)

 ${\displaystyle \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

 ${\displaystyle \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 ${\displaystyle \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 ${\displaystyle \mathbf {{\mbox{λ}}_{\mbox{R}}} ;}$ e incluso de aquellos elementos de ${\displaystyle \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 Fig.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 Fig. 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 I0, J0, 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 vH, vB, vK, vR, vT y vP a los vectores asociados a las respectivas matrices H, B, K, R, T y P e Id 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 Fig.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 Fig.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 Fig.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 Fig.5b). 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 Fg.5b), 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 Fg.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

- 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×105. 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×105, 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×105 y 2.877×105 para los respectivos procedimientos, siendo el valor exacto 3×105.

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.

## 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.

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).

### Document information

Published on 29/03/18
Accepted on 01/03/18
Submitted on 21/09/17

Volume 34, Issue 1, 2018
DOI: 10.23967/j.rimni.2018.03.002