Estudio de la flexin en vigas rectas con

inercia variable

Trabajo Fin de Mster

Departamento de Mecnica de Medios Continuos

y Teora de Estructuras

Juan Jos Osorno Gil

Supervisado por Prof. Juan Carlos Mosquera Feijo

Resumen
El presente trabajo estudia el comportamiento a flexin de vigas rectas de
seccin doble T con inercia variable en las que el canto vara con la longitud,
mientras que las dimensiones de las alas son constantes.
Se presenta un compendio de diversas formulaciones existentes sobre vigas
rectas de canto variable, identificando las hiptesis y simplificaciones que cada una
adopta. Posteriormente se hace un estudio comparativo de la respuesta segn cada
modelo analizado para el clculo de la matriz de rigidez de un elemento viga
no prismtico obtenida a partir de las funciones de forma.
Finalmente se valida una matriz de rigidez prctica, aproximada y sencilla para
su tratamiento en un cdigo de elementos finitos en el entorno de programacin
MATLAB[1] y se aplica a un anlisis numrico de naves 2-D compuestas por
elementos no prismticos.

Tabla de Contenido
1.

Introduccin ...................................................................................................... 4

2.

Revisin bibliogrfica ...................................................................................... 6

3.

2.1

Ecuaciones diferenciales de la elstica de una viga a flexin ......................... 6

2.2

Consideraciones del elemento viga ....................................................................... 9

2.3

Karabalis, D. L., Beskos, D. E. (1983) ................................................................ 12

2.4

Eisenberger, M. (1991) ........................................................................................... 16

2.5

Aristizbal-Ochoa, J. D. (1993) ............................................................................ 18

2.6

Al-Gahtani, H. J. (1996) ........................................................................................ 21

2.7

Saka, M. P. (1997) ................................................................................................... 23

2.8

Mtodo de la fuerza unidad con SAP2000 ......................................................... 26

Formulacin de las funciones de forma ..................................................... 29

3.1

4.

Funciones de forma................................................................................................. 29

Comparacin de resultados .......................................................................... 34

4.1

Aplicacin prctica .................................................................................................. 37

5.

Conclusiones .................................................................................................... 39

6.

Vas abiertas de investigacin ..................................................................... 40

7.

Referencias ...................................................................................................... 41

8.

Anexo ................................................................................................................. 43
8.1

Integracin de la ecuacin diferencial de la viga Euler-Bernoulli ............... 43

8.2

Solucin de la formulacin dbil .......................................................................... 45

8.3

Mtodo de Galerkin para deducir las ecuaciones de la viga .......................... 47

8.4

Cdigos de las formulaciones................................................................................ 50

1.

Introduccin

Los marcos con elementos no prismticos son preferidos en el diseo de

estructuras de acero en donde quiera que los requerimientos arquitectnicos
permitan su presencia. No solo proporcionan una mejor distribucin de tensiones
sino que tambin ofrecen un diseo ms liviano[2]. En la mayora de las
estructuras de ingeniera civil conformadas por vigas no prismticas el ancho de la
seccin transversal permanece constante mientras la altura vara lineal o no
linealmente (usualmente parablicamente) con la longitud.
Los elementos no prismticos son usados en muchas estructuras tales como
naves industriales, puentes y edificios de varias alturas. Debido a la facilidad de
construccin es muy prctico el uso de elementos con seccin doble T con canto
linealmente variable con la longitud[3]. Adems, la optimizacin del peso propio, el
incremento de la estabilidad, la flexibilidad en la fabricacin y el diseo e incluso
para satisfacer consideraciones arquitectnicas, son algunas ventajas que los
perfiles laminados no pueden ofrecer.
El uso de la va analtica para la resolucin de las condiciones de equilibrio de
las vigas con inercia variable es la nica manera de conocer cmo se comporta la
viga frente a flexin como una unidad estructural. Los mtodos numricos de
clculo discretizan la viga y asocian las condiciones de equilibrio a rebanadas, de
tal manera que los resultados obtenidos, aun siendo vlidos, no representan el
comportamiento global de la viga. En definitiva, la va numrica solo proporciona
soluciones cuantitativas, mientras que la va analtica permite adems obtener una
visin cualitativa del problema.
Muchos investigadores han abordado el problema de la flexin en vigas con
inercia variable sometidas a diferentes condiciones de contorno. Se han
desarrollado varios mtodos incluyendo soluciones de forma cerrada y tcnicas
numricas. Karabalis y Beskos[4] desarrollaron un mtodo basado en matrices de
rigidez y masa para vigas de acho constante y canto variable linealmente.
Eisenberger[5] desarroll una matriz de rigidez para algunos casos particulares
de anlisis de flexin de vigas no prismticas. Al-Gahtani[6] us el concepto de la
integral de contorno para encontrar la deformada, la distribucin de esfuerzos
cortantes y momentos en vigas no prismticas con condiciones cualesquiera de
contorno en ambos extremos. Saka[2], usando la matriz de rigidez obtenida por
Just (1977), desarroll un algoritmo basado en el mtodo del criterio de ptimo
para obtener el diseo ptimo de prticos metlicos con elementos no prismticos.
El presente trabajo estudia el comportamiento a flexin de vigas rectas con
inercia variable en las que solo el canto vara con la longitud. Se presenta un
compendio de varias formulaciones existentes para el clculo de la matriz de
rigidez elstica de vigas con canto variable y se hace un estudio comparativo de la
respuesta de cada modelo identificando las hiptesis y simplificaciones de cada
una.
Finalmente se presenta la matriz de rigidez obtenida mediante el mtodo de la
fuerza unidad para un elemento viga no prismtica. Se utiliza tanto el software
SAP2000[7] como las funciones de forma. Las formulaciones fueron desarrolladas
4

mediante cdigos de programacin utilizando el entorno MATLAB[1]. En el

captulo de resultados se presentan las matrices de rigidez obtenidas por las
diferentes formulaciones para un mismo elemento viga no prismtico. En los
Anexos se incluyen los cdigos de programacin utilizados.

2.

Revisin bibliogrfica

Una viga plana de directriz recta cargada por fuerzas transversales, se deforma
y adopta una configuracin llamada la elstica o curva de la deflexin de la viga[8].
En este captulo se determina la ecuacin de la curva de la viga sometida a flexin.
Posteriormente se presentan las formulaciones analizadas que resuelven, a partir
de las ecuaciones diferenciales de curva de deflexin, la matriz de rigidez a flexin
y a carga axial de un elemento viga no prismtico.

2.1 Ecuaciones diferenciales de la elstica de una viga a

flexin
La mayora de los procedimientos para encontrar deflexiones de vigas estn
basados en las ecuaciones diferenciales de la elstica y sus relaciones asociadas.
Consecuentemente se empezar por deducir la ecuacin bsica de la curva de
deflexin de una viga.

(a)

(b)
Figura1. Curva de deflexin de una viga en voladizo

Se considera una viga en voladizo con una carga concentrada actuando hacia
arriba en el extremo libre como se muestra en la Figura1(a). Bajo la accin de la
carga, la directriz de la viga se deforma en una curva, como se muestra en la
Figura1(b). Los ejes de referencia tienen su origen en el extremo fijo de la viga, con
el eje
en direccin a la derecha y el eje  en direccin hacia arriba.
Se asume que el plano
 es un plano de simetra de la viga y que todas las
cargas actan ste plano.

La deflexin  es el desplazamiento en la direccin  de cualquier punto del eje

de la viga. Para obtener la ecuacin de la curva de deflexin, se expresa la deflexin
 en funcin de la coordenada
. La deflexin  en cualquier punto  de la curva
se muestra en la Figura2(a). El punto  es localizado a una distancia
del origen
y el punto  a una distancia 
del primer punto. La deflexin en el segundo
punto es  + .
6

En la viga flectada se produce adems una rotacin en cada punto, , que es el

ngulo entre el eje
y la tangente a la curva deformada, como se muestra en la
Figura2(b). El ngulo de rotacin en el punto  es +  .

(a)

(b)
La relacin entre 
y  est dada por

Figura2. Elemento deformado de una viga a flexin.

 = 

en la que es el radio de curvatura, , quese expresa por la ecuacin

=

1 
=
 

(1)

Por convenio de signos, la curvatura es positiva cuando el ngulo de rotacin

aumenta en sentido antihorario segn se avanza en la direccin positiva del eje x.
La pendiente de la curva de deflexin es la primera derivada /
de la
expresin para la deflexin . Es decir, es el incremento  en la deflexin dividido
entre el incremento 
en la distancia a lo largo del eje
. Debido a que  y 
son
infinitesimalmente pequeos, la pendiente /
es igual a la tangente del ngulo
de rotacin , en consecuencia:
(
)
= tan ,


= arctan

(
)


(2)

Vigas con ngulos infinitesimales

Dado que la mayor parte de las vigas y columnas de las estructuras presentan
pequeos cambios en la forma durante la vida de servicio, los ngulos de rotacin,
las deflexiones y las curvaturas son muy pequeos[8]. Por lo tanto se pueden hacer
algunas aproximaciones matemticas:



As, la curvatura (1) resulta:

(3)

1 
=
 

(4)

Adems, cuando es pequeo,tan ( ) , luego de (2):

tan =

Al derivar (5) con respecto a
se obtiene:

  (
)
=








(5)

(6)

Combinando (4) con (6), la ecuacin de la curvatura queda:

=

1  (
)
=




(7)

Si el material de la viga es hookeano, linealmente elstico, la curvatura es:

1 (
)
=
 (
)

(8)

 (
) (
)
=


(
)

(9)

En la expresin anterior, M es el momento flector y EI es la rigidez a flexin de

la seccin transversal de la viga. Al combinar (7) con (8) se obtiene la ecuacin
diferencial bsica de la curva de deflexin de una viga:

A partir de las relaciones entre el momento flector (
), el esfuerzo cortante

(
) y la intensidad (
) de la carga distribuida, se obtiene:

= (
)


(
)
= (
)


(10)
(11)

Al derivar ambos lados de la ecuacin obtenemos:


 

(
)  ! =
=






8

(12)


 

(
)
= (
)
!=








La ecuacin de la deflexin de la viga de inercia variable sometida a flexin

podr obtenerse a partir de la resolucin analtica o numrica de cualquiera de las
ecuaciones (9) y (12).

Rigidez Axial
Se considera una barra sometida a fuerzas de traccin T, como se muestra en la
Figura 3.

Figura 3. Elemento barra sometido a fuerzas de traccin T.

De la ley de Hooke unidimensional se tiene

"# = $#

(13)

y de la relacin deformacin-movimiento:
$# =

%(
)


(14)

Del equilibrio de fuerzas para cargas aplicadas solo en los extremos, se tiene
&"# = ' = ()*
+,*+-

(15)

Sustituyendo (14) en (13), reemplazando en (15) y derivando con respecto a
se

obtiene la ecuacin diferencial que gobierna el comportamiento elstico-lineal del
elemento barra de la Figura 3:

%(
)
.&(
)
/=0




En la expresin anterior, % es la funcin del desplazamiento longitudinal en la

direccin
cualquier punto del elemento y &(
) es el rea de la seccin transversal,
la cual se considera variable a lo largo del eje
.

2.2

Consideraciones del elemento viga

Para estudiar la respuesta a flexin de vigas rectas con inercia variable

mediante diversas formulaciones se considerar un elemento viga de acero con
seccin en doble T doblemente simtrica en el que solo el canto vara con la
longitud, el ancho de las alas y los espesores del alma y las alas se consideran
constantes, como se muestra en la Figura 4.

Figura 4. Elemento viga no prismtico.

Propiedades geomtricas del elemento viga

Para los ejemplos que siguen en el estudio comparativo, se consideran los
siguientes valores de los parmetros geomtricos:
+1 = 0.02 .
45 = 0.30 .

Espesor del alma.

 = 0.60 .

Altura del alma en el extremo inicial.

Ancho de la seccin.

+5 = 0.02 .

 = 0.40 .

Espesor de las alas.

:< = 0.64 .

Altura total de la seccin en el extremo inicial.

Altura del alma en el extremo final.

:; = 0.44 .
= = 5 .

Altura total de la seccin en el extremo final.

 = 3 10@ A/

B=0

Longitud total del elemento.

Mdulo de elasticidad del acero.
Angulo del elemento respecto del eje horizontal.

Expresin del momento de inercia

Para obtener una relacin del momento de inercia y del rea con la longitud del
elemento viga no prismtico se considera una seccin intermedia como se muestra
en la Figura 5.

Figura 5. Seccin genrica del elemento viga no prismtico.

Canto y rea de la seccin:

:# = :; +

C:< :; D
=

&# = :# +1 + 2C45 +1 D+5

10

(16)

Momento esttico respecto de la fibra inferior de la seccin:

