Otro MEF: El elemento Triangular de 3 Nodos.

Another FEM: the Triangular element with 3 nodes.

Autor: F. Martín Chica (Octubre de 2017).
[email protected]

SUMMARY
The procedure for every finite element that is shown in this article use the Mix Method which is consisten-
ted by this way: 1) by the assembly of elemental matrix that it only depend on the mechanic and geometri-
cal features of the element. You achieve a system of lineal ecuations coupled with the unknown for displa-
cement, tension and deformation of the nodes of the structure, and 2) In this system, you can add three
kind of conditions by the Method of Lagrange Multipliers, cinematics restrictions of fastening, statics con-
ditions of balance of the whole elements and statics conditions of the border. The resolution of this sys-
tem show you the valor of the whole unknown: 1) Displacement of the nodes, 2) tension on that nodes, 3)
deformation, 4) reactions, 5) Lagrange multipliers that enforce compliance the balance in all the isolated
element, and 6) ditto that enforce compliance the statics conditions of the border. Formulation used does
not require any integration type.
Keywords: matrix, finite element, balance, compatibility, constitutive, restriction.

RESUMEN
Con carácter general, es decir, para cualquier tipo de elemento finito elegido, el procedimiento que aquí se
presenta utiliza el Método Mixto que, en esencia, consiste en lo siguiente: 1) mediante el ensamblaje de
unas determinadas Matrices Elementales -que sólo dependen de las propiedades mecánicas y geométricas
del elemento - se consigue un sistema lineal de ecuaciones acopladas en las incógnitas de desplazamiento,
tensión y deformación de todos los nodos de la estructura, y 2) a dicho sistema se le añaden tres tipos de
condiciones mediante el Método de los Multiplicadores de Lagrange: restricciones cinemáticas de susten-
tación, condiciones estáticas de equilibrio de todos los elementos y condiciones estáticas de borde. La
resolución de dicho sistema devuelve directamente el valor de todas las incógnitas del problema: 1) despla-
zamientos de todos los nodos, 2) tensiones en todos ellos, 3) deformaciones de éstos, 4) reacciones, 5) Mul-
tiplicadores de Lagrange para que se cumpla el equilibrio de cada elemento aislado y 6) id. id. para que se
cumplan las condiciones estáticas de borde. La formulación utilizada no requiere, pues, ningún tipo de inte-
gración.
En este trabajo aplicaremos todo lo anterior al elemento Triangular de 3 Nodos (también conocido como
Triángulo de Turner y como elemento CST), incluyéndose al final una serie de ejercicios que muestran
cómo la metodología que aquí se propone proporciona los mismos resultados que el tradicional MEF -si se
obvian las pequeñas discrepancias numéricas de precisión- puesto que ambos procedimientos parten de
idénticas hipótesis.
Palabras clave: matriz, elemento finito, equilibrio, compatibilidad, constitutivas, restricciones.

Introito.- ¡Ah, el Método de los Elementos Finitos!. En verdad que desde que se principia a penetrarlo no
es extraño que de forma involuntaria y, por tanto, inevitable, afluya de inmediato a la memoria el famoso so-
neto de D. MIguel, y sobre todo su primer cuarteto, exclamando con él:
¡Voto a Dios, que me espanta esta grandeza
y que diera un doblón por describilla,
porque ¿a quién no sorprende y maravilla
esta máquina insigne, esta riqueza?!
¿o acaso no es ingenioso el invento de tratar de calcular el valor que unas desconocidas magnitudes tomarán
en un punto genérico (x,y) del plano a partir de otros valores… ¡también desconocidos!?. ¿No recuerda esto al
1
pez que se muerde la cola, a la serpiente que trata de engullirse a sí misma, al movimiento continuo o, peor
aún, a un círculo vicioso?.
Ciñéndonos por ahora exclusivamente al estudio del comportamiento mecánico de un material elástico en el
espacio 2D sometido a una determinada excitación contenida en su plano, el primer cimiento en el que se sus-
tenta el MEF es una función ¡desconocida! f(x,y) -función de interpolación- con la que se pretende calcular el
desplazamiento, la tensión y la deformación en cualquier punto en función -valga aquí la repetición de la pala-
bra- de los valores -ojo, ¡de nuevo desconocidos!- que tales magnitudes toman en otros puntos del plano.
Todo ello parece, sí, una aporía; pero hete aquí que, ¡oh maravilla!, en socorro de ella aparece milagrosamente
el príncipe encantado, el séptimo de caballería, el ángel salvador; esto es, la relación fundamental que propor-
ciona la teoría de la elasticidad y que despeja totalmente la senda al deshacer el nudo gordiano, al romper tal
círculo vicioso: la deformación en un punto no es más que la derivada del desplazamiento de dicho punto. ¡Me-
nos mal!. Primer problema resuelto.
El segundo cimiento del método no es más que el inapelable PTV: puesto que la energía ni crear se puede ni
tampoco destruir, necesariamente ha de colegirse que el trabajo que realizan las excitaciones debe ser idéntico
al que efectúan las tensiones! ¡Maravilloso principio! Segundo problema resuelto.
Con esos dos firmes pilares, a los que se les ha de añadir la imprescindible Ley de Hooke, el resto del proce-
so consistirá simplemente en obligar a que se cumplan unas determinadas relaciones entre las excitaciones y
las tensiones -ecuaciones de Equilibrio-, otras relaciones entre las deformaciones y los desplazamientos -
relaciones de Compatibilidad-, unas relaciones entre las tensiones y las deformaciones -relaciones Constitu-
tivas- y, por último, unas determinadas condiciones de contorno, con lo que todo parece quedar atado y bien
atado. Pero, ¡oh, sorpresa!, resulta que, debido a las especiales características del Elemento Triangular de 3
Nodos que nos ocupa, no será en absoluto necesario la utilización del primer cimiento comentado anterior-
mente, esto es, la función de interpolación; de ahí el enigmático estrambote que cierra el poema:
Y luego, incontinente,
caló el chapeo, requirió la espada,
miró al soslayo, fuese y no hubo nada.