E &; F; = 45 +5 .:#

+5
+5
:#
/ + +1 C:# 2+5 D + 45 +5
2
2
2

(17)

Profundidad de la fibra baricntrica:

FG =

&; F;
&

(18)

Momento de inercia respecto del eje de flexin:

# =

I245 +5J + +1 C:# 2+5 D K

J

12

+ 45 +5 .:#


+5
+ 45 +5 . FG /
2



+5
:#
FG / + +1 C:# 2+5 D . FG /
2
2

(19)

Combinando las ecuaciones (16) hasta (19) se escriben las expresiones del rea y
del momento de inercia en funcin de
,
0, =
&(
) = 245 +5 + +1 

(  )
!
=

(20)

(
) = (+1 =J J + 645 =J  +5 + 1245 =J  +5 + 845 =J +5J 3+1 = J
+ 3+1 =  

1245 =  +5
+ 1245 =   +5
1245 =  +5
+ 1245 =  +5

+ 3+1 =J
 6+1 = 
 + 645 = +5
 + 3+1 = 

1245 =  +5
 + 645 = +5
 +1 J
J + 3+1 
J 3+1  
J
+ +1 J
J )/12=J

(21)

Los valores de &(
) y de (
) en los extremos O y P son:

En el extremo O:

En el extremo P:

& = 0,02  ,

& = 0,024  ,

 = 0,00063627 R
 = 0,001514 R

Formulaciones analizadas
Se han analizado en este trabajo diversas formulaciones para la obtencin de la
matriz de rigidez del elemento viga no prismtico. Se describen a continuacin.
Para todos los casos se obtuvo la matriz de rigidez del elemento de la Figura 4.
La primera formulacin es la presentada por Karabalis y Beskos (1983)[4] en la
cual se present un mtodo numrico para el anlisis esttico, dinmico y de
estabilidad de estructuras planas compuestas de vigas con canto variable. A partir
de la ecuacin diferencial de una viga a sometida a flexin obtuvo las funciones de
desplazamiento para la construccin de las matrices de rigidez de los elementos no
prismticos siguiendo el mtodo de los elementos finitos.
11

En 1991, Eisenberger[5] present las matrices de rigidez de elementos no

prismticos comunes incluyendo el efecto del esfuerzo cortante. Las rigideces
fueron formuladas a partir de los coeficientes de flexibilidad del elemento. Los
trminos de la matriz de flexibilidad fueron obtenidos usando el mtodo de la carga
unidad y fueron presentados para los dos elementos no prismticos ms usados:
con variacin lineal y con variacin parablica del canto.
En 1993 Aristizabal-Ochoa[9] propuso un algoritmo para evaluar la respuesta
esttica, de estabilidad y de vibracin de vigas y columnas no-prismticas.
Aplicando el principio de la viga conjugada dedujo los coeficientes bsicos que
componen la matriz de flexibilidad, la cual, una vez invertida, da lugar a los
coeficientes de rigidez a flexin a partir de los cuales se obtienen todos los
elementos de la matriz de rigidez.
Al-Gahtani (1996)[6] propuso un mtodo para obtener las expresiones cerradas
para los componentes de la matriz de rigidez y fuerzas y momentos de
empotramiento perfecto para elementos no prismticos. Utilizando las ecuaciones
diferenciales gobernantes de la teora de la viga Bernoulli-Euler y el mtodo de la
integral de contorno obtuvo las expresiones para las rigideces axial, torsional y de
flexin.
M. P. Saka (1997)[2] present un algoritmo para el diseo ptimo de marcos
metlicos compuestos por elementos prismticos y/o no prismticos. Utiliz la
matriz de rigidez para el elemento no prismtico obtenida por Just (1977) y luego
us el criterio de optimalidad para relacionar las variables de diseo (el canto)
sometidas a las restricciones de desplazamiento y de esfuerzos.
S. Z. Al-Sadder y H. Y. Qasrawi (2004)[10] presentaron una solucin analtica y
una matriz de rigidez para cualquier elemento viga-columna no prismtico con
conexiones semirrgidas en las uniones sometido a una fuerza axial de compresin
o tensin y a una carga generalizada. Para resolver la ecuacin diferencial
ordinaria de cuarto orden con coeficientes variables utilizaron una aproximacin
por series de potencias. No se tiene constancia de los resultados obtenidos ni de su
validez.

2.3

Karabalis, D. L., Beskos, D. E. (1983)

En [4] se presenta una metodologa de elementos finitos para el anlisis

esttico, de vibracin a flexin libre y de estabilidad de estructuras planas
linealmente elsticas conformadas por vigas no prismticas. Estas se asumen con
seccin transversal arbitraria pero un eje vertical de simetra, el ancho constante
mientras que el canto vara con la longitud.
El mtodo se apoya en los conceptos de rigidez a flexin, rigidez a esfuerzo axial,
rigidez geomtrica y matrices de masa consistentes para un elemento viga de ancho
constante y canto linealmente variable. Se obtienen expresiones generales para las
matrices de rigidez a flexin y a esfuerzo axial a partir de funciones de
desplazamiento, las cuales son soluciones exactas de las ecuaciones diferenciales
pertinentes. Tambin se obtienen expresiones explcitas para una matriz de rigidez
a flexin aproximada de un elemento no prismtico, cuyo canto vara linealmente,
12

con seccin en doble T o en cajn, obtenidas a partir de una funcin cbica de

desplazamiento.
Construccin de las matrices del elemento.

(a)

(b)
Figura 6. Geometra y convencin de signos de un elemento viga no prismtico genrico.

Se considera un elemento viga genrico de plano medio, no prismtico, de

longitud = hecho de material homogneo isotrpico y elstico lineal de mdulo de
elasticidad , como el que se muestra en la Figura 6(a). Tanto el rea de su seccin
transversal &(
) como el momento de Inercia S = (
) varan a lo largo de la
coordenada longitudinal
. En el supuesto de que la carga es aplicada solo en los
nudos 1 y 2 y que la deformacin de la viga est descrita mediante la teora de
Bernoulli-Euler, la ecuacin que gobierna la flexin es:

y la deformacin axial:


 (
)
(
)
!=0




&(
)

(22)

%(
)
=0


(23)

En las expresiones anteriores, (
) es el desplazamiento vertical y %(
) el

desplazamiento axial del centroide de una seccin genrica de la viga. La matriz de
rigidez TUV se calcula empleando el MEF.
Rigidez a flexin
Siguiendo la convencin de signos de la Figura 6 e Integrando la ecuacin (22)
cuatro veces se obtiene:
# S

# S

1
Y
(
) = W + W
+ WJ X X
Y [ + WR X X
Y [
(Y)
(Y)
Z Z

Z Z

W , W , WJ  WR son constantes de integracin. La evaluacin de la deflexin

(
) y la pendiente (
) = (
)/
en los puntos nodales 1 y 2 resulta en una
13

expresin para el campo de desplazamientos (
) en

desplazamientos nodales y las rotaciones  ,  ,  ,  como:

trminos

de

los

(
) = \] ] ]J ]R ^\ ,  ,  ,  ^_ = \]^\`^

El superndice ' representa transposicin, \`^ es el vector de

desplazamientos nodales y T]V_ es el vector de las funciones de forma ]; (O = 1,2,3,4)
dados explcitamente como:

b
a(
) + c(
)
:
:
= cd
ad =b
] =
+
a(
) +
c(
)
:
:

b
]J = a(
) c(
)
:
:
cd
ad
]R =
a(
) c(
)
:
:
] = 1

siendo:

1
b=X

,
(
)
d S

b=X

1

[ ,
ad = X X
(
)



(
)

d S

cd = X X

Z Z

Z Z

# S

# S


[
(
)

1
Y
a(
) = X X
Y [ , c(
) = X X
Y [
(Y)
(Y)
Z Z

: = ad cd b

Z Z

Siguiendo el procedimiento estndar del mtodo de los elementos finitos se

puede calcular la matriz de rigidez T
V a partir de la siguiente expresin
d

T
V = X\]^_ \]^(
)

Z

Donde las primas indican diferenciacin respecto de
. En consecuencia, de

acuerdo con la convencin de signos de la Figura 6(b), se escribe la relacin general
nodal fuerza-desplazamiento para la flexin del elemento viga no prismtica de la
forma:

en la que


j 


f  g = i 


i
J


h
R





J

R

J

J

JJ

RJ

R

m

R
l 
f g

JR
l 

RR
k 


=
JJ
=
J
=
J
=

14

b
:


:
ad
=
:

=

=
J
=
J
=



R
=
R
=
RJ
=
JR


=

cd
:

R
=
R
=

(= cd )
,
:

RR
=

(=ad cd )
:

Rigidez axial
La integracin de la ecuacin (23) da como resultado
#

%(
) = ( + ( X

1
Y
&(Y)

en la cual ( y c son constantes. La evaluacin de los desplazamientos axiales

%(
) en los nudos 1 y 2 permite expresar %(
) en trminos de los desplazamientos
nodales % y % de la Figura 6(b) como
%(
) = \]o ]p ^\% % ^_ = q] s q`s
_

]p =

P(
)
Pd

(24)

v
v
qus es el vector de desplazamientos nodales axiales y qfs es el vector de las
funciones de forma fo y fp dados explcitamente como

siendo

P(
)
]o = 1 x
y,
Pd
#

1
P(
) = X
Y ,
&(Y)
Z

Pd = X
Z

1


&(
)

Siguiendo el procedimiento estndar del mtodo de los elementos finitos se

puede calcular la matriz de rigidez axial T
V a partir de la expresin
T
V

= X\]^_ \]^&(
)

Z

(25)

En consecuencia, utilizando las ecuaciones (24) - (25) y la convencin de signos

de la Figura 6(b) se puede escribir la relacin nodal fuerza-desplazamiento para la
deformacin axial del elemento no prismtico tipo viga general de la forma

siendo

uuuu

uuuu
{
 %
y }% ~
z  | = x 

uuuu
uuuu
{






uuuu
uuuu

uuuu
 =
 =
 =

15


Pd

2.4

Eisenberger, M. (1991)

En [5] se presentan los trminos exactos de la matriz de rigidez para elementos

no prismticos incluyendo las deformaciones por cortante.

Deduccin de la matriz de rigidez

Se obtiene la matriz de rigidez invirtiendo la de flexibilidad correspondiente.
Para el elemento de la Figura 7 la matriz de flexibilidad de los desplazamientos en
el extremo del elemento es:
RR
 = 0
0

0
oo
po

0
op
pp

Los trminos ;< se obtienen mediante el mtodo de la carga unidad. Invirtiendo

esta matriz, se obtiene la rigidez correspondiente a los grados de libertad 4 a 6
segn la convencin de signos de la Figura 7. La matriz de rigidez ser:
Figura 7. Elemento genrico tipo viga de un entramado plano.

1
,
RR

Y los trminos son:

con

&=

& 0
= 0
0 W

pp
,
c

0
W
:

W=

pp
,
c

:=


c = oo pp op

oo
,
c

La matriz de rigidez del elemento , se obtiene imponiendo el equilibrio en las

configuraciones deformadas que corresponden a los desplazamientos unitarios
asociados al mtodo de la rigidez. Usando los trminos en , resulta:
j
i
= i
i
i
h

&

O

0
W + =
: + 2W= + =

&
0
0

0 W =
&
0

en la que = es la longitud del elemento viga.

16

0
W m
l
: W=l
0
l
l
W
k
:

Coeficientes de flexibilidad del elemento no prismtico

A continuacin se dan los trminos explcitos de la matriz de flexibilidad para el
elemento no prismtico cuyo canto vara linealmente. Se incluye el efecto del
esfuerzo cortante para lo cual se considera el factor de forma para cortante ], que
para secciones doble T por lo general est en el rango de 1.1 a 1.2[8].
Para aplicar este mtodo al elemento viga no prismtico con seccin transversal
doble T, se supone como un elemento equivalente con seccin transversal
rectangular de ancho constante y canto variable linealmente, conservando las
mismas reas y momentos de inercia en los extremos inicial y final de la viga de la
Figura 4.

Se denominan < y el canto de la seccin en los extremos P y del elemento,

respectivamente. El ancho de la viga es 4 y = es la longitud.

El resultado es una viga como se muestra en la Figura 8 con las siguientes

dimensiones
En el extremo :

4 = 0.032 ,

En el extremo P:

4 = 0.032 ,

= 0.617 ,

< = 0.828 ,

& = 0.024    = 0.001514 

&< = 0.02   < = 0.000634 

Figura 8. Elemento no prismtico con canto linealmente variable.

Los trminos de la matriz de flexibilidad son:

ln I K
=

RR =
+
4 <

oo =

siendo

<
2

x(1
+
b)
ln
+

1.5y
!
J

<
2<
4C< D
12=J

op =

6=
,
4<

pp =

6=C< + D
4< 

]C< D
b=
12a=

17

2.5

Aristizbal-Ochoa, J. D. (1993)

En [9] se presenta un mtodo basado en el principio clsico de la viga conjugada

por el cual todas las matrices caractersticas (de rigidez, de masa, de
amortiguamiento, geomtrica y elstica incremental as como el vector de cargas
nodales equivalentes) pueden ser obtenidos a partir de los coeficientes de rigidez
bsicos de la viga (cuatro coeficientes para un caso bidimensional (2-D) u ocho para
un caso tridimensional (3-D) sin usar polinomios aproximados de interpolacin.
Tambin se muestra que el anlisis de estructuras porticadas, con elementos
prismticos y no prismticos bajo cualquier condicin de carga y sustentacin,
puede ser realizado con la ayuda del mtodo de la viga conjugada y un adecuado
mtodo de integracin (Cuadratura Gaussiana, por ejemplo).

Modelo analtico

Figura 9. Modelo 2-D de viga propuesto por Aristizbal. Numeracin de grados de libertad.

Se considera un elemento viga 2-D que conecta los puntos A y B, como se

muestra en la Figura 9, con seis grados de libertad. Se asume que:

la viga est hecha de un material homogneo elstico lineal con un mdulo

de elasticidad ;

la seccin transversal de la viga puede variar de cualquier manera siempre y

cuando la directriz sea una lnea recta;
la viga est cargada en un plano a lo largo de uno de los ejes principales de
la seccin transversal con momentos de inercia  e  y las reas de seccin
transversal & y & en los extremos A y B, respectivamente;
todas las deformaciones son pequeas de modo que se puede aplicar
superposicin lineal.

Los grados de libertad y la convencin de signos de la viga AB se muestran en la

Figura 10(b). La notacin de las acciones nodales y de deformacin se muestra en
la Figura 10(a).

18

Figura 10. (a) Momentos, cortantes y esfuerzos axiles en los extremos A y B. (b) Rotaciones,
desplazamientos transversales y axiales en los extremos A y B.

Matriz de rigidez bsica

Cuando la viga AB est sometida a momentos en los extremos  y  los

extremos de la viga A y B experimentan sendos giros y , respectivamente
(Figura 11(a)). Aplicando el principio de la viga conjugada a la viga AB (Figura
11(b)) se obtienen las siguientes relaciones
=

(]  + ]  )=
,


(]  + ]  )=


(26)

que pueden reescribirse en forma matricial como:

= ]

=

 ]

] 

] 

(27)

Figura 11. (a) Momentos flectores aplicados en los extremos A y B. (b) Diagrama de curvaturas
correspondiente.

] , ] y ] son los coeficientes bsicos de flexibilidad a flexin de la viga AB
(normalizados con respecto a  /= ). Se pueden determinar analtica o
numricamente a partir de las integrales que los definen:
19

 (=
)


] = J X
=
#

]


(=
)
= JX


=
#



] = J X 

=
#

donde # es el momento de inercia de la seccin transversal a una distancia x

desde el extremo A.

1
]

=
 ]

= ] ] ]


Invirtiendo la matriz de flexibilidad de (27) se obtiene la matriz de rigidez:

]

]

Por tanto, los coeficientes de rigidez a la flexin de la viga AB son:

 =

]


= ] ] ]
 =

 =

]


= ] ] ]


]

= ] ] ]

Para los seis grados de libertad considerados, la matriz de rigidez de la viga AB

queda:
oo
j
i 0
i
i
i 0
TUV = i
op
i
i 0
i
i
h 0

0
 +  + 2 
=
 + 
=
0
(  +  + 2  )
=
 + 
=

0
 + 
=

op





0
0

0
oo
(  +  )
0
=

0
(  +  + 2  )
=
(  +  )
=
0
 +  + 2 
=
(  +  )
=

0
m
 + 
l
=
l
l

l
l
0
(  +  )l
l
=
l
l

k

Donde oo y op son las rigideces axiales las cuales se obtienen mediante

oo = op = X

1


&#

(28)

&# representa el rea de seccin transversal de la viga AB a una distancia x del

extremo A(Figura 9). La ecuacin (28) es genrica para determinar la rigidez axial
de elementos no prismticos tipo barra. Tambin se observa que los componentes
correspondientes a los grados de libertad de traslacin son obtenidos a partir del
equilibrio esttico y son expresados en trminos de los coeficientes de flexin
bsicos ,  y .

20

2.6

Al-Gahtani, H. J. (1996)

En [6] se presenta un procedimiento para deducir las expresiones de rigidez de

forma cerrada para elementos 3D no prismticos simtricos con directriz recta. El
procedimiento se basa en el mtodo de la integral de contorno, el cual no requiere
discretizacin del elemento. Se basa en la teora de la viga Bernoulli-Euler, la cual
es bastante adecuada siempre que la relacin de ahusamiento no sea muy grande.
Conduce a relaciones exactas entre fuerzas y desplazamientos en los extremos del
elemento.

Planteamiento del problema

Se considera una viga Bernoulli-Euler no uniforme de longitud = como se