1. Introducción
Una vez discretizado un dominio en el espacio 2D mediante elementos con cualquier tipo de geometría, dos
elementos contiguos compartirán el segmento lineal común a ambos, esto es, una simple infinidad de pun-
tos. Para reducir este problema continuo a otro discreto es necesario crear en dicho lado común una serie de
puntos que, a efectos de cálculo por el MEF, serán los únicos contactos que existirán entre ambos elemen-
tos. A dichos puntos le llamaremos nodos del elemento, signados por las letras i, j, k…, cada uno de los
cuales estará sometido, en general, a tres componentes de tensión (σx, σy y τ) y a otras tantas de deformación
(εx, εy y γ) como consecuencia de unas determinadas excitaciones externas que actúan sobre el dominio,
mientras que llamaremos nudos del dominio (signados por un número natural n: 1, 2, 3…) a los puntos
donde inciden uno o más nodos y sólo pueden sufrir desplazamientos. Así pues, las magnitudes asociadas a
un nodo cualquiera (p. e. el nodo j) de un elemento e serán las tres componentes de tensión que vendrán
signadas por σ xj , σ yj , y τ j , para tal elemento y las tres de deformación εxj , ε yj , y γ j , mientras que las magni-
tudes asociadas a un nudo n serán exclusivamente las dos componentes de desplazamiento un y vn. No ten-
drá, por tanto, sentido hablar de tensión o deformación en un nudo, puesto que habrá en él tantas tensiones
y deformaciones diferentes como nodos conecten con tal nudo; de ahí que se suelan tomar para estas magni-
tudes en tal punto la media de las de los nodos que en él inciden.

2. Objetivos
De acuerdo con la Ref.1, para el análisis de toda estructura formada por elementos discretos -sean éstos de
cualquier tipo de geometría y con cualquier número de nodos- será en general imprescindible el uso de tres
tipos de ecuaciones básicas o de gobierno: a) Ecuaciones de Equilibrio de fuerzas que inciden en cada nudo
del dominio, 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 nodos, y c) Ecuaciones Constitu-
tivas del material, que relacionan las deformaciones con las tensiones en los nodos. A estas ecuaciones se
2
han de añadir unas determinadas restricciones de contorno, siendo las más habituales las Ecuaciones que
eviten el movimiento de la estructura como sólido rígido, esto es, las Condiciones Cinemáticas de sustenta-
ción aplicadas a los nudos del dominio. El tratamiento de este último tipo de ecuaciones por el Método de
los Multiplicadores de Lagrange1 permite obtener, junto con los tres tipos anteriores, un sistema cuadrado,
compatible y determinado 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 exclusivas de desplazamientos (Sistema Reducido del Método Directo de la
Rigidez); sin embargo, en Método Mixto[1] se mantienen de forma explícita todas las ecuaciones menciona-
das -y, naturalmente, sus correspondientes incógnitas- lo que permitirá introducir cualquier tipo de restric-
ción a cualquiera de dichas incógnitas, 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.
No obstante lo anterior, las especiales características del elemento triangular que nos ocupa en este trabajo -
y concretamente la hipótesis de tensión constante2- permite simplificar el proceso antes descrito, por lo que,
por una parte, no será necesario la introducción de las ecuaciones que obliguen al equilibrio de cada ele-
mento aislado, toda vez que éste se ha de cumplir de forma automática, y por otra que las condiciones está-
ticas de contorno también están implícitamente descritas por tal hipótesis.
Es, pues, objetivo de esta investigación el cálculo de estructuras en el espacio 2D, discretizadas mediante
elementos triangulares de 3 nodos utilizando el sistema de ecuaciones formado exclusivamente por las cua-
tro primeras antes descritas.
3. Herramientas del método
El procedimiento que aquí se propone es un método matricial, por lo tanto las unidades operativas que ma-
neja son matrices Elementales agrupadas en cuatro tipos: a) matrices de Fuerzas Nodales o de Equilibrio, b)
matrices de Compatibilidad, c) matrices Constitutivas, d) matrices de Restricciones Cinemáticas (restric-
ciones en aquellos aparatos de apoyo que sustentan a la estructura).
3.1. Matrices Elementales de Fuerzas Nodales o de Equilibrio.- Son aquellas matrices que permiten
determinar las dos componentes de fuerza concentrada (X, Y), en cada uno de los tres nodos, en función de
las tensiones existentes en el elemento. En la Fg.1 se muestran las coordenadas de los nodos del Triángulo
de 3 Nodos, de espesor t, así como las componentes de sus lados -respeto a unos ejes XY cuyo origen puede
ser cualquier punto arbitrario-, entendidas éstas como segmentos orientados, y donde se ha tomado para las
tensiones el criterio convencional de signos y se ha supuesto que el elemento no contiene cargas volumétri-
cas ni de superficie.
Signando por i, j, y k a los nodos del elemento en sentido levógiro, las componentes de Fuerza Nodal Equi-
valente en el nodo i, en función de las tensiones (Fig. 1a) y de las proyecciones de los lados (Fig. 1b), po-
drán obtenerse en la forma habitual
  by c y   b c  t
X i  t  σ x  σ x   t  τ x  τ x     by  c y  σ x  (bx  cx ) τ 
 2 2  2 2 2
(1)
 bx c x    by cy  t
Yi  t  σ y  σ y   t  τ  τ      bx  c x  σ y  (by  c y ) τ 
 2 2  2 2 2

donde Xi e Yi son, respectivamente, la componente horizontal y vertical de la fuerza concentrada en el nodo

i, componentes que sustituyen a las tensiones distribuidas a lo largo de los lados i-j e i-k del elemento. Con-
cibiendo, como se ha dicho, los tres lados a, b y c (respectivamente opuestos a los vértices i, j, k del triángu-
lo) como vectores con las direcciones indicadas en la Fg.1b, sus proyecciones sobre los ejes coordenados

1
¡Ay, cuán poco agradecido estará este señor con todos nosotros!
2
Hipótesis, por otra parte, un tanto laxa que supondrá una inevitable fuente de error, puesto que parece poco exigente
admitir que en un medio continuo los efectos de una determinada excitación sean idénticos en cualquiera de sus pun-
tos, por muy cercanos que se les considere: “The Turner triangle is not a good performer for structural stress
analysis”[2]§15.1.
3
 σy Y bx ax

Yk
Xk
τ yk
k k
a
by b ay
σx τ τ σx
Xj
j j
yj cy
i τ Xi i c
Yj
τ yi
Yi cx
σy
O xi xk xj
X
a) b)
Fig. 1.- Triángulo de Turner. a) Tensiones, b) Coordenadas de los nodos y proyecciones de los lados.

vendrán dadas por los segmentos orientados

ax = xk - x j , ay = yk - y j
bx = xi - xk , by = yi - yk (2)
cx = x j - xi , cy = y j - yi

pero sumando m.a.m. las anteriores se tiene que

ax +bx  cx = 0 y ay +by  cy = 0 (3)

y, por tanto, se cumple que

by  cy = -ay y bx  cx = -ax (4)

Llamando σe al vector columna que contiene las componentes de tensión existentes en todo punto del ele-
mento, dicho vector vendrá expresado como
T
σ e  σ x σ y τ  (5)

por lo que, de acuerdo con la (4), las expresiones (1) se podrán escribir en la forma matricial más simple y
mnemotécnica como
 X i   t  a y 0 a x  σ e (6)
 Yi  2  0 a x a y 

Signando por Fie , Fje , y Fke a cada uno de los vectores que contienen a las dos componentes de fuerza nodal
de los nudos i, j, y k, respectivamente, la anterior vendrá dada para el nodo i por
Fie  I e σ e (7)
donde
t  a y 0 a x 
Ie  (8)
2  0 a x a y 

es la matriz de Fuerza Nodal o matriz de Equilibrio asociada al nodo i, que, como se ha visto, permite cal-
cular las componentes de dicha fuerza en función de las proyecciones de los lados y de las tensiones a las
que está sometido el elemento.
Procediendo de forma similar para los nodos j, y k, las expresiones para las correspondientes matrices de
Equilibrio para los tres nodos vendrán dadas (por simple rotación de las designaciones a, b y c de los seg-
mentos orientados) por

4
 t  a y 0 a x  t  by 0 bx  t   c y 0 cx 
Ie  Je  Ke  (9)
2  0 a x  a y  2  0 bx by  2  0 cx  c y 

Las expresiones dadas en (9) son las unidades operativas que llamaremos matrices Elementales de Fuerzas
Nodales o de Equilibrio, 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 in-
cógnitas de tensiones de todos los elementos de la estructura, que obligará a que se cumpla el equilibrio
estático en cada uno sus nudos, siendo el vector de términos independientes de este sistema las acciones que
actúan sobre todos ellos, incluidos los que pertenezcan a los aparatos de apoyo[1]. Llamando Fe al vector
columna que describe todas las componentes de Fuerzas Nodales del elemento y He a la matriz -a la que
llamaremos matriz de Equilibrio del elemento3- que relaciona tales componentes con las tensiones del mis-
mo indicadas en (5), puede escribirse que
 Ie 
 
Fe = H e σ e siendo He   J e  (10)
 K e 

ecuación matricial, cuyo desarrollo será, de acuerdo con todo lo anterior, de la forma
 Xi   -a y 0 ax 
Y   0 a x -a y 
 i    σ x 
 X j  t  -by 0 bx   
    σy (11)
 Yj  2  0 bx -by   
τ
Xk   -c y 0 cx   
   
 Yk   0 cx -c y 

3.2. Matrices Elementales de Compatibilidad.- Los valores de las deformaciones y de los desplazamien-
tos no pueden ser independientes entre sí, por lo que necesariamente ha de existir un conjunto de ecuacio-
nes 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 expresión
ε e(x,y)  B (x,y)
e e
d (x,y) (12)

donde ε e(x,y) y de(x,y) 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 triangular que nos ocupa, viene dada por la conocida expresión
 y  yk 0 y k  yi 0 yi  y j 0 
1  j
Be(x,y)  0 xk  x j 0 xi  xk 0 x j  xi  (13)
2A  x  x y  y x  x y  y x  x y  y 
 k j j k i k k i j i i j

la cual, según las relaciones dadas en (2) y omitiendo por comodidad el subíndice en la designación de la
matriz, puede venir en la forma
a 0  by 0  c y 0 
1  y
Be  0 ax 0 bx 0 cx  (14)
2 A  a  a b b c  c 
 x y x y x y

matriz de compatibilidad4 que relaciona las deformaciones εe del elemento con los desplazamientos de de
los nudos en donde inciden sus nodos, mediante la ecuación matricial
ε e = Be d e (15)

3
Tensor de divergencia[2]§ 14.3.1
4
Operador de gradiente simétrico[2]§ 14.3.1
5
donde
ε e   ε x εv γ 
T
y d e   ui vi u j v j uk vk 
T
(16)

ecuación que, desarrollada, será

 ui 
v 
 ε x  1  a y 0  b y 0 c y 0  i 
u
ε y    0 ax 0 bx 0 cx   j  (17)
 γ  2A  a x a y b x  b y c x  v
   c y   j 
u
 k
 vk 
Particionando la matriz anterior se obtienen las submatrices
a b c
1  y
0 
1  y
0 
1  y
0 
ie  0 ax  je  0 bx  ke  0 cx  (18)
2 A  a x a y  2 A  bx by  2 A  cx  c y 
     

que llamaremos matrices Elementales de Compatibilidad asociadas a cada uno de los nodos del elemento.
La primera de ellas determina las deformaciones de éste debidas a los desplazamientos del nodo i, la segun-
da las deformaciones originadas por los desplazamientos del j y la tercera las originadas por los desplaza-
mientos del nodo k. De forma semejante a las matrices de Fuerzas Nodales descritas en el apartado anterior,
las dadas en (18) no son más que las unidades operativas con cuyo pertinente ensamblaje se obtendrá otra
matriz, matriz de Compatibilidad B de la Estructura, originándose un segundo conjunto de ecuaciones li-
neales simultáneas y homogéneas en las incógnitas de desplazamiento d de todos los nudos (incluidos los
que pertenezcan a aparatos de apoyo[1]) y de las deformaciones ε de todos los elementos que discretizan el
dominio, ecuaciones que obligarán al cumplimiento de la compatibilidad entre éstas y aquéllos. La expre-
sión de la matriz de compatibilidad B puede venir de modo más compacto, de acuerdo con (18) en forma
B =  i e je k e  (19)