muestra en la Figura 12. La ecuacin de rigidez esttica puede escribirse como:
 = U{

donde

U es la matriz de rigidez de tamao 12x12;

 es una matriz columna que contiene las fuerzas axiales (  y  ) ,

momentos torsores (R y Z ) , esfuerzos cortantes y momentos flectores
asociados a la flexin en el plano
 ( , @ , p y  ) , y los esfuerzos
cortantes y momentos flectores asociados a la flexin en el plano

[ (J ,  , o y  );

{ es una matriz columna que contiene los movimientos axiales, de torsin y

de flexin correspondientes.

Figura 12. Convencin de signos para la viga no prismtica

Se asume que:

la viga est hecha de un material homogneo, istropo elstico lineal con un

mdulo de elasticidad  y un mdulo cortante a;

el eje centroidal (eje-
) y las direcciones de los ejes principales (ejes  y [)

son las mismas para todas las secciones transversales;

la deformacin y el acoplamiento entre los comportamientos de torsin y de

flexin son insignificantes.

21

Adems, el rea &(
) y el momento de inercia (
) de la seccin transversal

estn dados por:
&(
) = &Z I1 +

(

K ;
=

(
) = Z I1 +

(

K
=

&Z y Z son propiedades de la seccin transversal en el origen; ( es una constante

que representa el grado de ahusamiento; y * es un nmero real que depende de la
forma de la seccin transversal.
Para obtener la matriz de rigidez, las ecuaciones diferenciales que gobiernan
los comportamientos axial y de flexin de la viga se consideran separadamente.

Rigidez axial
La ecuacin diferencial que rige el comportamiento de una barra cargada
axialmente est dada por:

%
&(
) + (
) = 0,
(0, =)




(29)

En la expresin anterior, % es el desplazamiento axial, y es la fuerza axial

distribuida. Las condiciones de contorno son:
%
(0) = 


%
(=) = 
%(=) = { ) &


%(0) = { ) &

Para obtener la llamada formulacin inversa se multiplican ambos lados de

(29) por una funcin de peso % y se integra dos veces a lo largo de la longitud de la
viga para obtener
d


% d
% d
%

X % &(
)

+ % &(
) %&(
)
+ X % (
) 
= 0





Z

Z

U
U

U

Operando se llega a la matriz de rigidez a efecto axial[6]:

donde

U=

U = U =

U = U

&Z
(
+
=
log(1 + ()

&Z
( (1 *)
+
T(1 + () 1V
=

Rigidez a flexin

U = U = U

,,

; U = U = U ,,

*=1

*>1

La ecuacin diferencial que gobierna la flexin en el plano
 est dada por

22



x(
)
y + (
) = 0,
(0, =)





(
) es la flecha y (
) es la carga vertical distribuida.

El giro (
), el momento (
) y el esfuerzo cortante (
) estn relacionados con
la variable primaria (
) de la forma
(
) =


,


(
) = (
)

Las condiciones en los extremos son:

= { , = {p ,  = p,  = 






(
) =

en
= 0;



x(
)  y




= {@ , = { ,  = ,  = @ en
= =.

Siguiendo el mismo procedimiento que para la matriz de rigidez axial, esto es,
multiplicando la ecuacin diferencial por una funcin de peso e integrando
cuatro veces a lo largo de la longitud de la viga, finalmente se llega a[6]:

U =

2.7

U@
Up@
U@@

U,
m
Up, l
U@, l
l
U, k

Z
Z
T4(1 + ()(3 + 3( + (  )V
; Up =  T2(1 + ()(3 + ()V
J
=
=
Z
Z
= U ,
U, =  T2(1 + () (3 + 2()V ; Upp =
T4(1 + ()V
=
=
Z
T2(1 + () (3 + 2()V ; U@@ = U
Up@ = Up ,
Up, =
=
Z
T4(1 + ()J V
U@, = U, ,
U, =
=

siendo

U@

U Up
j
Upp
i
U=i
i
h
O

Saka, M. P. (1997)

En [2] se plantea el diseo ptimo de marcos metlicos mediante el tratamiento

del canto de la seccin en cada nudo como una variable de diseo. Entre nudos
consecutivos, el canto se asume linealmente variable. El ancho y el espesor de las
alas y del alma se consideran constantes, como se muestra en la Figura 13.

23

Figura 13. Elemento viga no prismtico con seccin en doble T.

Matriz de rigidez del elemento no prismtico

Se considera un material homogneo, isotrpico y linealmente elstico de
mdulo de elasticidad E. Los elementos de la matriz de rigidez al esfuerzo axial y a
flexin del elemento no prismtico se basan en las funciones de desplazamiento, las
cuales son las soluciones exactas de las ecuaciones diferenciales gobernantes
pertinentes[2].

Figura 14. Fuerzas y deformaciones en los extremos del elemento no prismtico.

La matriz de rigidez del elemento no prismtico, mostrado en la Figura 14, se

obtiene mediante el mtodo de los elementos finitos. La matriz tiene la siguiente
forma en el sistema de coordenadas locales:

donde

U
j 
i 0
i 0
UC:; , :< D = i
U
i R
i 0
h 0

U
UJ
0
Uo
Up

UJJ
0
UoJ
UpJ

URR
0
0

U = URR = R =

#
Z

O
Uoo
Upo

m
l
l
l
l
l
Upp k



U =  TWJ
b + 2WJ WR  + WR
 V

UJ = TWJ WJ b + WJ WR  + WR WJ  + WR WR  V
24



b + 2 WJ WR  + WR
 V
UJJ = TWJ

Uo = TWJJ WJ b + WJJ WR  + WRJ WJ  + WRJ WR  V

UoJ = TWJJ WJ b + WJJ WR  + WRJ WJ  + WRJ WR  V


Uoo = TWJJ
b + 2WJJ WRJ  + WRJ
 V

Up = TWJR WJ b + WJR WR  + WRR WJ  + WRR WR  V

UpJ = TWJR WJ b + WJR WR  + WRR WJ  + WRR WR  V

Upo =  TWJR WJJ b + WJR WRJ  + WRR WJJ  + WRR WRJ  V


Upp = TWJR
b + 2WJR WRR  + WRR
 V

E es el mdulo de elasticidad; los dems trminos son:

b = b b ,

a = a a ,

siendo

WJ

 =   ,


,
:

WR =

b
,
:

b a
,
WJJ = WJ ,
WRJ = WR ,
:
a b
=
,
: = a  c b + b ;
:

WR =

WRR

 =   ,

b = X

#Z

c = X X
#Z



,


b = X

a = X X




,


#Z



,


 = X

#Z




,


c = X X
#




;


c 
:
 c
WJR =
:

WJ =

b = b  b 



,


a = X X
#

c = c c

 = X




,


 = X

#Z

 

,






El rea del elemento prismtico en la seccin
se expresa as:

siendo:

&# =  + 

 = +1 C:; 2+5 D + 245 +5 ;  = +1

El momento de inercia se expresa como:

siendo

 = :; +5 ,

 = ( + 
) (J + R
)

 =

:< :;
,

J =

25

 X

:< :;

+1  45 +5
+
,
12
2

R =

+1 
12

 



2.8

Mtodo de la fuerza unidad con SAP2000

SAP2000[7] admite una variacin del momento de inercia para los elementos
no prismticos de la siguiente manera:

siendo

S (
) =  I1 K +  I K
=
=

(30)

* = 1 para variacin lineal;

* = 2 para variacin parablica;

* = 3 para variacin cbica.

Para perfiles en doble T con variacin lineal del canto, la ley de variaciones es
parablica de segundo grado.

 y  , son los momentos de inercia en los extremos inicial y final del elemento
viga, respectivamente.
Desplazamiento unitario

El significado fsico de la matriz de rigidez puede ponerse de manifiesto

imponiendo a la estructura el siguiente estado de deformacin: al grado de libertad
P se le impone una deformacin de valor unidad < = 1 mientras que los restantes
grados de libertad se mantienen fijos ; = 0, O P [11]. Para imponer este estado de
deformacin hay que aplicar unas fuerzas sobre los distintos grados de libertad,
cuyos valores se trata de determinar.
K
F
j

i
F
=i

i

F  hK 

K

K

K

K  = 0
K
m

l
l = 1 = K
l

Knnk = 0 K

La ecuacin de equilibrio en la situacin de deformacin impuesta es:

Se observa que las fuerzas a efectuar sobre los grados de libertad corresponden a
la columna P de la matriz de rigidez. Por lo tanto se puede decir que la columna P de
la matriz de rigidez son las fuerzas que es necesario aplicar sobre los grados de
libertad de la estructura, para imponer un estado de deformacin compuesto por un
desplazamiento unidad en el grado de libertad P, y cero en todos los dems.
K
F


K = F


K F 

Las fuerzas que hay que aplicar sobre los distintos grados de libertad son:

26

Sobre el grado de libertad P hay que aplicar una fuerza exterior de valor
< = U<< , que es precisamente la fuerza necesaria para producir el < = 1.
Esta fuerza corresponde al trmino de la diagonal de K.
En todos los dems grados de libertad aparecen unas fuerzas de reaccin,
precisamente para que no se muevan, de valor ; = U;< .

Por lo tanto, un trmino cualquiera U;< de la matriz K representa la fuerza de

reaccin exterior que aparece en la direccin O, al aplicarse un desplazamiento
unitario en la direccin P (< = 1), siendo cero todos los dems desplazamientos
 = 0, O P) .

As, se calcularon con ayuda de SAP2000[7] los 6 estados de carga, asociados a

los grados de libertad buscados. El elemento de seccin variable se ha definido con
seccin creciente en sentido del eje x positivo.

(a)

(b)

Figura 15. Seccin de los extremos inicial (a) y final (b) viga no prismtica con seccin en doble T.

(a)

(b)
Figura 16. Perspectiva (a) y alzado (b) elemento viga no prismtica con seccin en doble T.

27

En el anlisis solo se consideran las deformaciones axial y de flexin. Las

deformaciones por cortante se ignoran. El anlisis se realiz con la opcin de
subdivisin automtica de la viga, en diez elementos de igual longitud y presenta
los resultados para el elemento entero.

Otras formulaciones estudiadas

En [10] Al-Sadder y Qasrawi (2004) presentaron una matriz de rigidez secante
que ellos mismos denominaron exacta y un vector de reacciones de empotramiento
perfecto para un elemento viga-columna no prismtico general con conexiones
semirrgidas sometido a cualquier fuerza axial.
El mtodo se desarrolla mediante una solucin analtica para resolver la
ecuacin diferencial ordinaria de cuarto orden con coeficientes variables de un
elemento viga-columna no prismtico usando el mtodo de las series de potencias.
Al final presentan los elementos de la matriz de rigidez axial y a flexin en funcin
de unas constantes U a UR ,  a o , F y F las cuales dicen son presentadas en otro
artculo suyo1. Pero dicho artculo no ha sido publicado desde entonces y por tanto
no se ha podido analizar su trabajo.

Al-Sadder, S.Z., Qasrawi, H.Y. (2004), Exact Secant Stiffness Matrix for Non-prismatic
Beam-Column Member, Submitted to International Journal of Computers and Structures.

28

3.

Formulacin de las funciones de forma

3.1

Funciones de forma

Matriz de rigidez basada en la teora de Euller-Bernouli

En lo que sigue se consideran solamente deformaciones por flexin.
La ecuacin diferencial que gobierna el comportamiento elstico-lineal de una
viga sometida a flexin est fundamentada en que las secciones transversales
planas y perpendiculares al eje centroidal longitudinal de la viga en el estado sin
deformar permanecen planas y perpendiculares al eje longitudinal despus de que
ocurra la flexin. Sin embargo es una suposicin razonable que conduce a
ecuaciones que predicen con bastante precisin el comportamiento para la mayora
de las vigas[12].
Considerando el esquema del elemento viga de la Figura
desplazamiento transversal en una seccin genrica
como:

17 se asume el

(
) = ,J
J + ,R
 + ,o
+ ,p

(31)

La funcin de desplazamiento cbica (31) satisface la ecuacin diferencial de la

viga; es apropiada porque hay cuatro grados de libertad (un desplazamiento
transversal ; y una pequea rotacin ; en cada nudo).

(a)

(b)

Figura 17 (a) Elemento viga con desplazamientos, rotaciones, fuerzas y momentos nodales positivos.
(b) Convencin de signos para esfuerzos cortantes y momentos flectores.

Se expresa  como una funcin de los grados de libertad nodales  ,  ,  y

 de la siguiente manera:
(0) =  = ,p ;

(0)
=  = ,o


(=) =  = ,J =J + ,R = + ,o = + ,p
(=)
=  = 3,J = + 2,R = + ,o


(32)

Donde = /
para la pequea rotacin asumida. Resolviendo las

ecuaciones (32) para ,J hasta ,p en trminos de los grados de libertad nodales y
sustituyendo en (31), se tiene:
2C  D  +  J
3C  D 2 +  
 = x
+
y

x
+
y
+ 
+ 
=J
=
=
=
29

(33)

En forma matricial, (33) se expresa como:

siendo

= TAV\^




\^ =
,



AJ =

1
(2
J 3
 = + =J ),
=J
1
Ao = J (2
J + 3
 =),
=

TAV = TAJ

AR

Ao

Ap V

1 J
(
= 2
 = +
=J )
=J
1
Ap = J (
J =
 = )
=

AR =

(34)

AJ , AR , Ao y Ap se denominan funciones de forma para un elemento viga,

conocidas como las funciones de interpolacin cbica de Hermite.
Deformaciones axiales

 , 
'

 , 

# , ]#

 , 

# , ]#
'

Figura 18. Elemento viga sometido a fuerzas de traccin . Desplazamientos nodales positivos y
fuerzas en ejes locales.

Denominando \
, ^ a los ejes locales, se asume que el campo del

desplazamientos segn
es lineal. Los valores en los extremos del elemento viga
son # y # . As, la funcin del desplazamiento se define como
% = , + ,

(35)

En general, el nmero total de coeficientes , es igual al nmero total de grados

de libertad asociados al elemento, que en este caso es 2: los desplazamientos
nodales. En forma matricial, la (35) se puede expresar como
,
% = T1
V }, ~


Para expresar % como una funcin de los desplazamientos nodales # y # se
resuelve para , y , ; de (35) se tiene
%(0) = # = ,

%(=) = # = , = + #

30

(36)

resolviendo de las ecuaciones (36) para , se llega a

, =

# #
=

(37)

sustituyendo (36) y (37) en (35) se obtiene

% =

# #
!
+ #
=

en forma matricial, (38) se expresa como

% = 1

siendo

=

#

= TA
= #

A = 1 ,
=

=

A =

(38)

A V #
#
(39)

llamadas las funciones de forma. Para un elemento, el desplazamiento en

cualquier punto interior se expresa en funcin de los desplazamientos nodales
mediante la ecuacin
\]^ = TAV\: ^

Donde \]^ es un vector de componentes de desplazamiento en cualquier punto;

es un vector de desplazamientos nodales; y TAV es la matriz de las funciones de
interpolacin o funciones de forma (Figura 19).

\: ^

Cada funcin de forma A; en las ecuaciones (34) y (39) representa la deformada

del elemento cuando :; = 1 mientras los otros desplazamientos nodales son iguales
a cero[13]. Las funciones de forma corresponden a las deformadas reales de un
elemento viga prismtico (ignorando la deformacin por cortante); se pueden usar
para obtener la matriz de rigidez de un elemento no prismtico.
A = 1

=

A =

=

AR =

J



2
+

=
=

Ap =

J


=
=

J


AJ = 2 I K 3 I K + 1
=
=

J


Ao = 2 I K + 3 I K
=
=

Figura 19. Funciones de forma del elemento viga de Euler-Bernoulli.

31

Relaciones deformacin-desplazamiento y tensin-deformacin

Se considera la relacin entre deformacin axial y desplazamiento
%(
, )


(40)

(
)


(41)

$# (
, ) =

donde % es el desplazamiento horizontal, que en la geometra deformada se

relaciona con el desplazamiento vertical mediante
% = 

sustituyendo (41) en (40) se obtiene

$# (
, ) = 

 



(42)

sustituyendo (34) y (39) en (40) y (42) se obtiene

T$# V =

T$# V = 

1
TV = =
0

siendo

 AJ



A


A
qs,


 AR



 Ao


0

\$# ^ = TVqs

 Ap
qs,


1
=

\$# ^ = TVqs
0

6 2

2 3

6 2

2 3

. 1/
. 2/ 0
. 1/
. 1/
= =
= =
= =
= =

La relacin esfuerzo-deformacin est dada por

donde

\"# ^ = T:V\$# ^

T:V =

&(
)
0

0
(
)

Finalmente la matriz de rigidez axial y a flexin del elemento en ejes locales

est dada por
d

U = XTV_ T:VTV 

Z

(43)

En la matriz de rigidez (43) se asume que la viga es esbelta; esto es, la relacin
de las dimensiones longitud, =, y canto, , es grande. En este caso la deflexin
debida a la flexin proporcionada mediante el uso de la matriz de rigidez de la
ecuacin (43) es bastante adecuada. Sin embargo para vigas cortas y de gran canto
la deformacin por esfuerzo cortante puede ser significativa y tener una
contribucin del mismo orden de magnitud a la deformacin total de la viga. Esto
puede observarse en las expresiones de las contribuciones de flexin y cortante a la
deflexin de la viga, donde la contribucin de la flexin es del orden de (=/)J
32

mientras que la contribucin del cortante es solo de orden de (/). Una regla
general para vigas de seccin transversal rectangular, es que para una longitud
igual a al menos ocho veces la altura de la seccin, la deflexin por cortante
transversal es menor que el 5% de la deflexin por flexin.

33

4.

Comparacin de resultados

A continuacin se presentan los resultados de la matriz de rigidez a flexin y a

esfuerzo axial del elemento viga no prismtico utilizando las diferentes
formulaciones mostradas en el apartado 3.
Para el estudio comparativo, se han considerado los siguientes valores de los
parmetros geomtricos:
+1 = 0.02 .
45 = 0.30 .

Espesor del alma.

 = 0.60 .

Altura del alma en el extremo inicial.

+5 = 0.02 .

Ancho de lasalas.

 = 0.40 .

Espesor de las alas.

:< = 0.64 .

Altura total de la seccin en el extremo inicial.

:; = 0.44 .
= = 5 .

 = 3 10@ A/

B=0

Altura del alma en el extremo final.

Altura total de la seccin en el extremo final.

Longitud total del elemento.
Mdulo de elasticidad del acero.
Angulo del elemento respecto del eje horizontal.

Karabalis, D. L., Beskos(1983)

13164
0
0
13164
0
0
j 0
288.84 618.49
0
288.84 825.7m
i
l
0
618.49
1910
0
618.49 1182.4 l

= i
0
0
13164
0
0 l
i13164
i 0
288.84 618.49
0
288.84
825.7 l
h 0
825.7
1182.4
0
825.7
2946.1 k

Eisenberger, M. (1991)
13773
j 0
i
0
= i
i13773
i 0
h 0

0
0
13773
262.57
560.58
0
560.58
1774.8
0
0
0
13773
262.57 560.58
0
752.29
1028.1
0

Aristizabal-Ochoa, J. D., (1993)

13164
j 0
i
0
U=i
13164
i
i 0
h 0

0
0
13164
288.84
618.49
0
618.49
1910
0
0
0
13164
288.84 618.49
0
825.7
1182.4
0
34

0
262.57
560.58
0
262.57
752.29

0
752.29 m
l
1028.1 l
0
l
752.29l
2733.3 k

0
288.84
618.49
0
288.84
825.7

0
825.7 m
l
1182.4 l
0 l
825.7l
2946.1 k

Al-Gahtani,H. J. (1996)
13444
j 0
i
0
U=i
13444
i
i 0
h 0

0
0
13444
291.09
620.36
0
620.36
1908.8
0
0
0
13444
291.09 620.36
0
835.1
1193
0

0
291.09
620.36
0
291.09
835.1

0
835.1 m
l
1193 l
0 l
835.1l
2982.5 k

Saka, M. P.(1997)

13164
0
0
13164
0
0
j 0
296
634
0
296 846 m
i
l
0
634 1959
0
634 1212 l
U (:; , :< ) = i
0
0
13164
0
0 l
i13164
i 0
296 634
0
296 846l
h 0
846 1212
0
846 3019 k

Fuerza unidad (SAP2000)

13163.56
0
0
13163.56
j
0
293.62 616.63
0
i
0
623.34
1914.44
0
U=i
13163.56
0
0
13163.56
i
i
0
285.15 625.09
0
h
0
824.22 1189.22
0

Funciones de forma de Hermite

13200
j 0
i
0
U = i
i13200
i 0
h 0

0
0
13200
303.24
652.99
0
652.99
1992.6
0
0
0
13200
303.24 652.99
0
863.23
1272.3
0

0
0
285.94 824.60m
l
617.58
1182.12 l
0
0
l
294.41
833.07 l
832.66
2961.41 k
0
303.24
652.99
0
303.24
863.23

0
863.23 m
l
1272.3 l
0
l
863.23l
3043.8 k

Las matrices de rigidez obtenidas presentan valores muy similares. A

continuacin se presentar la diferencia de los resultados respecto de la matriz
obtenida a partir de las funciones de forma.

+5

0.30

0.40 ;

:(
)

:; = 0.44

5.0

<

+1

Figura 20. Viga no prismtica con seccin en doble T.

35

+5
0.30

0.60

Para calcular la precisin de las formulaciones analizadas se realiz el clculo

de la deflexin vertical de la viga en voladizo mostrada en la Figura 20. La altura
de la seccin en el extremo libre permanece constante y la seccin apoyada vara
entre 0.64, 0.84 1.04, 1.24, 1.44, 1.64 1.84 y 2.04 m para radios de ahusamiento de
1.5, 2, 2.5, 3, 3.5, 4, 4.5 y 5 respectivamente.
Las formulaciones analizadas fueron Karabalis-Beskos, Aristizbal-Ochoa,
Eisenberger, Al-Gahtani, Saka, SAP2000 y la de las funciones de forma. Los
resultados se presentan en la Figura 21. Expresan la diferencia porcentual relativa
- = ` `/` 100% en el clculo de la deflexin vertical en el extremo libre debido
a una carga aplicada = 70 UA para diferentes grados de ahusamiento. El smbolo
` representa el valor exacto de la deflexin calculada mediante el mtodo de las
funciones de forma.
En la Figura 21 se observa que a excepcin de la formulacin de Eisenberger,
las dems presentan resultados similares para el clculo de la deflexin mxima
del extremo libre de la viga de la Figura 20. Esto se debe a que, adems de
considerar un factor de forma para el cortante, dicha formulacin solo considera
como parmetros geomtricos el ancho de la seccin y los cantos de los extremos del
elemento viga; por lo que no es posible considerar un rea y un momento de inercia
correspondiente a una seccin en doble T.
45,00
40,00

% Diferencia

35,00
Karabalis-Beskos/Aristizabal
Einseberger
Saka
SAP2000
Al-Gahtani

30,00
25,00
20,00
15,00
10,00
5,00
0,00
1,5

2,5

3 dj/di 3,5

4,5

Figura 21. Efecto del radio de ahusamiento en el error del clculo de la deflexin vertical en el
extremo libre de la viga mostrada en la Figura 20.

En la Figura 22 se presenta la diferencia en los resultados para las

formulaciones ms precisas. Se puede observar que para valores de ahusamiento
pequeos todas las formulaciones ofrecen resultados muy similares a los obtenidos
a partir de las funciones de forma. La formulacin presentada por Al-Gahtani
presenta resultados muy diferentes para valores de ahusamiento a partir de 2. Esto
puede deberse a que los coeficientes de la matriz de rigidez estn dados nicamente
en trminos del momento de inercia y el rea de la seccin transversal en el
extremo inicial y de una constante que relaciona el radio de ahusamiento del
elemento. Los datos de la seccin transversal del extremo final no son tenidos en
cuenta.

36

Las formulaciones de Karabalis-Beskos y Aristizabal presentan los resultados

ms parecidos a los obtenidos mediante las funciones de forma para radios de
ahusamiento pequeos pero la diferencia en los resultados es mayor a medida que
se aumenta el radio de ahusamiento del elemento. Los resultados hallados a travs
del programa SAP2000 son los ms parecidos a los obtenidos a partir de las
funciones de forma para todos los radios de ahusamiento.

14,00

KarabalisBeskos/Aristizabal
Saka

12,00

% Diferencia

10,00
8,00
6,00
4,00
2,00
0,00
1,5

2,5

3
dj/di

3,5

4,5

Figura 22. Efecto del radio de ahusamiento en el error del clculo de la deflexin vertical en el
extremo libre de la viga mostrada en la Figura 20 para las formulaciones ms precisas.

4.1

Aplicacin prctica

Con la matriz de rigidez obtenida por la formulacin Aristizbal-Ochoa se ha

elaborado un cdigo en el entorno MATLAB[1] para analizar naves 2D a dos aguas,
compuestas por elementos no prismticos de seccin en doble T en los que solo el
canto vara linealmente con la longitud y el ancho es constante.
A continuacin se presentan los resultados obtenidos para una nave sometida a
cargas puntuales verticales de valor 70 kN en los nudos de la cubierta y una carga
distribuida de 3.5 kN en direccin horizontal en los elementos 1 y 2, como se
muestra en la Figura 23.
En los resultados presentados en la Figura 24(a) se observa que para el clculo
de las deflexiones de los grados de libertad 7, 8 y 9, que corresponden al nodo ms
alto de la nave, se obtuvo una diferencia inferior al 15% para los desplazamientos,
las deflexiones y las rotaciones de los grados de libertad correspondientes a los tres
nodos del tejado. El cdigo de programacin se adjunta en un anexo.

37

Figura 23. Nave 2-D analizada mediante SAP2000

(a)

(b)

Figura 24. Resultados obtenidos para las deflexiones de la nave de la Figura 23 mediante (a)
MATLAB (b) SAP2000.

38

5.

Conclusiones

Del trabajo realizado para el presente documento se extraen las siguientes

conclusiones:

Se ha compilado un conjunto de formulaciones analticas y numricas del

comportamiento de un elemento no prismtico tipo viga. Se ha realizado un
anlisis comparativo de aquellas, atendiendo a las hiptesis y aplicabilidad
de cada uno.

Se ha seleccionado el mtodo de las funciones de forma como procedimiento

para deducir la matriz de rigidez del elemento no prismtico. Se ha
empleado dicha matriz para analizar, mediante el mtodo de los elementos
finitos, la respuesta de una nave a 2 aguas.

De todas las formulaciones analizadas las que presentan la menor diferencia

para el anlisis matricial de elementos no prismticos con seccin en doble
T, en los que solo el canto vara linealmente con la longitud y el ancho es
constante, son las propuestas por Karabalis-Beskos, Saka y AristizbalOchoa.

La formulacin presentada por Eisenberger no es fiable para calcular la

matriz de rigidez de elementos no prismticos con seccin en doble T, sobre
todo si el radio de ahusamiento es muy grande (a partir de   = 1,5).

As mismo, la formulacin presentada por Al-Saddery Qasrawi no puede

usarse dado que no se conocen los valores de las constantes U a UR ,  a o ,
F y F que anuncian en su publicacin.

39

6.

Vas abiertas de investigacin

Se plantean las siguientes lneas por el inters que suscitan para su aplicacin
profesional:

Estudio de vigas no prismticas trapeciales, en las que el canto vara

linealmente pero la cara superior se mantiene horizontal. Anlisis del efecto
de las tensiones tangenciales y de la variacin de las tensiones normales.

Obtencin de una matriz de rigidez de la viga anterior en la que adems se

tienen en cuenta las deformaciones por esfuerzo cortante.

Obtencin de una matriz de rigidez de la viga anterior en la que adems se

tiene en cuenta el efecto acoplado del esfuerzo axil.

40

7.

Referencias

1. MATLAB, Ver.R2011a. (2011). User's Manual. s.l. : The MathWorks, Inc, Natick, Mass, USA.
2. Saka, M. P. (1997). Optimum design of steel frames with tapered members. s.l. : Computers
and Structures Vol. 63, No 4, pp. 797-811.
3. Hadidi, A., Azar, B. F. y Marand, H. Z. (2014). Second-Order Nonlinear Analysis of Steel
Tapered Beams. s.l. : Advances in Mechanical Engineering Volume 2014, Article ID 237983.
4. Karabalis, D. L. y Beskos, D. E. (1983). Static, dynamic and stability analysis of structures
composed of tapered members. s.l. : Computers & Structures Vol. 16, No. 6, pp.731-748.
5. Eisenberger, M. (1991). Stiffness matrices for non-prismatic members including transverse
shear. s.l. : Computers & Strictures Vol. 40, No. 4, pp. 831-835.
6. Al-Gahtani, H. J. (1996). Exact stiffnesses for tapered members. s.l. : Journal of Structural
Engineering, vol. 122, no. 10, pp. 12341239,.
7. SAP2000, Ver.15. (2011). Analysis Reference Manual. s.l. : Computers and Structures, Inc.,
Berkeley, Calif, USA.
8. Gere, J. M. y Goodno, B. J. (2009). Mechanics of Materials. s.l. : Cengage Learning, Stamford,
CT, USA.
9. Aristizabal-Ochoa, J. D. (1993). Statics, stability and vibration of non-prismatic beams and
columns. s.l. : Journal of Sound and Vibration 162(3), 441-455.
10. Al-Sadder, S. Z. y Qasrawi, H. Y. (2004). Exact secant stiffness matrix for nonprismatic
beam-columns with elastic semirigid joint connections. s.l. : Emirates Journal for Engineering
Research, 9 (2), 127-135.
11. Celigeta, J. T. (1998). Curso de Anlisis Estructural. s.l. : S.A. Eunsa. Ediciones Universidad
de Navarra.
12. Logan, D. L. (2007). A First Course in the Finite Element Method, Fourth Edition. s.l. : Nelson,
Toronto, Ontario, Canada.
13. Ghali, A., Neville, A. M. y Brown, T. G. (2009). Structural Analysis. A unified Classical and
Matrix Approach. s.l. : Spon Press, Taylor & Francis, 270 Madison Avenue, New York, NY 10016,
USA.
14. Fertis, D. G. (2006). Nonlinear Structural Engineering. s.l. : Springer, Heidelberg, Berlin,
Germany.
15. Reddy, J. N. (2006). An Introduction to the Finite Element Method. s.l. : McGraw-Hill. 1221
Avenue of the Americas, New York, NY 10020.

41

42

8.

Anexo

8.1 Integracin de la ecuacin diferencial de la viga

Euler-Bernoulli
La Figura 25 representa la configuracin deformada de una viga cantilver no
prismtica con una carga concentrada en el extremo libre. En esta figura  es la
deflexin vertical del elemento viga en cualquier punto
y es su rotacin en
cualquier punto
. Se consideran adems las siguientes relaciones
 =


,


 =


,


 = tan ,

= ,tan()

(a)

(b)

Figura 25. (a) Viga en voladizo no prismtica con una carga concentrada en el extremo libre.

En coordenadas rectangulares
y  la teora de Euler-Bernoulli para flexin se

puede escribir como[14]:
(b) Elemento infinitesimal de la viga.

T1 + ( ) V

#
# #

(44)

donde # es el momento flector producido por la carga en la viga, # es el

mdulo de elasticidad del material y # es el momento de inercia de la seccin
transversal.

Debido a que la carga de la viga puede ser arbitraria y # y # pueden ser

variables, se reescribe la ecuacin (44) de la siguiente manera


T1 + ( ) V

#
  ](
)b(
)

(45)

siendo ](
) la funcin del momento de inercia que representa la variacin de #

con  como valor de referencia y b(
) es la funcin del mdulo de elasticidad que
43

representa la variacin de # con  como valor de referencia. Si  y  son

constantes, entonces ](
) = b(
) = 1.00.

Para integrar la ecuacin (45) se hacen los siguientes cambios de variables:

 = y  = . De sta manera de la ecuacin (44) se obtiene

T1 +  V

siendo

= (
)

(
) =

(46)

#
# #

Se reescribe la ecuacin (46) de la siguiente forma

/

T1 +  V

= (
)

(47)

Multiplicando ambos lados de (47) por 
e integrando una vez, se obtiene

X

T1 +  V

= X (
) 

(48)

Se puede integrar la ecuacin (68) haciendo las siguientes sustituciones

= tan

(49)

 = sec  

(50)

Usando el elemento viga mostrado en la Figura 25(b) y aplicando el teorema de

Pitgoras, se obtiene

= T
 +   V


 = 
 +   ,

 
= x1 + . / y = T1 + (tan ) V = T1 +  V




en consecuencia

()
=



= 1/T1 +  V


(51)

y de la ecuacin (69), se tiene

sin = cos = /T1 +  V

(52)

Susituyendo las ecuaciones (49) y (50) en la ecuacin (48) y haciendo uso de

las ecuaciones (60) y (61), se obtiene

44

sec  

1 +

= X (
)

(53)

utilizando identidades triginomtricas, la ecuacin (53) se reduce a la

siguiente expresin
X cos  = X (
)

(54)

La integracin de la ecuacin (54) da como resultado

sin = (
) + W

donde la funcin (
) representa la integracin de (
).

(55)

La ecuacin (55) se puede reescribir en trminos de y  usando la

ecuacin (52). Lo que da como resultado

T1 +  V


= (
) + W

T1 + () V

= (
) + W

(56)

siendo W una constante de integracin la cual puede ser determinada de las

condiciones de contorno. Resolviendo la ecuacin (56) para  (
) se obtiene la
siguiente ecuacin
 (
) =

(
) + W

1 T(
) + WV

(57)

La integracin de la ecuacin (57) conlleva a la deflexin (
) del elemento. En

consecuencia
#

(
) = X
Z

() + W

1 T() + WV

(58)

Esto muestra que cuando # /# # es conocido y es integrable la ecuacin de

Euler-Bernoulli se puede resolver directamente para  (
) . Para resolver la
ecuacin (58) se puede usar una integracin numrica como la regla de Simpson.

8.2

Solucin de la formulacin dbil

Se considera el problema de encontrar la funcin (
) que satisface la ecuacin

diferencial de la flexin elstica de vigas (Ecuacin (59)) (bajo la hiptesis de la viga
Euler-Bernoulli para la cual las secciones planas y perpendiculares al eje de la viga
antes de la deformacin permanecen planas y ortogonales al eje despus de la
deformacin), donde  denota la deflexin transversal de la viga, = es la longitud,
(
) es la rigidez a flexin y (
) es la carga transversal distribuida.
45


 
(
)
! + (
) = 0,





0, =

(59)

Debido a que la ecuacin tiene una derivada de cuarto orden, se debe integrar
dos veces por partes para distribuir el grado de derivacin por igual entre la
variable dependiente  y la funcin de peso u*[15]. En este caso, % debe ser dos
veces diferenciable y satisfacer la forma homognea de las condiciones de contorno
esenciales.
Multiplicando (59) por la funcin de peso % e integrando en el dominio queda:
d

 (
)
X % x  (
)
! + (
)y 
= 0




Z

Integrando por partes


 (
)
 = x  (
)
!y 
,





Se llega a

% = % ,

% =

%


(60)


 (
)
(
)
!





% 
 (
)

(
)
X x.
/
(
)
(
)
! + (
)% y 
+ x%
!y = 0









 Z
Z
d

Integrando otra vez por partes


 (
)
 = x
(
)
!y 
,





Se llega a

%=

%
,


% =

= (
)

 %



(61)

 (
)



 %
 (
)

 (
)
%
 (
)

X x  (
)
+
(
)%
y


+
x%
(
)

(
)
y = 0 (62)
!











 Z
d

De la ecuacin (62) se deduce que la especificacin de  y /
constituyen las

ecuaciones de contorno esenciales y la especificacin de

 
(
)  ! (Fuerza cortante) y






 
!



(Momento lector)

(63)

constituyen las condiciones de contorno naturales para la teora de la viga

Euler-Bernoulli.
Si se considera una viga una viga empotrada en el extremo izquierdo y sujeta a
un momento flector en
= =:
46

(0) = 0,

. /
= 0,

#Z

(
)   !
= Z ,

 #d



 
x
(
)  !y
= 0 (64)



#d

donde Z es el momento flector. Dado que  y /
estn especificadas en

= 0, se requiere que la funcin de peso % y su derivada % /
cumplan la forma
homognea de las condiciones de contorno esenciales
%
/
=0

#Z

(65)

 %  (
)
%

+
(
)%



=0
!

Z

 


#d

(66)

% (0) = .

Las dos condiciones de contorno restantes en (64) son condiciones de contorno

naturales, las cuales no ponen ninguna restriccin en % y sus derivadas. En
consecuencia, la ecuacin (62) queda
d

X (
)
Z

Por ltimo, las condiciones esenciales de contorno incluyen la variable

dependiente y su derivada. En cualquier punto del contorno slo se puede
especificar una de las dos condiciones (esenciales o bien naturales), es decir estn
contrapuestas y son excluyentes entre s la flecha y el cortante en el mismo punto.
Anlogamente sucede con el giro y el momento flector.

8.3 Mtodo de Galerkin para deducir las ecuaciones de

la viga
Viga Prismtica
Se parte de la ecuacin diferencial de la flexin:


R 
+ (
) = 0,

R

0, =

(67)

(
) es la sobrecarga repartida (positiva hacia arriba) y (
) es la variable

dependiente del problema.
El mtodo de los residuos ponderados aplicado directamente a la ecuacin
diferencial puede usarse para desarrollar las ecuaciones del elemento finito. En el
mtodo de los residuos ponderados se escoge una funcin de prueba para aproximar
la variable independiente del problema definido por la ecuacin diferencial[12]. Se
escoge una funcin de prueba que se aproxime a la variable de campo (
) como

 (
) = E ,; ; (
) + Z (
)

;

En la cual ; (
) y Z (
) son funciones preestablecidas tales que la funcin de

prueba  (
) satisfaga las condiciones de contorno del problema. En cambio, los
coeficientes ,; han de ser tales que la funcin de prueba satisfaga la ecuacin
diferencial en algn sentido (bien exactamente o bien en sentido integral
ponderada).
47

En general (
)no satisfar la ecuacin diferencial de dominio que rige el

problema. Su sustitucin en ella conlleva un residuo :
 = mnimo

El mtodo de los residuos ponderados de Galerkin debe cumplir que la integral,

extendida al dominio, de los residuos multiplicados por unas funciones de
pesodebe ser nula. En otras palabras, el mtodo establece que el error o residuo
es ortogonal a las funciones de peso:
 = 0

(68)

Se define como residuo la propia expresin diferencial (67)

= 

R (
)
+ (
)

R

(69)

Utilizando el mtodo de Galerkin, se elige como funcin de interpolacin para la

variable de campo (
) en la ecuacin diferencial, la expresada en trminos de los
splines cbicos de Hermite:


(
) = TAV\{^ = TAV f g



En general, esta sustitucin conlleva un residuo no nulo. Aplicando el criterio

de Galerkin, se escogen las funciones de forma A; para que hagan las veces de las
funciones de peso . Entonces, para cada funcin cbica de Hermite se tiene que:
A;  = 0 = 

R (
)
+ (
) A; (
) = 0 (O = 1,2,3,4)

R

(70)

Integrando una vez por partes el primer trmino de se tiene

 =

% =  A; (
),

R (
)

,

R

% = 

J (
)

J

A; (
)

=  A; (
)



d

R (
)
J (
)
J (
)
(
)

(
)
X 
A
=

A


X
A; (
) 
(O = 1,2,3,4)
!
;
;

R

J Z

J

Integrando otra vez por partes la expresin anterior se obtiene

 =

% =  A; (
),

J (
)

,

J
% = 

 (
)



A; (
)

=  A; (
)



48

 R (
)
X 
A; (
)
=

R
Z

 J (
)
A; (
)
!

J Z
d

A; (
)

  (
)
  (
)
!
+

X
A; (
)


 Z


d

(71)

que expresadas con notacin abreviada son las cuatro ecuaciones siguientes
d

X  (
)A; (
)

Z

= X  (
)A; (
)
+ TA; (
) (
) A; (
) (
)VdZ

Z

(72)

una vez ms, la integracin por partes introduce las condiciones de contorno
 (
)
 =

,



% =  A; (
),

% = 

(
)


A; (
)

=  A; (
) 



 R (
)
(
)
d
X 
A; (
)
=  CA; (
) (
)DZ (A; (
) (
))dZ + (A; (
) (
))dZ X
A; (
)

R




Z

dado que

 (
) = TA V\{^ = TV\{^

(
) = TAV\{^
(
)


= TA V\{^

sustituyendo las funciones de forma de Hermite y sus derivadas en (72) y en

(70) se obtiene
d

X A; (
)TV
\{^
Z

+ X A; (
)
A; (
)5 (
) + A; (
)(
) = 0 (O = 1,2,3,4)

d

(73)

La ecuacin(73) est formada por cuatro ecuaciones (una para cada A; =

A , A , AJ y AR ). Usando la relacin TA; V = TV, estas ecuaciones se pueden expresar
en forma matricial como
d

XTV TV 
\{^ = XTAV_  
+ TA V_ 5 (
) + TAV_ (
)

Z

(74)

El primer trmino del segundo miembro representa las fuerzas nodales

equivalentes a la sobrecarga repartida. El trmino entre corchetes del segundo
miembro engloba los esfuerzos en los extremos del elemento
\A ^ 5 (
) + \A^ (
) =
d

(0)
5 (0)
(=)
5 (=)

El primer miembro de (74) es la matriz de rigidez.

d

TU V = X\^  TV 

Z

49

8.4

Cdigos de las formulaciones

Karabalis, D. L., Beskos(1983)

function [s]=Karabalis_Beskos_1983(E,L,h1,h2,tf,tw,bf)
%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin endoble-T, doblemente simtrica)
%
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis: [s]=Karabalis_Beskos_1983(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de Variables:
% x representa la posicin de la seccin a lo largo de la barra
%
% Geometra:
% Ala superior: bftf
% Alma:
twhw( =h_seccion_total - 2*tf)
% Ala inferior: bftf
%
% altura del ALMA en{el extremo inicial (x=0): h1
%
{el extremo final
(x=L): h2
%
% altura TOTAL DE LA SECCIN en{el extremo inicial (x=0): D1 Di
%{el extremo final
(x=L): D2 Dj
% E
Mdulo de elasticidad
% L
Longitud del elemento
% beta
ngulo de giro para ir de los ejes GLOBALES a los ejes
%
LOCALES;positivo si el eje local est en sentido
%
antihorario desdeel eje global
%
% -Variable de salida
%
s - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%% Propiedades Geomtricas de la Seccin
clearall; closeall
symsxz positive
tw=0.02;
bf=0.30;
tf=0.02;
h1=0.40;
h2=0.60;
L=5;
E=3e6;
beta=0;
Ax = 2*bf*tf + tw*(h1 - (x*(h1 - h2))/L);

50

Ix=simplify(tw*L^3*h1^3 + 6*bf*L^3*h1^2*tf + 12*bf*L^3*h1*tf^2 +

8*bf*L^3*tf^3
3*tw*L^2*h1^3*x
+
3*tw*L^2*h1^2*h2*x
12*bf*L^2*h1^2*tf*x + 12*bf*L^2*h1*h2*tf*x - 12*bf*L^2*h1*tf^2*x +
12*bf*L^2*h2*tf^2*x
+
3*tw*L*h1^3*x^2
6*tw*L*h1^2*h2*x^2
+
6*bf*L*h1^2*tf*x^2
+
3*tw*L*h1*h2^2*x^212*bf*L*h1*h2*tf*x^2+6*bf*L*h2^2*tf*x^2-tw*h1^3*x^3+3*tw*h1^2*h2*x^33*tw*h1*h2^2*x^3+tw*h2^3*x^3)/(12*L^3);
I1=subs(Ix,0); % valor de IG en el extremo 1
I2=subs(Ix,L); % valor de IG en el extremo 2
%% funciones de forma a flexin
g=eval(int((1/Ix),x,0,L));
h=eval(int((x/Ix),x,0,L));
GL1=eval(int((1/Ix),x,0,z));
GL=eval(int(GL1,z,0,L));
HL1=eval(int((x/Ix),x,0,z));
HL=eval(int(HL1,z,0,L));
D=GL*h-HL*g;
%% Matriz de Rigidez a Flexin
s11=E*g/D;
s12=-E*h/D;
s41=-E*GL/D;
s42=E*HL/D;
s22=E*(L*h-HL)/D;
s44=E*(L*GL-HL)/D;
s_flex=[s11 s12
s12 s22 -s12
-s11 -s12 s11
s41 s42 -s41

-s11 s41
s42
-s41
s44];

%% Funciones de forma axiales

jx=eval(int(1/Ax,x,0,x));
jL=eval(int(1/Ax,x,0,L));
f5=1-(jx/jL);
f6=(jx/jL);
f_axial=[f5,f6];
%% Matriz de Rigidez Axial
f=diff(f_axial,x);
s_axial=eval(int(f'*f*E*Ax,x,0,L));
%% Matriz de Rigidez a Flexin y Axial s* en ejes Locales
s_local=[s_axial(1,1),0,0,s_axial(1,2),0,0;
0,s_flex(1,1),s_flex(1,2),0,s_flex(1,3),s_flex(1,4);
0,s_flex(2,1),s_flex(2,2),0,s_flex(2,3),s_flex(2,4);
s_axial(2,1),0,0,s_axial(2,2),0,0;
0,s_flex(3,1),s_flex(3,2),0,s_flex(3,3),s_flex(3,4);
0,s_flex(4,1),s_flex(4,2),0,s_flex(4,3),s_flex(4,4)];
%% Matriz de rotacin para pasar de ejes Locales a ejes Globales
r= [cos(beta),sin(beta),0,0,0,0
-sin(beta),cos(beta),0,0,0,0
0,0,1,0,0,0
0,0,0,cos(beta),sin(beta),0
0,0,0,-sin(beta),cos(beta),0
0,0,0,0,0,1];

51

r(abs(r)<100*eps)=0;
%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES
s=r'*s_local*r;
end

Eisenberger, M. (1991)
function [Sm]=Eisenberger_1991(E,L,hj,hk,b)
%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin endoble-t, doblemente simtrica),
%
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis:
%
[Sm]=Eisenberger_1991(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de Variables:
% x representa la posicin de la seccin a lo largo de la barra
%
% Geometra:
% Ancho:
b
%
% % altura TOTAL DE LA SECCIN en {el extremo inicial (x=0): hj
%
{el extremo final
(x=L): hk
% E
Mdulo de elasticidad
% L
Longitud del elemento
% beta
ngulo de giro para ir de los ejes GLOBALES a los ejes
%
LOCALES; positivo si el eje local est en sentido antihorario
resp.eje global
%
% -Variable de salida
%
Sm - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%% Propiedades Geomtricas de la Seccin
clearall; closeall
tic;
symsxpositive;
b=0.032;
L=5;
E=3e6;
G=8.4e5
f=1.2
beta=0;

%Ancho de la seccin rectangular

hj=0.617;
hk=0.828;

%Altura de la seccin en el extremo inicial

%Altura de la seccin en el extremo final

%Factor de Forma por cortante

52

%% Trminos de la Matriz de Flexibilidad

g=f*E*(hj-hk)^2/(12*G*L^2);
F44=L/(E*b)*log(hj/hk)/(hj-hk);
F55=12*L^3/(E*b*(hj-hk)^3)*((1+g)*log(hj/hk)+2*hk/hj-hk^2/(2*hj^2)1.5);
F56=6*L^2/(E*b*hj^2*hk);
F66=6*L*(hj+hk)/(E*b*hj^2*hk^2);
H=F55*F66-F56^2;
A=1/F44;
B=F66/H;
C=-F56/H;
D=F55/H;
%% Matriz de Rigidez Sm en ejes Locales
Sm_local=[A 0 0 -A 0 0
0 B C+B*L 0 -B -C
0 C+B*L D+2*C*L+B*L^2 0 -C-B*L -D-C*L
-A 0 0 A 0 0
0 -B -C-B*L 0 B C
0 -C -D-C*L 0 C D];
%% Matriz de rotacin
r= [cos(beta),sin(beta),0,0,0,0;
-sin(beta),cos(beta),0,0,0,0;

0,0,1,0,0,0 ;
0,0,0,cos(beta),sin(beta),0;
0,0,0,-sin(beta),cos(beta),0;
0,0,0,0,0,1];
r(abs(r)<100*eps)=0;
%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES
Sm=r'*Sm_local*r
end

Aristizabal-Ochoa, J. D. (1993)
function [K]=Arisizabal_1993(E,L,h1,h2,tf,tw,bf,beta)
%%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin en doble-t, doblemente simtrica),
%
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis:
%
[K]=Aristizabal_1993(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de Variables:
% x representa la posicin de la seccin a lo largo de la barra
%

53

% Geometra:
% Ala superior: bftf
% Alma:
twhw( =h_seccion_total - 2*tf)
% Ala inferior: bftf
%
% altura del ALMA en {el extremo inicial (x=0): h1
%
{el extremo final
(x=L): h2
%
% altura TOTAL DE LA SECCIN {el extremo inicial (x=0): D1 Di
%
{el extremo final
(x=L): D2 Dj
% E
Mdulo de elasticidad
% L
Longitud del elemento
% beta:
ngulo de giro para ir de los ejes GLOBALES a los ejes
LOCALES;
%
positivo si el eje local est en sentido antihorario desde eleje
global
%
% -Variable de salida
%
s - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%
% Propiedades Geomtricas de la Seccin
clearall; close all
tic
symsxpositive
tw=0.02;
bf=0.30;
tf=0.02;
h1=0.40;
h2=0.60;
L=5;
E=3e6;
beta=0;
A = 2*bf*tf + tw*(h1 - (x*(h1 - h2))/L);
Ixx=simplify(tw*L^3*h1^3 + 6*bf*L^3*h1^2*tf + 12*bf*L^3*h1*tf^2 +
8*bf*L^3*tf^3
3*tw*L^2*h1^3*x
+
3*tw*L^2*h1^2*h2*x
12*bf*L^2*h1^2*tf*x + 12*bf*L^2*h1*h2*tf*x - 12*bf*L^2*h1*tf^2*x +
12*bf*L^2*h2*tf^2*x
+
3*tw*L*h1^3*x^2
6*tw*L*h1^2*h2*x^2
+
6*bf*L*h1^2*tf*x^2
+
3*tw*L*h1*h2^2*x^212*bf*L*h1*h2*tf*x^2+6*bf*L*h2^2*tf*x^2-tw*h1^3*x^3+3*tw*h1^2*h2*x^33*tw*h1*h2^2*x^3+tw*h2^3*x^3)/(12*L^3);
I1=subs(Ixx,0); % valor de IG en el extremo 1
I2=subs(Ixx,L); % valor de IG en el extremo 2
%% coeficientes bsicos de Flexibilidad
f11=(I1/L^3)*int((L-x)^2/Ixx,x,0,L);
f22=(I1/L^3)*int(x^2/Ixx,x,0,L);
f12=-I1/L^3*int(x*(L-x)/Ixx,x,0,L);
f55=int(1/(E*A),x,0,L);
%% Coeficientes de Rigidez a flexin:
k11=eval((E*I1/L)*f22/(f11*f22-f12^2));
k22=eval((E*I1/L)*f11/(f11*f22-f12^2));
k12=eval(-(E*I1/L)*f12/(f11*f22-f12^2));
k55=eval(1/f55);
k56=-k55;
% coeficientes auxiliares

54

C1=(k11+k12)/L;
C2=(k22+k12)/L;
C3=(k11+k22+2*k12)/L^2;
%% Matriz de rigidez en ejes locales:
K_local=[k55
0
0 C3 C1
0 -C3
0 C1 k11
0 -C1
k56
0
0 -C3 -C1
0 C2 k12
0 -C2

0 k56 0
0;
C2;
k12;
0 k55
0
0;
0 C3 -C2;
k22];

%% Matriz de rotacin
r=
[cos(beta),sin(beta),0,0,0,0
;
sin(beta),cos(beta),0,0,0,0;
0,0,1,0,0,0 ;
0,0,0,cos(beta),sin(beta),0;
0,0,0,-sin(beta),cos(beta),0 ;
0,0,0,0,0,1];
r(abs(r)<100*eps)=0;

%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES

K=r'*K_local*r;
end

Al-Gahtani,H. J. (1996)
function [K]=Al_Gahtani_1996(E,L,h1,h2,tf,tw,bf)
%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin en doble-t, doblemente simtrica),
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis:
%
[s]=Al-Gahtani_1996(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de Variables:
% x representa la posicin de la seccin a lo largo de la barra
%
% Geometra:
% Ancho:
b
%
% % altura TOTAL DE LA SECCIN en {el extremo inicial (x=0): hj
%
{el extremo final
(x=L): hk
% E
Mdulo de elasticidad
% L
Longitud del elemento
% beta
ngulo de giro para ir de los ejes GLOBALES a los ejes
LOCALES
%
positivo si el eje local est en sentido antihorariodesde eleje
global
%
% -Variable de salida

55

%
Sm - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%% Propiedades Geomtricas de la Seccin
clearall; close all
symsx positive
tw=0.02;
bf=0.30;
tf=0.02;
h1=0.40;
h2=0.60;
Di=0.44;
Dj=0.64;
L=5;
E=3e6;
beta=0;
% Altura de la seccin en el extremo
Di=h1+2*tf;
inicial
Dj=h2+2*tf;
% Altura de la seccin en el extremo final
Dx=Di+x*(Dj-Di)/L;
% Altura de la seccin a una distancia x
desde el extremo inicial
A=Dx*tw+(bf-tw)*2*tf;
% rea de la seccin a una distancia x
desde el extremo inicial
AY=bf*tf*(Dx-tf/2)+tw*(Dx-2*tf)*(Dx/2)+bf*tf*tf/2; % Momento esttico
desde la fibra inferior de la seccin
Yg=AY/A;
% Centro de gravedad de la seccin desde la fibra
inferior
IYx= (2*bf*tf^3+tw*(Dx-2*tf)^3)/12 + bf*tf*(Dx-tf/2-Yg)^2 + tw*(Dx2*tf)*(Dx/2-Yg)^2 + bf*tf*(tf/2-Yg)^2; %Momento de inercia respecto
del c.d.g sobre el eje Y-Y
I1=subs(IYx,0);
% momento de Inercia en el extremo 1
I2=subs(IYx,L);
% momento de Inercia en el extremo 2
A1=subs(A,0);
% Area en el extremo 1
A2=subs(A,L);
% Area en el extremo 2
%% Coeficientes de la Matriz de Rigidez
c=L/(h2-h1)*1/100;
K(1,1)=E*A1/L*c/log(1+c); %para n=1
K(1,2)=0;
K(1,3)=0;
K(1,4)=-K(1,1);
K(2,2)=E*I1/L^3*(4*(1+c)*(3+3*c+c^2));
K(2,3)=E*I1/L^2*(2*(1+c)*(3+c));
K(2,4)=0;
K(2,5)=-K(2,2);
K(2,6)=E*I1/L^2*(2*(1+c)^2*(3+2*c));
K(3,3)=E*I1/L*(4*(1+c));
K(3,4)=0;
K(3,5)=-K(2,3);
K(3,6)=E*I1/L*(2*(1+c)^2);
K(4,4)=K(1,1);
K(4,5)=0;
K(4,6)=0;
K(5,5)=K(2,2);
K(5,6)=-K(2,6);
K(6,6)=E*I1/L*(4*(1+c)^3);
fori=1:6;
j=1:6;
K(j,i)=K(i,j);

56

end
K_local=K;
%% Matriz de rotacin
r= [cos(beta),sin(beta),0,0,0,0 ;
-sin(beta),cos(beta),0,0,0,0;

0,0,1,0,0,0 ;
0,0,0,cos(beta),sin(beta),0;
0,0,0,-sin(beta),cos(beta),0 ;
0,0,0,0,0,1];
r(abs(r)<100*eps)=0;
%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES
Sm=r'*K_local*r;
end

Saka, M. P. (1997)
function [K]=Saka_1997(E,L,h1,h2,tf,tw,bf)
%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin en doble-t, doblemente simtrica),
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis:
%
[s]=Saka_1997(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de variables:
% x representa la posicin de la seccin a lo largo de la barra
%
% Geometra:
% Ancho:
b
%
% % altura TOTAL DE LA SECCIN en {el extremo inicial (x=0): hj
%
{el extremo final
(x=L): hk
% Elast Mdulo de elasticidad
% L
Longitud del elemento
% beta
ngulo de giro para ir de los ejes GLOBALES a los ejes
LOCALES;
%
positivo si el eje local est en sentido antihorario desde eleje
global
%
% -Variable de salida
%
Sm - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%% Propiedades Geomtricas de la Seccin
symsx;
tw=0.02;
bf=0.30;
tf=0.02;
h1=0.40;
h2=0.60;

57

l=5;
Elast=3e6;
beta=0;
Di=h1+2*tf;
Dj=h2+2*tf;
%% Area de la seccin transversal
psi1=tw*(Di-2*tf)+2*bf*tf;
psi2=tw*(Dj-Di)/l;
Ax=psi1+psi2*x;
%% Momento de Inercia de la seccin transversal
phi1=Di-tf;
phi2=(Dj-Di)/l;
phi3=(tw*(Di-tf)/12+bf*tf/2);
phi4=tw*phi2/12;
Ix=(phi1+phi2*x)^2*(phi3+phi4*x);
%% Trminos para obtener los elementos de la Matriz de Rigidez
g=eval(int(1/Ix,x));
g1=subs(g,0);
g2=subs(g,l);
h=eval(int(x/Ix,x));
h1=subs(h,0);
h2=subs(h,l);
G=eval(int(int(1/Ix,x),x));
G1=subs(G,0);
G2=subs(G,l);
H=eval(int(int(x/Ix,x),x));
H1=subs(H,0);
H2=subs(H,l);
E=eval(int(x^2/Ix,x));
E1=subs(E,0);
E2=subs(E,l);
g21=g2-g1;
h21=h2-h1;
G21=G2-G1;
H21=H2-H1;
E21=E2-E1;
g22=l*g2*h1-l*g1*h2;
D=G21*h21-H21*g21+g22;
C31=-h21/D;
C41=g21/D;
C32=(H21-l*h2)/D;
C42=(l*g2-G21)/D;
C33=-C31;
C43=-C41;
C34=(l*h1-H21)/D;
C44=(G21-l*g1)/D;

58

%%Matriz de Rigidez en ejes locales

K(1,1)=Elast/eval(int(1/Ax,x,0,l));
K(2,2)=Elast*(C31^2*g21+2*C31*C41*h21+C41^2*E21);
K(3,2)=Elast*(C32*C31*g21+C32*C41*h21+C42*C31*h21+C42*C41*E21);
K(3,3)=Elast*(C32^2*g21+2*C32*C42*h21+C42^2*E21);
K(5,2)=Elast*(C33*C31*g21+C33*C41*h21+C43*C31*h21+C43*C41*E21);
K(5,3)=Elast*(C33*C32*g21+C33*C42*h21+C43*C32*h21+C43*C42*E21);
K(5,5)=Elast*(C33^2*g21+2*C33*C43*h21+C43^2*E21);
K(6,2)=Elast*(C34*C31*g21+C34*C41*h21+C44*C31*h21+C44*C41*E21);
K(6,3)=Elast*(C34*C32*g21+C34*C42*h21+C44*C32*h21+C44*C42*E21);
K(6,5)=Elast*(C34*C33*g21+C34*C43*h21+C44*C33*h21+C44*C43*E21);
K(6,6)=Elast*(C34^2*g21+2*C34*C44*h21+C44^2*E21);
K(4,4)=K(1,1);
K(4,1)=-K(1,1);
fori=1:6;
j=1:6;
K(i,j)=K(j,i);
end
K_local=K
%% Matriz de rotacin
r= [cos(beta),sin(beta),0,0,0,0 ;

-sin(beta),cos(beta),0,0,0,0;

0,0,1,0,0,0 ;
0,0,0,cos(beta),sin(beta),0;
0,0,0,-sin(beta),cos(beta),0;
0,0,0,0,0,1];
r(abs(r)<100*eps)=0;
%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES
K=r'*K_local*r;
End

Mtodo de las funciones de forma (Hermite)

function [K]=Funciones_de_forma(E,L,h1,h2,tf,tw,bf)
%-------------------------------------------------------------% -Determina la matriz de rigidez LINEAL (6x6) para el elemento viga
%
2-D no prismtico (canto variable linealmente, ancho constante),
%
seccin en doble T, bisimtrica),
% -La variacin del canto es lineal, CRECIENTE.
%
% -Numeracin local de grados de libertad:
% {u_1 v_1 theta_1 u_2 v_2 theta_2}
%
^2
^5
%
|
|
%
1--->@|-------------------|@--->4
%
3
6
%
% -Sintaxis:
%
[K]=Funciones_de_Forma(E,L,h1,h2,tf,tw,bf,beta)
%
% -Descripcin de Variables:
% x representa la posicin de la seccin a lo largo de la barra
%
% Geometra:
% Ancho:
b
%
% % altura TOTAL DE LA SECCIN en {el extremo inicial (x=0): hj
%
{el extremo final
(x=L): hk
% E
Mdulo de elasticidad

59

% L
Longitud del elemento
% beta
ngulo de giro para ir de los ejes GLOBALES a los ejes
LOCALES;
%
positivo si el eje local est en sentido antihorario desde el
eje global
%
% -Variable de salida
%
Sm - matriz de rigidez elemental (6x6) en ejes GLOBALES
%--------------------------------------------------------------------%% Parmetrosgeomtricos
clearall
closeall
symsx
tw=0.02;
bf=0.30;
tf=0.02;
h1=0.40;
h2=0.60;
Di=0.44;
Dj=0.64;
L=5;
E=3e6;
beta=0;
% Canto de la seccin en el extremo
Di=h1+2*tf;
inicial
Dj=h2+2*tf;
% Canto de la seccin en el extremo final
Dx=Di+x*(Dj-Di)/L;
% Canto a una distancia x desde el
extremo inicial
Ax=Dx*tw+(bf-tw)*2*tf;
% rea a una distancia x desde el extremo
inicial
AY=bf*tf*(Dx-tf/2)+tw*(Dx-2*tf)*(Dx/2)+bf*tf*tf/2; % Momento esttico
desde la fibra inferior de la seccin
Yg=AY/Ax;
% Centro de gravedad de la seccin desde la fibra
inferior
IYx= (2*bf*tf^3+tw*(Dx-2*tf)^3)/12 + bf*tf*(Dx-tf/2-Yg)^2 + tw*(Dx2*tf)*(Dx/2-Yg)^2 + bf*tf*(tf/2-Yg)^2;
%Momento de inercia respecto
del eje Y-Y
I1=subs(IYx,0);
% Momento de Inercia en el extremo 1
I2=subs(IYx,L);
% momento de Inercia en el extremo 2
A1=subs(Ax,0);
% Area en el extremo 1
A2=subs(Ax,L);
% Area en el extremo 2
%% Funciones de Forma cbicas de Hermite
B1=-1/L;
B2=1/L;
B3=6/L^2*(2*x/L-1);
B4=2/L*(3*x/L-2);
B5=-6/L^2*(2*x/L-1);
B6=2/L*(3*x/L-1);
%% Matriz de Interpolacin de Deformaciones
B=[B1,0,0,B2/L,0,0
0,B3,B4,0,B5,B6];
D=[E*Ax,0
0,E*IYx];
%% Matriz de Rigidez en ejes locales

60

K_local=eval(int(B'*D*B,x,0,L));
%% Matriz de rotacin
r= [cos(beta),sin(beta),0,0,0,0
-sin(beta),cos(beta),0,0,0,0
0,0,1,0,0,0
0,0,0,cos(beta),sin(beta),0
0,0,0,-sin(beta),cos(beta),0
0,0,0,0,0,1];
r(abs(r)<100*eps)=0;
%% Matriz de rigidez del elemento (6x6) en EJES GLOBALES
K=r'*K_local*r;
end

61

Cdigo para nave 2_D

function f=Nave_2D
clear all;
close all;
%clc;
%% Datos de Entrada
E=3e6;
%Modulo de Elasticidad
Num_elem=4;
%Nmero de elementos de la nave
%% Coordenadas:
%Elemento1: Primer Pilar
X1=[0 0];
Y1=[0 8.5];
%Elemento 2: Primer dintel
X2=[0 5];
Y2=[8.5 10];
%Elemento 3: Segundo Dintel
X3=[5 10];
Y3=[10 8.5];
%Elemento 4: Segundo Pilar
X4=[10 10];
Y4=[8.5 0];
%% Matrices de Rigidez en coordenadas locales
X=[X1 X2 X3 X4];
Y=[Y1 Y2 Y3 Y4];
for i=1:Num_elem
L(i)=sqrt((X(2*i)-X(2*i-1))^2+(Y(2*i)-Y(2*i-1))^2);
syms x
tw=0.02;
bf=0.30;
tf=0.02;
h1=[0.4 0.6 0.4 0.6];
h2=[0.6 0.4 0.6 0.4];
A = 2*bf*tf + tw*(h1(i)- (x*(h1(i)-h2(i)))/L(i));
I=(tw*L(i)^3*h1(i)^3+6*bf*L(i)^3*h1(i)^2*tf+12*bf*L(i)^3*h1(i)*tf^2+8*
bf*L(i)^3*tf^3-3*tw*L(i)^2*h1(i)^3*x+3*tw*L(i)^2*h1(i)^2*h2(i)*x12*bf*L(i)^2*h1(i)^2*tf*x+12*bf*L(i)^2*h1(i)*h2(i)*tf*x12*bf*L(i)^2*h1(i)*tf^2*x+12*bf*L(i)^2*h2(i)*tf^2*x+3*tw*L(i)*h1(i)^3*
x^26*tw*L(i)*h1(i)^2*h2(i)*x^2+6*bf*L(i)*h1(i)^2*tf*x^2+3*tw*L(i)*h1(i)*h
2(i)^2*x^2-12*bf*L(i)*h1(i)*h2(i)*tf*x^2+6*bf*L(i)*h2(i)^2*tf*x^2tw*h1(i)^3*x^3+3*tw*h1(i)^2*h2(i)*x^33*tw*h1(i)*h2(i)^2*x^3+tw*h2(i)^3*x^3)/(12*L(i)^3);
I1=subs(I,0);
% valor de IG en el extremo 1
I2=subs(I,L(i));
% valor de IG en el extremo 2
A1=subs(A,0);
% valor de IG en el extremo 1
A2=subs(A,L(i));
% valor de IG en el extremo 2
Ax=A1*(1-x/L(i))+A2*(x/L(i));
n=3;
Ix=(I1^(1/n)*(1-x/L(i))+I2^(1/n)*(x/L(i)))^n; %Variacin parablica de
momento de inercia
% Coeficientes bsicos de Flexibilidad_Aristizabal-Ochoa

62

f11=(I1/L(i)^3)*int((L(i)-x)^2/Ix,x,0,L(i));
f22=(I1/L(i)^3)*int(x^2/Ix,x,0,L(i));
f12=(-I1/L(i)^3)*int(x*(L(i)-x)/Ix,x,0,L(i));
f55=int(1/(E*Ax),x,0,L(i));
% Coeficientes de Rigidez a flexin:
k11=eval((E*I1/L(i))*f22/(f11*f22-f12^2));
k22=eval((E*I1/L(i))*f11/(f11*f22-f12^2));
k12=eval(-(E*I1/L(i))*f12/(f11*f22-f12^2));
k55=eval(1/f55);
k56=-k55;
% coeficientes auxiliares
C1=(k11+k12)/L(i);
C2=(k22+k12)/L(i);
C3=(k11+k22+2*k12)/L(i)^2;
% Matrices de Rigidez en coordenadas locales
k_local(i).n=[k55
0
0 k56 0
0;
0 C3 C1
0 -C3
C2;
0 C1 k11
0 -C1 k12;
k56
0
0 k55
0
0;
0 -C3 -C1
0 C3 -C2;
0 C2 k12
0 -C2 k22];
end
%% Matriz de Rigidez Global
gdl_tot=(Num_elem+1)*3;
K=zeros(gdl_tot,gdl_tot);
for i=1:Num_elem;
gdl=[3*i-2 3*i-2+1 3*i-2+2 3*i-2+3 3*i-2+4 3*i-2+5];
%Grados de
libertad correspondientes a cada elemento
tetax=asind((Y(2*i)-Y(2*i-1))/L(i));
%Calcula
los ngulos teta x y teta y para cada elemento a partir de las
coordendas
tetay=acosd((Y(2*i)-Y(2*i-1))/L(i));
%Matriz de transformacin de ejes locales a ejes globales
T= [cosd(tetax) cosd(tetay) 0 0
0
0
-cosd(tetay) cosd(tetax) 0 0
0
0
0
0
1 0
0
0
0
0
0 cosd(tetax) cosd(tetay) 0
0
0
0 -cosd(tetay) cosd(tetax) 0
0
0
0 0
0
1];
K1=T'*k_local(i).n*T;
DeltaK1=zeros(gdl_tot,gdl_tot);
DeltaK1(gdl,gdl)=K1;
K=K+DeltaK1;
%Matriz de Rigidez en ejes globales
end
%% Clculo de Desplazamientos y Reacciones:
a=[1 2 3 13 14 15]';
desplazamiento=0
b=[4 5 6 7 8 9 10 11 12]';
desplazamientos a calcular
Kbb=K(b,b);
Kab=K(a,b);
w=3.5;
F=70;

%Grados de libertad con

%Grados de libertad con

%Carga de Viento
%Fuerza vertical aplicada en cubierta

63

P=[0 0 0 w*L(1)/2+w*L(2)/2 -F w*L(1)^2/12+w*L(2)^2/12 w*L(2)/2 -F w*L(2)^2/12 0 -F 0 0 0 0]'; %Cargas externas aplicadas

Pb=P(b);
Db=Kbb\Pb;
Pa=Kab*Db;
D=zeros(gdl_tot,1);
D(b)=Db;
num=floor(1:gdl_tot);
disp('
gdl
Movimientos');
Movimientos=[num' D]
P=zeros(gdl_tot,1);
P(a)=Pa;
disp('
gdl
Reacciones');
Reacciones=[num' P]
end

64