Obsérvese la curiosa relación existente entre las unidades matriciales de Equilibrio dadas en (9) y sus res-
pectivas de Compatibilidad mostradas en (18): cada una de las primeras no son más que las transpuestas de
las respectivas segundas multiplicadas por el volumen del elemento, esto es
T T T
I e   ve i e  , J e   ve je  y K e   ve k e  (20)

o, lo que es igual
T T T
I e   V e i e  , J e   V e je  y K e   V e k e  (21)

siendo
 ve 0 0
Ve   0 ve 0 (22)
 
 0 0 ve 

una matriz diagonal de dimensión 3×3 donde los valores de tal diagonal son el volumen del elemento, esto
es ve = t A (A=área del elemento), por lo que en adelante se le llamará matriz elemental de Volumen5, donde
ve = t A es el volumen del triángulo. De acuerdo con esto, la matriz de Equilibrio dada en (10) puede expre-
sarse en la forma
  V e ie T   ie T  V e T   ie T 
          
H   V j     je   V e      je    V e 
 e e T 
e T T T T

     
(23)
 e e T  e T e T  e T
  V k     k   V     k  

pero teniendo en cuenta (19) y que Ve es diagonal, de la anterior se desprende que

5
Puesto que este es un nuevo personaje totalmente imprevisto, habrá que bautizarlo de alguna forma, ¿no?
6
 T T
H e   Be  V e por lo que  H e   V e B e (24)

y, puesto que Ve no puede ser singular, se tiene finalmente que

-1 T
Be   V e   He  (25)

La circunstancia de que se cumplan las relaciones (24) no ha de extrañar, pues ello no es más que la conse-
cuencia inmediata de la hipótesis adoptada para tensiones y deformaciones y del principio de conservación
de la energía; según este principio se ha de cumplir que el trabajo Te que realizan las Fuerzas Nodales apli-
cadas en los nodos del elemento durante los desplazamientos de de éstos ha de ser igual a la energía elástica
We almacenada por las tensiones σe durante las deformaciones εe del mismo, es decir, según este principio
se ha de cumplir que
1 e T e 1 T
Te =  F   d  = W e =  σ e  ε e  dve (26)
2 2 v e

pero de la (10) y la (15), y teniendo en cuenta que las tensiones y deformaciones son constantes en todo el
volumen del elemento, sustituyendo en la anterior y simplificando, podrá escribirse que
T T T
 σ e   H e  d e = σ e   ve Be  d e (27)

o, de acuerdo con (22)

T T T
 σ e   H e  d e = σ e   V e B e  d e (28)

de donde, por identificación de términos, se desprende que necesariamente

T
 H e   V e Be (29)

“relación de contragradiencia” para el elemento triangular de 3 nodos que, como puede comprobarse, es
idéntica a la dada en (24). Esta “relación de contragradiencia” es de suma importancia por cuanto ella per-
mite, como se muestra en (25), obtener la matriz de Compatibilidad B de la estructura en función de la de
Equilibrio H (o viceversa, según (24)), con lo que bastará con ensamblar las matrices dadas en (9) para
conseguir H, de lo que se desprende que, tal y como ya se indicó al principio de esta exposición, no es en
absoluto necesaria la determinación de las matrices de compatibilidad (18) ni, por tanto, la función de inter-
polación para el elemento de Turner, lo que simplificará el proceso de cálculo.
3.3. Matriz Elemental Constitutiva.- Como es bien sabido, la relación existente entre las tensiones y de-
formaciones para cualquier punto del elemento e viene dada por la ley de Hooke generalizada
σ e = De ε e (30)
siendo
 f1 f2 0 1 1
De  E  f 2 f1 0  ; con f1  ; f 2   f1 ; f 3  ; (31)
 0 0 
f3  1   2
2(1  )

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 (30) suele acarrear algún proble-
ma 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 ε'e  E ε e dicha ecuación vendrá dada
en la forma σ e = K e E ε e y, por tanto, la ecuación constitutiva adquiere la forma definitiva

σ e = K e ε'e (32)

Siendo
e
 f1 f2 0
K   f2 f1 0 (33)
 0 0 f 3 

7
unidad matricial 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 elementos de la estructura en estudio.
Para introducir en la (15) dicho cambio de variable, bastará con multiplicar esta ecuación por E, con lo que
la ecuación de Compatibilidad para el elemento quedará en la forma
ε'e = B e d'e (34)

donde ε'e  E ε e y d'e  E d 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 d', tensiones σ y deformaciones ε', cuya matriz de coeficientes será cuadra-
da, pero singular, puesto que aún no se han introducido las condiciones de sustentación que eviten el mo-
vimiento 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 tratamien-
to 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, 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 lineales o no, etc.) en este trabajo nos ceñiremos ex-
clusivamente 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 en el espacio 2D)
R n d'n  0 (35)

siendo d'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

Rn  
1 0
, ó R n   0 1 , ó R n  1 0 (36)
 0 1

donde el subíndice designará el nudo restringido, correspondiendo la primera de ellas a un apoyo fijo (nuli-
dad 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 (36) proporcionará una matriz R, Matriz de Restricciones Cinemáticas, que genera un primer gru-
po de restricciones en las incógnitas d' que, junto con las tres anteriores, formarán el sistema matricial
 0 H 0 R T   d'   v H 
 B 0 -I 0   σ   0 
 0 -I K 0   ε'    0  (37)
R 0 0 0   λ R   0 

que ahora ya es compatible y determinado, siendo respectivamente d' , σ, ε' y λR los vectores columna de
desplazamientos de todos los nudos (multiplicados por E), incluso apoyos, tensiones en todos los elemen-
tos, deformaciones en éstos (multiplicadas por E) y Multiplicadores de Lagrange asociados a las restriccio-
nes de los apoyos (reacciones con signo 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.
3.5. Matrices Elementales de Restricciones Estáticas internas.- El sistema (37) tiene, pues, solución;
pero para cualquier otro elemento finito que no sea el triángulo de Turner, dicha solución será inadmisible,
puesto que no está asegurado aún que se cumpla el equilibrio estático de cada uno de los elementos que
discretizan a la estructura. Sin embargo para el elemento que nos ocupa dicho equilibrio está intrínsecamen-
te asegurado, debido a la hipótesis de tensión constante adoptada, como se indicó anteriormente, y su inclu-
sión en el sistema no sería más que una redundancia puesto que se añadirían unas ecuaciones idénticamente
nulas, por lo que serán innecesarias dichas matrices elementales.

8
3.6. Matrices Elementales de Restricciones Estáticas externas o perimetrales.- Es evidente que, para un
correcto análisis de un continuo, es imprescindible además describir el estado tensional de su perímetro, lo
que, de forma semejante a lo expuesto en el apartado anterior, ha de ser obligado para todo tipo de elemento
finito. Pero, nuevamente debido a la dichosa hipótesis, el triángulo de 3 nodos no permite dicha condición,
ya que sería una contradicción imponer a los nodos de un elemento situados en el contorno del dominio
unas tensiones diferentes a las tensiones de éste. He aquí otra inevitable fuente de error. Las matrices ele-
mentales de Restricciones Estáticas para el elemento que nos ocupa no pueden, por tanto, introducirse en la
matriz del sistema, por lo que éste tendrá la configuración definitiva dada por la expresión (37).

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 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 del
sistema 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 elemen-
tos por i, j, y k, de acuerdo con la Fig.1 y partiendo de un arreglo de tantas filas como nudos n y tantas co-
lumnas como elementos e tenga la estructura, el ensamblaje de las matrices elementales dadas en (9) se
podrá abordar por filas o por columnas. En el primer caso se ubicará en la fila f asociada al nudo f las sub-
matrices Ie, Je o Ke, correspondientes a los elementos 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 la anteriores submatrices de dicho ele-
mento en la correspondiente fila (nudo) donde inciden cada uno de los cuatro nodos del mismo.6
4.2. Ensamblaje de la Matriz de Compatibilidad (B).- Evidentemente, para la consecución de la matriz
de compatibilidad, podría procederse al ensamblaje de las matrices elementales dadas en (18) en un arreglo
de tantas filas como elementos y tantas columnas como nudos, de forma recíproca a lo indicado en el apar-
tado anterior, pero ello no será necesario, como ya se ha indicado. En efecto, 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 B vendrá dada directamente por HT, sin necesidad, pues, de
proceder a un nuevo ensamblaje para la obtención de B7; pero en el triángulo de Turner tal relación se trans-
forma en la expresión dada por (29), por lo que la (25) para toda la estructura será
B  V -1H T (38)
expresión que evita el ensamblaje de la matriz de Compatibilidad B, pero que obliga al ensamblaje de la
matriz V, aunque éste es muy simple ya que V es un arreglo diagonal por bloques de dimensiones iguales al
número de elementos, siendo los elementos de tal diagonal cada una de las matrices Ve dadas en (22) para
cada elemento que discretiza al dominio.
4.3. Ensamblaje de la Matriz Constitutiva (K).- De la misma forma que la matriz V, el ensamblaje de
esta matriz es inmediato, pues no es más que un arreglo diagonal por bloques, de las mismas dimensiones
que aquélla, ubicando la expresión (33) para cada elemento 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 (36) (según sea el tipo de apoyo con el que se modele la susten-
tación) y en la columna correspondiente al nudo restringido.
4.5. 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

6
Los ejemplos que más adelante se presentan aclarará mejor dicho ensamblaje.
7
O, recíprocamente, ensamblar B y obtener H por transposición de B.
9
a las ecuaciones de equilibrio, que contendrá las componentes de acción que actúan en los nudos de la es-
tructura, esto es, si n es el número de nudos de la estructura, el vector vH vendrá dado por
v H =  a1 a 2 … an 
T
(39)

donde as (s = 1,2…n) es el vector de acciones que actúa en el nudo s.

5. Sistema Completo
Sustituyendo, pues, la (38) en la expresión (37) el Sistema Completo adopta la forma
d' σ ε' λR
   
(2 n ) (3e ) (3e ) (?)

Ecuaciones de equilibrio (2n )   0 H 0 R T   d'   v H  (40)

" de compatibilidad (3e )   V -1 H T 0 -I 0  σ    0 
" constitutivas (3e )   0 -I K 0   ε'   0 

Restricciones cinemáticas (?) 
 R 0 0 0   λ R   0 

en donde, para facilitar su lectura, se han colocado a su izquierda la denominación de cada uno de los blo-
ques de ecuaciones con sus respectivas dimensiones y en la parte superior las incógnitas y dimensiones
asociadas a cada columna. Las interrogaciones que figuran en el esquema anterior indican que el número de
ecuaciones e incógnitas de restricciones dependerá del número y tipo de aparatos de apoyo que sustentan a
la estructura.
Como puede observarse la matriz de coeficientes no es simétrica -como así ocurre en el Método Mixto[1]
para estructuras de barras- debido a la nueva “relación de Contragradiencia”, aunque, si se desea, puede
conseguirse fácilmente dicha simetría pre-multiplicando el segundo y tercer bloque de (40) por la matriz V,
con lo que se obtendrá el sistema simétrico, compatible y determinado

 0 H 0 R T   d'   v H 
 HT 0 -V 0  σ    0 
 0 Ms = v
-V V K 0   ε'   0 
ó (41)

 R 0 0 0   λ R   0 

siendo M la matriz de coeficientes, s el vector solución y v el vector de términos independientes.

6. Sistema Reducido
Aunque, por compleja que sea la estructura, la capacidad de los actuales (y sobre todo futuros8) ordenadores
puede ser presumiblemente suficiente para la resolución del Sistema Completo (41), éste puede reducirse si
se desea eliminado las deformaciones, ya que éstas no suelen ser magnitudes fundamentales en el análisis9.
En efecto, puesto que V no puede ser singular, despejando el vector de deformaciones ε' del segundo blo-
que de (40) y sustituyéndolo en el tercero se tiene que
ε' = V -1H T d', - Iσ + K V -1H T d' = 0 (42)
y, por tanto, el Sistema Reducido vendrá dado por

 0 H R T   d'   v H 
 K V -1H T -I 0  σ    0  (43)
 R 0 0   λ R   0 
 

que será cuadrado, compatible y determinado, aunque no simétrico. No obstante, puesto que K tampoco
puede ser singular, puede conseguirse la simetría. En efecto, pre-multiplicando el segundo bloque de (43)
por la inversa de K éste quedará en la forma

8
“Las tecnologías cuánticas” (National Geographic). Eduardo Arroyo, 2017, RBA Coleccionables, S.A.U. :“El desa-
rrollo de la computación cuántica anuncia una auténtica revolución en el campo de la tecnología.”
9
Lo que en absoluto ha de significar la renuncia a su cálculo, pues éste puede hacerse a posteriori.
10
 V -1H T d' - K -1σ = 0 (44)
y volviendo a pre-multiplicar la anterior por la matriz V se tiene
H T d' - V K -1 σ = 0 (45)
por lo que el Sistema Completo quedará modificado y en la forma siguiente

 0 H R T   d'   v H 
 HT - V K -1 0  σ    0  (46)
R 0 0   λ R   0 
 

Sistema Reducido que, por ser simétricas V y K, vuelve a ser simétrico, compatible y determinado, habién-
dose eliminado, pues, 3e 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 eleva-
do si se quiere conseguir una solución medianamente aceptable.
7. Elementos repetitivos
Aunque no es el objeto principal de este trabajo, se proponen aquí unas ideas que pueden ser útiles para la
simplificación del proceso de cálculo. En efecto, puesto que para conseguir tal precisión será necesario un
gran número de elementos para modelizar el dominio en estudio es muy normal que entre ellos se den rela-
ciones de semejanza y, por tanto, se puede evitar el cálculo de las matrices de Fuerzas Nodales para cada
uno de ellos, pues bastará con elegir uno cualquiera y a partir de éste determinar dichas componentes para
todos los demás, sin necesidad, pues, de describir las coordenadas de cada uno de los nodos de cada ele-
mento finito que discretizan a la estructura. Sea, pues, el elemento eP de la Fg.2 el que se ha elegido como
prototipo, cuyas matrices de Equilibrio están dadas en (9), y sea ex otro elemento simétrico al anterior res-
pecto al eje horizontal OX, cuya relación de semejanza se ha supuesto por brevedad que es la unidad.

ax bx bx ax
Y
yk
k k
a a
ay b by by b ay
ey ep
j j
cy i yj cy
c i c
cx
yi

O xi xk xj
cx cx X
cy c i i
c cy

j j
exy by by ex
b b
ay a aa ay
k k
ax bx bx ax

Fig. 2.- Relaciones de simetrías

        
Puesto que las proyecciones de los lados sobre el eje horizontal de ambos triángulo son iguales y del mismo
signo y las verticales de ex son de signo contrario respecto a las de ep bastará con cambiar de signo estas
últimas en las expresiones (9) para obtener las matrices de Fuerza Nodales del elemento ex, con lo que se
tendrá para este elemento
t a y 0 ax  t b 0 bx  t  c 0 cx 
I ex    ; J ex   y  y K ex   y (47)
2 0 a x a y 2 0 bx b y 2  0 c x c y 

11
De forma contraria, el elemento ey, simétrico al ep respecto al eje vertical, mantiene paridad de signos res-
pecto al prototipo para las proyecciones verticales, mientras que las horizontales son de signo contrario, por
lo que para obtener las matrices de Fuerzas Nodales para este elemento habrán de cambiarse los signos de
estas proyecciones en (9). Finalmente, el elemento exy es doblemente simétrico al ep respecto a ambos ejes,
por lo que las proyecciones son todas de signo contrario. De acuerdo con todo lo anterior, las matrices de
Fuerzas Nodales para ey vendrán dadas por
t  a y 0 a x  t  b 0 bx  t  c 0 cx 
Iey  ; J ey   y y K ey   y (48)
2  0 a x  a y  2  0 bx by  2  0 c x c y 

y para el elemento exy

t a y 0 a x  t b 0 bx  t c 0 cx 
I exy    ; J exy   y  y K exy   y (49)
2 0  a x a y  2 0  bx b y  2 0  c x c y 

Para facilitar el proceso de cálculo pueden crearse las siguientes matrices:

 1 0 1  1 0 1  1 0 1
Sx    Sy    S xy    (50)
 0 1 1  0 1 1   0 1 1

la primera de ellas cambia los signos de las proyecciones verticales, por lo que le llamaremos matriz de
simetría respecto al eje OX -de ahí el subíndice-, la segunda cambia los de las proyecciones horizontales,
matriz de simetría respecto al eje OY y la tercera cambia los signos de ambas proyecciones, matriz de doble
simetría. De esta forma las expresiones (47), (48) y (49) pueden escribirse, respectivamente, como
I ex = Ie  S x ; J ex = J e  S x y K ex = K e  Sx (51)

Iey = Ie  S y ; J ey = J e  S y y K ey = K e  S y (52)

I exy = I e  S xy ; J exy = J e  S xy y K exy = K e  S xy (53)

expresiones que permiten obtener las matrices de Equilibrio de los elementos ex, ey y exy respectivamente, en
función de las del prototipo dadas en (9), y en donde el asterisco es un operador que hace multiplicar las
matrices elemento a elemento.10 Si un elemento e´, además de una relación de simetría, tiene un factor de
semejanza r respecto al prototipo ep, esto es si a x = r a x bastará con multiplicar las expresiones (50) por el
factor r para obtener sus matrices de Equilibrio en función de las del prototipo.

8. Ilustrativo ejemplo 1
El ejercicio que se muestra en la Fig.3 está extraído de la Ref. 3§6.5, con los valores t = 1 in., E = 30×106 psi.
y ν = 0.3. Se trata de calcular los desplazamientos nodales y las tensiones de la placa con una discretización
de dos elementos triangulares de tensión constante sometida a una tracción continua en su lado derecho de
valor 1000 psi, por lo que las acciones nodales serán las que se muestran en los nudos 3 y 4 del dibujo.
20

2 3 5000
i k
1 j
10
2
j
k i 5000
1
4

Fig.3.- Voladizo a tracción simple. Discretización, cotas y acciones en los nudos.

10
Así lo entiende el programa Mathematica.
12
A continuación se expone el desarrollo completo del cálculo con el programa Mathematica 10.0, evitando
por comodidad de escritura subíndices y superíndices en los Input a introducir, habiéndose sombreado di-
chos inputs para facilitar su lectura.
8.1. Introducción de datos generales.- De acuerdo con el dibujo (n=4, e=2) y la expresión (31)11:

8.2. Coordenadas de los nodos del elemento prototipo y proyecciones de los lados.-Tomando como referen-
cia el elemento 1, las coordenadas de sus nodos respecto a un sistema de eje horizontal OX y vertical OY
en sentido levógiro con origen en el nodo i (nudo 2), así como las proyecciones de los lados del triángulo
según (2) vendrán dadas por

8.3.- Matrices elementales de Volumen, de Simetrías y de Equilibrio.- Llamando A1 al área del elemento 1,
v1 a su volumen y V1 a la matriz elemental de volumen dada en (23), teniendo en cuenta las (50) y que la
relación de semejanza de los dos elementos es la unidad, se tendrá12

con lo que las matrices de Equilibrio para el elemento prototipo vendrán dadas, según (9), por

y puesto que el elemento 2 es doblemente simétrico al 1, aplicando las expresiones (53)

se obtienen las matrices de Equilibrio para dicho elemento.

8.4. Matriz de Equilibrio H de la Estructura.- Partiendo de un arreglo de tantas filas como nudos (n=4) y
tantas columnas como elementos (e=2), a la vista de la Fig.3 y según el Apartado 4.1, la matriz H será

que para que Mathematica la lea y calcule su transpuesta, se ha de escribir

8.5. Matriz de Volumen V de la Estructura.- Según lo indicado en el Apartado 4.2., la matriz de Volumen
para toda la Estructura viene dada por un arreglo diagonal de dimensiones igual al de elementos, siendo los
términos de dicha diagonal las matrices elementales de Volumen de cada elemento que configura al domi-
nio, ya calculadas en el Apartado 8.3. Así pues

11
Al módulo de elasticidad E se le ha designado como Ee por estar protegida la letra E para Mathematica.
12
Aunque aquí sólo es necesaria la matriz Sxy por ser el elemento 2 doblemente simétrico respecto a OX y OY, se
introducen las tres por generalizar.
13
8.6. Matriz Constitutiva K de la Estructura.- Según lo indicado en el Apartado 4.3. y teniendo en cuenta la
(33), se tendrá

8.7. Matriz de Restricciones cinemáticas R.- Puesto que los nudos 1 y 2 son apoyos fijos, de acuerdo con el
Apartado 4.4. y la primera de (36) esta matriz tendrá la siguiente configuración

que, puesto que los elementos nulos han de ser matrices con dimensiones compatibles con las filas y co-
lumnas correspondientes, para Mathematica se transforma en13

donde, además, se incluye el cálculo de la matriz transpuesta de R.

8.8. Vectores de términos independientes.- Llamando vH, vB, vK y vR a los vectores asociados a las respecti-
vas matrices H, HT, K y R, se habrá de escribir:
- para el primero (vH), teniendo en cuenta los vectores de acción as de todos los nudos, se tendrá

- para el segundo y el tercero, puesto que de acuerdo con (23), las filas del vector V son 3e

- y para el último, puesto que son 4 los desplazamientos restringidos

8.9. Sistema Completo.- Con las matrices y vectores anteriores se tiene el sistema (41), cuya solución se
obtiene en la forma

con un Número de Condición de la matriz de coeficientes

de valor 3658.73, lo que indica que el sistema no está mal condicionado.

8.10. Resultados.- Extrayendo del vector solución s cada uno de los vectores incógnitas y deshaciendo el
cambio de variable para las magnitudes cinemáticas, se tiene

13
Con la notación N22 quiere expresarse que se trata de una matriz nula de dimensiones 2×2, aunque no se utilizan
subíndices. (Recuérdese que para Mathematica no se usan por comodidad de escritura).
14
expresiones que devuelven los valores de todas las incógnitas ivolucradas en el problema y que se muestran
en las Tablas siguientes:
Tabla 1.- Desplazamientos d de los nodos
u1 v1 u2 v2 u3 v3 u4 v4
0 0 0 0 609.581 4.163 663.704 104.083 ×10-6

Tabla 2.- Tensiones σ en los elementos

σ 1
x σ 1y τ1 σ x2 σ y2 τ2

1004.8 301.441 2.402 995.196 -1.201 -2.402

Tabla 3.- Deformaciones ε de los elementos

ε1x ε1y γ1 εx2 εy2 γ2

30.479 0 0.208 33.185 -9.992 -0.208 ×10-6

Tabla 4.- Reacciones (-λ)

X1 Y1 X2 Y2
-5000 -3002.4 -5000 3002.4

resultados que concuerdan con los del original, obviando las pequeñas diferencias de precisión.
Nota: si no se hubiese hecho el cambio de variable para eliminar el factor E, el Sistema Completo y su solución
vendrían dadas por

y al resolverlo, Mathematica emitiría el mensaje: ”LinearSolve:: luc: Result for LinearSolve of badly conditioned
matrix […] may contain significant numerical errors”, esto es, los resultados no son fiables debido a que el nuevo
Número de condición

pasa a valer 8.62×109, con lo que el sistema estárá muy mal condicionado. No obstante, en este ejemplo tan simple, la
solución que se obtiene con este sistema coincide con la calculada anteriormente.

9. Ilustrativo ejemplo 2
El ejercicio que se muestra en la Fig.4 está extraído de la Ref. 4§20-5.6, con los valores t = 1 in., E = 1 psi. y ν
= 0.25 en el que se pretende calcular los desplazamientos de los nudos y las tensiones de la placa, la cual
está empotrada en su lado inferior, con las cargas y cotas indicadas en la figura.

15
 1/2 1 1/2
4
7
j i i j
1 k
1 5 k
q=1
11

2
k k 6
2 i j j i
5 8
i j j i
2 3 k k 7
1

k 4 8 k
42 2 j i j
i
3 9
6
1 1

a) b)
Fig. 4.- Placa empotrada en la base.
a) Cotas y acciones sobre la placa. b) Modelización, discretización y Acciones Nodales de la mitad derecha.

9.1. Introducción de datos generales.- De acuerdo con el dibujo (n=9, e=8) y la expresión (31):

9.2. Coordenadas de los nodos del elemento prototipo y proyecciones de los lados.-Tomando como prototi-
po el elemento 1, las coordenadas de sus nodos respecto a un sistema de eje horizontal OX y vertical OY en
sentido levógiro con origen en el nodo i (nudo 1), las proyecciones de los lados del triángulo según (2) ven-
drán dadas por

9.3.- Matrices elementales de Volumen, de Simetrías y de Equilibrio.- Llamando Ae al área del elemento 1,
ve a su volumen y Ve a la matriz elemental de volumen dada en (23) y teniendo en cuenta las (50) se tendrá

con lo que las matrices de Equilibrio para los elementos 1 y 7 vendrán dadas, según (9), por

y puesto que los elementos 2 y 8 son doblemente simétricos al 1; el 3 y el 5 simétricos respecto al eje OY, y
el 4 y 6 simétricos respecto al eje OX, aplicando las expresiones (51), (52) y (53) se habrá de escribir, res-
pectivamente,

con lo que se obtienen las matrices de Equilibrio para todos los elementos del dominio a partir de las del
prototipo.

16
9.4. Matriz de Equilibrio H de la Estructura.- Partiendo de un arreglo de tantas filas como nudos y tantas
columnas como elementos, de acuerdo con la Fig.4 y el Apartado 4.1, el ensamblaje de la matriz H será

que para que Mathematica la lea y calcule su transpuesta, se ha de escribir

9.5. Matriz de Volumen V de la Estructura.- Según lo indicado en el Apartado 4.2., la matriz de Volumen
para toda la Estructura viene dada por un arreglo diagonal de dimensiones igual al de elementos, siendo los
términos de dicha diagonal las matrices elementales de Volumen de cada elemento que configura al domi-
nio ya calculadas en el Apartado 9.3. Así pues, puesto que el volumen de todos los elementos es el mismo,
tal matriz será

9.6. Matriz Constitutiva K de la Estructura.- Según lo indicado en el Apartado 4.3. y teniendo en cuenta la
(33), se tendrá, de forma similar a lo anterior

9.7. Matriz de Restricciones cinemáticas R.- Puesto que los nudos 1 y 2 son apoyos móviles con rodamien-
to vertical, mientras que el 3, 6 y 9 son fijos debido al empotramiento del lado inferior de la placa, de
acuerdo con el Apartado 4.4. y las expresiones (36), el ensamblaje de la matriz de Restricciones será

expresión que, puesto que los elementos nulos ha de ser matrices con dimensiones compatibles con las filas
y columnas correspondientes, para Mathematica se transforma en14

14
Con la notación n12 quiere expresarse el vector fila nulo de dimensión 1×2.
17
donde, además, se incluye el cálculo de la matriz transpuesta de R.
9.8. Vectores de términos independientes.- Llamando vH, vB, vK y vR a los vectores asociados a las respecti-
vas matrices H, HT, K y R, se habrá de escribir:
- para el primero (vH), teniendo en cuenta los vectores de acción as de todos los nudos, se tendrá

- para el segundo y el tercero, puesto que de acuerdo con (23), las filas del vector V son 3e

- y para el último, puesto que son 8 los desplazamientos restringidos

9.9. Sistema Completo.- Con las matrices y vectores anteriores se tiene el sistema (41), cuya solución se
calcula en la forma

con un Número de Condición de la matriz de coeficientes

de valor 15.2633, lo que indica que el sistema no está mal condicionado15.

9.10. Resultados.- Extrayendo del vector solución s cada uno de los vectores incógnitas y deshaciendo el
cambio de variable para las magnitudes cinemáticas, se tiene

expresiones que devuelven los valores de todas las incógnitas ivolucradas en el problema y que se muestran
en las Tablas siguientes:

15
Si no se hace el cambio de variable el condicionamiento sería el mismo, debido a que aquí se ha tomado E=1.
18
 Tabla 1.- Desplazamientos de los nudos
u1 v1 u2 v2 u3 v3 u4 v4 u5 v5 u6 v6 u7 v7 u8 v8 u9 v9
0 -1.892 0 -0.939 0 0 0.242 -1.933 0.147 -0.919 0 0 0.498 -2.047 0.361 -1.063 0 0

Tabla 2A.- Tensiones en los elementos 1, 2, 3 y 4

σ 1
x σ 1
y τ 1
σ x2 σ y2 τ2 σ x3 σ y3 τ3 σ x4 σ y4 τ4

-0.012 -0.017 0.022 -0.098 -0.978 0.008 -0.094 -0.962 0.008 -0.245 -0.98 0.059

Tabla 2B.- Tensiones en los elementos 5, 6, 7 y 8

σ 5
x σ 5
y τ 5
σ x6 σ y6 τ6 σ x7 σ y7 τ7 σ x8 σ y8 τ8

0.003 -1.013 -0.007 -0.034 -0.992 -0.003 -0.055 -1.077 0.009 -0.245 -0.98 0.059

Tabla 3A.- Deformaciones de los elementos 1, 2, 3 y 4

ε 1
x ε 1
y γ1
εx2 εy2 γ2 εx3 εy3 γ3 εx4 εy4 γ4

0.242 -1.014 0.055 0.147 -0.953 0.02 0.147 -0.939 0.02 0 -0.912 0.147

Tabla 3B.- Deformaciones de los elementos 5, 6, 7 y 8

ε 5
x ε 5
y γ5
εx6 εy6 γ6 εx7 εy7 γ7 εx8 εy8 γ8

0.256 -1.014 -0.019 0.214 -0.984 -0.006 0.214 -1.063 0.217 0 -0.919 0.147

Tabla 4.- Reacciones

X1 X2 X3 Y3 X6 Y6 X9 Y9
0.01 0.096 0.119 0.452 -0.059 0.98 -0.166 0.568

valores que concuerdan con los del original, obviando las pequeñas discrepancias de precisión. Como puede
comprobarse, las acciones y reacciones cumplen con las ecuaciones de equilibrio del modelo.

8. Conclusiones.- Vds. dirán.

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.
Será continuado…

Referencias:
[1] F. Martín Chica, “Método Mixto de la Solicitaciones y los Corrimientos”, (2014).
[2] 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).
[3] Daryl L. Logan, “A First Course in the Finite Element Method”, Fourth Editicion. University of Wisconsin-
Platteville, (2007).
[4] A. Ghali y Neville, “Análisis Estructural”, Diana Técnico, S.A. México, D.F. (1983).

19