Matlab Para Dinamica de Estructuras RevB

PROHIBIDA SU REPRODUCCION TOTAL O PARCIAL

UNIVERSIDAD TÉCNICA FEDERICO SANTA MARÍA Apuntes de MatLab para Dinámica de Estructuras USM Campus Santiago II Semestre 2012

Apuntes de MatLab para Dinámica de Estructuras – II Semestre 2012 – USM Campus Santiago – RevB // Ricardo Bustamante H.


NOCIONES BÁSICAS

3

NOMENCLATURAS Y COMANDOS GRÁFICOS PROGRAMACIÓN DECLARACIÓN DE FUNCIÓN DECLARACIÓN DE FOR SIMPLE DECLARACIÓN DE FOR ANIDADA DECLARACIÓN DE WHILE DECLARACIÓN DE IF DECLARACIONES INTERACTUANDO

3 11 15 15 15 16 16 17 17

ANÁLISIS ESTRUCTURAL

18

CASO DE ANÁLISIS MATRIZ DE RIGIDEZ VIGAS COLUMNAS MATRIZ DE RIGIDEZ GLOBAL CONDENSACIÓN ESTÁTICA

18 18 19 20 20 26

DINÁMICA ESTRUCTURAL

28

SUBRUTINA ‘ODE45’ FACTORES DE AMPLIFICACIÓN DINÁMICA (FAD) RESPUESTA ANTE CARGA “RAMPA-ESCALÓN” RESPUESTAS A PARTIR DE UN REGISTRO DE ACELERACIONES ESPECTRO LINEAL DE ACELERACIÓN ESPECTRO NO LINEAL DE ACELERACIÓN SUPERPOSICIÓN MODAL

28 31 37 39 43 54 65



Nociones básicas Nomenclaturas y comandos La interfaz al ejecutar MatLab es la que se detalla a continuación.

Se aprecia que aparece el indicador de comandos >>, el cual indica que el programa está listo para recibir instrucciones. MATLAB trabaja esencialmente con matrices numéricas rectangulares. La manera más fácil de entrar matrices pequeñas es enumerando los elementos de ésta de tal manera que:   

los elementos estén separados por blancos ó comas. los elementos estén cerrados entre corchetes, [ ]. muestre el final de cada fila con ; (punto y coma).

Así, algunos caracteres especiales de MatLab son:  

[ ] Se utilizan para formar vectores y matrices ( ) Define precedencia en expresiones aritméticas. Encierra argumentos de funciones en forma usual



  

, Separador de elementos de una matriz, argumentos de funciones y declaraciones en líneas con declaraciones múltiples ; Termina filas de una matriz, separador de declaraciones % Permite agregar un comentario comentario

Entonces un ejemplo sería ingresar una matriz de dimensiones 3x3, con números correlativos, como: [ 1 2 3; 4 5 6; 7 8 9 ]. El resultado, como lo entrega el programa, se indica a continuación: ans = 1 4 7

2 5 8

3 6 9

Ahora bien, se puede guardar esta matriz bajo un nombre. Para ello se tiene que ingresar el nombre de la variable seguida por una igualdad. Así, la matriz queda como: A=[ 1 2 3; 4 5 6; 7 8 9 ]. El resultado que entrega el programa es A = 1 4 7

2 5 8

3 6 9

Usualmente no se requiere que el programa entregue los resultados. Para no presentar el valor de la variable que se creó, se debe agregar un punto y coma (;) al final del comando. Como alternativa se puede crear la matriz sin la necesidad de puntos y comas como separadores. Así, la matriz A también se puede generar como: A =

[1 2 3 4 5 6 7 8 9];

Si se requiere una fila o columna de la matriz A, o tal vez un valor particular de la matriz o añadir una fila adicional, se realiza lo siguiente: 

se requiere el segundo vector columna: A(:,2) ans = 2 5 8



se requiere el primer vector fila: A(1,:) ans = 1 2 3



se requiere el valor del elemento en la fila 2 y columna 3: A(2,3) ans = 8



se requiere añadir una fila con los valores b=[ 10 11 12], redefiniendo A: A=[A;b]: A = 1 4 7 10

2 5 8 11

3 6 9 12



También se pueden generar vectores con instrucciones. Para ello, una alternativa es definir x=[1:5]. Esta instrucción permite generar un vector fila de 5 columnas. También se puede omitir los paréntesis cuadrados para simplificar el proceso. Así, se obtiene: x = 1 2 3 4 5

El vector transpuesto se genera agregando el carácter apóstrofe (‘) a la variable, es decir, x’ ans = 1 2 3 4 5

Se pueden crear vectores generando una secuencia, con un punto de inicio y uno de término. Un ejemplo es definir un vector cuyo valor inicial sea 5 y que disminuya su valor en una unidad hasta llegar a 1. De esta manera, x=5:-1:1 x = 5 4 3 2 1

Con anterioridad se evidenció la posibilidad de pedir valores específicos a la matriz. De igual manera, se pueden reasignar valores dentro de la matriz sin necesidad de cambiarla. Para ello se tienen las siguientes opciones:  cambiar un valor en específico de manera directa: A(3,3)=10 A = [1 2 3 4 5 6 7 8 10];



cambiar un valor en específico de manera directa: A(3, 3) = A(1, 3) + A(3, 1) A = [1 2 3 4 5 6 7 8 10];

Ahora bien, un índice también puede ser un vector. Si x y v son vectores, entonces x(v) es [x(v(1)),x(v(2)), ...,x(v(n))]. Para matrices, los índices de vectores permiten acceso a submatrices contiguas y no contiguas. Por ejemplo, suponga que B es una matriz 12 por 12. Entonces B(1:5, 2) especifica la submatriz 5 x 1, ó vector columna, que consiste de los primeros cinco elementos en la segunda columna de B. De igual manera, B(1:5, 9:12) es la submatriz 5 x 4 de las primeras cinco filas y las últimas cuatro columnas. Las operaciones suma (+) y resta (-) son definidas para las matrices siempre y cuando éstas tengan la misma dimensión. El resultado de la operación es por defecto almacenado en la variable ans, la cual se puede almacenar en otra variable. Así, definiendo A=[1 2 3;4 5 6]; y B=[6 5 4; 3 2 1]; se puede originar C=A+B. C = 7 7

7 7

7 7



La operación de multiplicación de matrices (*) está definida siempre que el número de columnas de la primera matriz sea igual al número de filas de la segunda matriz. De esta manera, D=A*B sólo es posible si D=A’*B. D = 18 27 36

13 20 27

8 13 18

La operación de multiplicación de matriz con vector es un caso especial de la multiplicación de matrices. En la operación de división de matrices, si A es una matriz cuadrada no-singular, entonces A\B y B/A corresponden a la multiplicación izquierda y derecha de B por el inverso de A, esto es, inv(A) * B y B * inv(A) respectivamente. Así, 

X = A\B es una solución a A * X = B



X = B/A es una solución a X * A = B

>> A = [3,4,7;5,2,-9;-1,13,3] A = 3 5 -1

4 2 13

7 -9 3

>> B = [6;1;8] B = 6 1 8 >> x=A\B x = 0.5123 0.5848 0.3034 >> x=inv(A)*B x = 0.5123 0.5848 0.3034

Notar que A\B es definido cuando B tiene la misma cantidad de filas que A. Si A es cuadrada, el método usado para su resolución es la Eliminación Gaussiana. El resultado es una matriz con las mismas dimensiones que B. Ahora bien, si A no es cuadrada, el procedimiento se factoriza utilizando la ortogonalización de Householder con pivoteo de columnas. Para dividir los valores de una matriz por un escalar se requiere utilizar un operador que le indique al arreglo que el denominador debe dividir cada valor. Así, si A=[1 2;3 4] y se requiere dividir cada valor por 2, la operación es: A=[1 2;3 4]./2



A = 0.5000 1.5000

1.0000 2.0000

Para generar matrices o arreglos con valores aleatorios existe un comando llamado rand(). Los valores están entre cero y uno. De esta manera se tiene: 

matriz de dimensiones 3x3 con rand(3) A = 0.8235 0.6948 0.3171



0.9502 0.0344 0.4387

0.3816 0.7655 0.7952

vector de dimensiones 3x1 con rand(3,1) A = 0.1869 0.4898 0.4456



matriz de nxn de ceros: zeros(n)

>> A=zeros(4) A = 0 0 0 0



0 0 0 0

0 0 0 0

0 0 0 0

0 0 0

0 0 0

0 0 0

matriz de mxn de ceros: zeros(m,n)

>> B=zeros(3,4) B = 0 0 0



regresa el número de filas y columnas de la matriz A ó B: size(A) ó size(B)

>> size(A) ans = 4

4

3

4

>> size(B) ans = >> [m,n]=size(B) m = 3 n = 4



matriz de nxn de unos: ones (n)

>> C=ones(4) C = 1 1

1 1

1 1

1 1



1 1



1 1

1 1

1 1

matriz de mxn de unos: ones(m,n)

>> D=ones(3,4) D =



1 1 1

1 1 1

1 1 1

1 1 1

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

matriz identidad de nxn: eye(n)

>> E=eye(4) E =



matriz identidad de mxn: eye(m,n)

>> F=eye(3,4) F = 1 0 0



0 1 0

0 0 1

0 0 0

matriz transpuesta: transpose(A)

>> A=rand(3) A = 0.8147 0.9058 0.1270

0.9134 0.6324 0.0975

0.2785 0.5469 0.9575

0.8147 0.9134 0.2785

0.9058 0.6324 0.5469

0.1270 0.0975 0.9575

>> B=transpose(A) B =



determinante de una matriz: det(A)

>> det(A) ans = -0.2767



inversa de una matriz: inv(A)

>> inv(A) ans =



-1.9958 2.8839 -0.0291



3.0630 -2.6919 -0.1320

-1.1690 0.6987 1.1282

valores propios de una matriz: eig(A)

>> eig(A) ans = -0.1879 1.7527 0.8399



vectores y valores propios de una matriz: [V,D]=eig(A)

>> [V,D]=eig(A) V = 0.6752 -0.7375 -0.0120

-0.7134 -0.6727 -0.1964

-0.5420 -0.2587 0.7996

D = -0.1879 0 0



0 1.7527 0

0 0 0.8399

obtiene la diagonal de una matriz: diag(D)

>> diag(D) ans = -0.1879 1.7527 0.8399



agrupar de forma creciente: sort(v)

>> x=rand(5,1) x = 0.6948 0.3171 0.9502 0.0344 0.4387 >> v=sort(x) v = 0.0344 0.3171 0.4387 0.6948 0.9502





obtener la suma de los términos de un arreglo: sum(x)

w=sum(x) w = 2.4353



obtener el máximo valor de un arreglo: max(x)

y=max(x) y = 0.9502



obtener el mínimo valor de un arreglo: min(x)

>> y=min(x) y = 0.0344



redondea los términos: round(x)

>> r=round(x) %redondea por defecto al entero más cercano ó 0 r = 1 0 1 0 0 >> r=round(10*x) %idem anterior pero ahora se multiplica el arreglo por 10 r = 7 3 10 0 4



Los operadores relacionales son < <= > >= == =~



menor que menor o igual a mayor que mayor o igual a igual a no igual a

Comando de limpieza de variables - clear - clear all - clc



Gráficos La estructura básica del gráfico se compone por un comando que crea una gráfica de vectores o columnas de matrices a través de ‘plot’. Una vez ingresado el comando y los vectores a representar gráficamente es posible agregar descripciones al gráfico, a través de títulos, encabezamientos de ejes, líneas segmentadas o texto mediante etiquetas. De igual manera se puede asignar un color específico al trazado y/o con algún estilo. Cabe destacar que por defecto el gráfico tiene una línea continua y de color azul. Para ello se enumeran las opciones disponibles: Colores y

amarillo

m

magenta

c

cyan

r

rojo

g

verde

b

azul

w

blanco

k

negro

Tipología línea .

punto

o

circulo

x

cruz

+

signo '+'

*

estrella

s

cuadrado

d

diamante

v

triangulo

-

línea solida

:

segmentado

-.

línea - punto

--

línea - línea

Se considerará el siguiente ejemplo para ilustrar el uso de la función plot y sus etiquetas. Sea x=1:8; y=2.^x; plot(x)

plot(y)

Si x es un vector, plot(x) produce una gráfica lineal de los elementos de x versus el índice de estos. Si especifica dos vectores como argumentos, plot(x, y) produce una gráfica de y versus x como se indica en la figura.



Para introducir las etiquetas es necesario colocar (,) después del vector y abrir una llamado a la etiqueta, ingresando el llamado entre (‘’) plot(y,’m’) plot(y,’mx’)

El gráfico de la izquierda grafica el vector ‘y’, agregándole color magenta a la línea trazada. El gráfico de la derecha mantiene el color magenta pero cambia la línea por cruces; es decir, añade una característica adicional insertando una nueva etiqueta de manera continua. plot(y,’d’)

plot(y,’k--’)

El gráfico de la izquierda grafica el vector ‘y’, indicando sólo los valores del vector con un símbolo de diamante. El gráfico de la derecha segmenta el gráfico y le cambia el color a negro.



Para agregar títulos al gráfico plot(y,'k.-'); title('Gráfico y(x)');

plot(y,'k.-'); title('Gráfico y(x)'); xlabel('Eje X'); ylabel('Eje Y');

De manera adicional se puede agregar cambios en los espesores de línea. Algunas de las etiquetas son: ,'LineWidth',2,

espesor de linea en 2

,'MarkerEdgeColor','r',

dato en color rojo

,'MarkerSize',30,

tamaño dato

Así, se puede generar gráficos como: plot(y,'k.','LineWidth',2,'MarkerEdgeColor','r','MarkerSize',15); title('Gráfico y(x)');



Para graficar de manera simultánea dos gráficos se puede utilizar un comando que permita agregar tantas funciones cómo se requiera. Para ello se definirá lo siguiente: >> >> >> >> >> >> >> >>

x=0:6; y1=x.^2; y2=sqrt(x); y3=x; plot(y1); hold on plot(y2); plot(y3); hold off



Programación Declaración de función La idea de es utilizar funciones predeterminadas y de uso recurrente agregando sólo una línea. Para los casos que se verán en Dinámica de Estructuras se implementarán para:   

Llamar una matriz de rigidez elemental de un elemento viga o columna Llamar un registro Llamar una función específica previamente creada y guardada, por ejemplo, una función que resuelva la ecuación diferencial de movimiento mediante la integración directa con aceleración constante

Por ejemplo, si se crea un archivo .m a través de una función (File->New->Function) se puede agregar una matriz de rigidez de un elemento viga en posición horizontal con 6 grados de libertad, como se detalla a continuación: function Kviga = Kviga(EA,EI,L) Kviga = [EA/L,0,0,-EA/L,0,0; 0,12*EI/L^3,6*EI/L^2,0,-12*EI/L^3,6*EI/L^2; 0,6*EI/L^2,4*EI/L,0,-6*EI/L^2,2*EI/L; -EA/L,0,0,EA/L,0,0; 0,-12*EI/L^3,-6*EI/L^2,0,12*EI/L^3,-6*EI/L^2; 0,6*EI/L^2,2*EI/L,0,-6*EI/L^2,4*EI/L]

De esta manera se define una matriz de rigidez elemental para una viga horizontal de propiedades ‘EA’, ‘EI’ y largo ‘L’. Para hacer uso de ella se debe llamar por el nombre ‘Kviga’, agregando las propiedades requeridas en su definición. Definiendo EA, EI & L, se obtiene: >> EA=10000; %EA de la viga >> EI=500000; %EI de la viga >> L=500; %L de la viga >> Kviga(EA,EI,L) Kviga = 1.0e+003 * 0.0200 0 0 0 0.0000 0.0120 0 0.0120 4.0000 -0.0200 0 0 0 -0.0000 -0.0120 0 0.0120 2.0000

-0.0200 0 0 0.0200 0 0

0 -0.0000 -0.0120 0 0.0000 -0.0120

0 0.0120 2.0000 0 -0.0120 4.0000

Declaración de FOR simple >> for i=1:101 x(i)=2*pi*(i-1)/100; end >> y=sin(x) >> plot(x,y)



Declaración de FOR anidada Se utiliza cuando se posee más de 1 parámetro al cual se requiere hacer un barrido. Por ejemplo, se puede crear una matriz de dimensiones ‘mxn’, creando cada término en función de las dimensiones de la matriz. Así, si m=4 y n=5, se tiene >> for i = 1:m for j = 1:n A(i, j) = 1/(i+j- 1); end end A A= 1.0000 0.5000 0.3333 0.2500

0.5000 0.3333 0.2500 0.2000

0.3333 0.2500 0.2000 0.1667

0.2500 0.2000 0.1667 0.1429

0.2000 0.1667 0.1429 0.1250

Declaración de WHILE Este comando permite a una o más instrucciones repetirse un número indefinido de veces, bajo el control de una condición lógica. La ventaja radica en que WHILE tiene una condición de parada definida, mientras que otros comandos, como IF, pueden no tener la condición de parada, lo que implica que nunca deja de buscar. Como ejemplo, suponer se desea iterar para encontrar un número tal satisfaga una suma. Así, >> e=1.0; while (1.0+e)>1.0001 e=e/2.0; end >> e e= 6.1035e-005



Declaración de IF Permite establecer casos particulares en los cuales si se cumple una afirmación se proceda a generar una rutina con el valor al caso que se cumple. if i==j rho=1; else rho=-1 end

Declaraciones interactuando Supongamos que se tienen las frecuencias de un marco de 5 pisos ordenadas como un vector columna llamado w, tal que w= 31.4159 100.5183 183.8185 276.7465 360.8302

Si se requiere realizar un análisis CQC (Complete Quadratic Combination, para combinar los resultados modales por la técnica de Combinación Cuadrática Completa descrita por Wilson, Der Kiureghian, y Bayo (1981)), se debe determinar el coeficiente de acoplamiento modal, cuya es expresión corresponde a

ij 



i

8 2 i   j  i 3 j 3

2

  j 2   4 2i j i   j  2

2

La expresión obliga a ocupar los comandos anteriormente vistos para generar los coeficientes de acoplamiento modal. for i=1:5 for j=1:5 if i==j rho=1; else rho=(8*bi^2*(w(i)+w(j))*(w(i)^3*w(j)^3)^0.5)/((w(i)^2w(j)^2)^2+4*bi^2*w(i)*w(j)*(w(i)+w(j))^2); end end



Análisis estructural Caso de análisis Suponga que deben calcular el desplazamiento de un edificio de 5 pisos. Ese edificio tiene ciertas cualidades que se detallan abajo. Lo importante es que deben conocer el desplazamiento de techo y los desplazamientos de entrepiso (drifts), por lo cual les interesa saber los desplazamientos desde 1 a 5. ¿Pero qué sucede? Hay más grados de libertad que contribuyen a la rigidez del edificio. Detallando los grados de libertad del marco de 5 pisos, donde las vigas son axialmente indeformables, los grados de libertad a considerar son:

Esa rigidez del edificio se puede expresar como una rigidez lateral (que es la que nos interesa porque es la que mueve toda la masa), una rigidez al giro y una rigidez a la deformación vertical. Pues de esta forma se entiende que el problema a resolver no son los 25 grados de libertad, sino que los 5 que relacionan el problema que debo solucionar. Así, sabiendo que tengo una matriz de grandes dimensiones, ordeno convenientemente los grados de libertad para luego aplicar condensación a los grados de libertad requeridos. ¿Qué significa eso? Que la matriz de 25x25 se reducirá a 5x5, de los cuales se hará una equivalencia en comportamiento dada su condición estática.

Matriz de Rigidez Se requiere obtener la matriz de rigidez de un edificio de 5 pisos compuesto por columnas y vigas. Las vigas se considerarán flexurablemente deformables y las columnas axial-flexurablemente deformables.



Para mayor facilidad en el desarrollo se consideraran primero los grados de libertad que representan el desplazamiento lateral de cada piso, para luego considerar el resto. De esta forma los primero cinco términos de la matriz de rigidez resultante serán los que se usarán para determinar la respuesta a nivel de piso debido al registro. La siguiente figura indica los grados de libertad globales considerados para el edificio de 5 pisos.

Vigas La matriz de rigidez elemental de la viga considera sólo los grados de libertad axiales, por lo que el elemento queda como lo indica la figura:

Considerando la matriz elemental para una viga horizontal y tomando en cuenta sólo los 4 grados de libertad locales indicados, la matriz de rigidez local del elemento viga se expresa como

6 3L   6 3L 3L 2 L2 3L L2  EI  kv  2 c 3 v   6 3L  L  6 3L   2 3L 2 L2  3L L



Columnas En relación a la matriz de rigidez de la columna, se consideran todos los grados de libertad de la matriz de rigidez elemental. De esta forma, se tiene que el elemento columna se expresa en coordenadas locales como:

Considerando los grados de libertad locales de la columna, la matriz de rigidez local del elemento columna se expresa como  Ec I c  12 h 3 i   0    Ec I c  6 2 hi kc    Ec I c  12 3 hi   0    EI  6 c 2c hi 

0

6

Ec I c hi 2

Ec Ac hi



0

0



Ec Ac hi

Ec I c hi 2

0

Ec I c hi 2

12

Ec I c hi 3

0

0

6

0 2

0

6

4

Ec Ac hi

Ec I c hi 3

Ec I c hi

0

0

12

0

Ec I c hi

6

Ec I c hi 2

Ec Ac hi 0

Ec I c  hi 2    0   EI  2 c c  hi  Ec I c  6 2  hi   0   Ec I c  4  hi  6

Se debe hacer la diferencia sobre la altura del primer piso con la del resto de los pisos, para lo cual se generarán matrices de rigidez para el primer piso y para el resto. Todas las vigas tienen las mismas propiedades, por lo cual no se hará diferencia con las matrices de rigidez de éstas. Los valores a utilizar son

 tonf Ec  2038.902  2  cm

   29000  ksi 

Ac  62 cm2  I c  9230 cm4 

I v  12250 cm4 

Lv  750 cm

  400  cm , para el primer piso hi    350  cm , para el resto delos pisos

Matriz de rigidez global Realizado lo anterior, será necesario establecer la relación entre los grados de libertad locales y globales. Para ello se enumeran los elementos, con la finalidad de establecer los vectores de



conectividad que permitan ensamblar la matriz de rigidez del sistema. La figura que se indica a continuación establece la numeración de los elementos. Se requieren 5 vectores de conectividad para las 5 vigas y 10 vectores de conectividad para las 10 columnas.

Utilizando lo visto en el capítulo 1, se definirán 2 funciones para llamar las matrices de rigidez de la viga y de la columna. Así, creando un nuevo archivo de extensión .m, se tiene Matriz de rigidez elemental viga 4x4 function kviga4x4 = kviga4x4(Ec,Iv,L) kviga_ad4x4=Ec*Iv/L^3*[12 6*L -12 6*L; 6*L 4*L^2 -6*L 2*L^2; -12 -6*L 12 -6*L; 6*L 2*L^2 -6*L 4*L^2];



Matriz de rigidez elemental columna function kcol = kcol(EI,EA,h) kcol = [12*EI/h^3,0,6*EI/h^2,-12*EI/h^3,0,6*EI/h^2; 0,EA/h,0,0,-EA/h,0; 6*EI/h^2,0,4*EI/h,-6*EI/h^2,0,2*EI/h; -12*EI/h^3,0,-6*EI/h^2,12*EI/h^3,0,-6*EI/h^2; 0,-EA/h,0,0,EA/h,0; 6*EI/h^2,0,2*EI/h,-6*EI/h^2,0,4*EI/h];

Notar que se consideran los parámetros Ec, Iv & L como los datos de entrada cuando se llame a la matriz de rigidez de la viga, mientras que para la columna los datos de entrada son EI, EA & h. La idea será generar que la matriz de rigidez global de la estructura se armará con las matrices de rigidez de cada elemento (método de rigidez directa). Esto se implementará mediante vectores que relacionarán el grado de libertad de cada elemento con una ubicación dentro de la matriz de rigidez global, de dimensiones 25x25. Si la viga 1 tiene los grados de libertad 6, 7, 8 & 9, el vector será p=[6:9]' y lo que hará será vincular, en el orden específico en el cual se ingresan los valores en el vector, la matriz de 4x4 de la viga del primer piso y ubicarla dentro de la matriz de dimensiones 25x25 entre las filas 6 y 9, y las columnas 6 y 9. Para las columnas se procede de la misma forma, salvo las columnas del primer piso, que requieren que se añadan grados de libertad adicionales. Esto sucede porque la matriz a la cual se llama tiene dimensiones de 6x6 y el vector p=[0;0;0;1;6;7] tendría sólo 3 grados de libertad, lo que genera un problema de dimensiones al generar la matriz. De esta forma, la manera más simple de generar una solución es añadiendo 5 grados de libertad en la base. De esta forma, se tienen 30 grados de libertad.



Así, Kg=zeros(30); %definición de una matriz de 25x25 para ir completando %Para el primer piso p=[6:9]'; %posicion de los GL de la viga Kg(p,p)=Kv+Kg(p,p); %aporte de la primer piso p=[26;27;28;1;6;7]; %posición de la columna izquierda Kg(p,p)=Kcp1+Kg(p,p); %aporte de la columna izquierda p=[26;29;30;1;8;9]; %posición de la columna derecha Kg(p,p)=Kcp1+Kg(p,p); %aporte de la columna derecha %Para el segundo piso p=[10:13]'; %posicion de los GL de la viga Kg(p,p)=Kv+Kg(p,p); %aporte del segundo piso p=[1;6;7;2;10;11]; %posición de la columna izquierda Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna izquierda p=[1;8;9;2;12;13]; %posición de la columna derecha Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna derecha %Para el tercer piso p=[14:17]'; %posicion de los GL de la viga Kg(p,p)=Kv+Kg(p,p); %aporte del tercer piso p=[2;10;11;3;14;15]; %posición de la columna izquierda Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna izquierda p=[2;12;13;3;16;17]; %posición de la columna derecha Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna derecha %Para el cuarto piso p=[18:21]'; %posicion de los GL de la viga Kg(p,p)=Kv+Kg(p,p); %aporte del cuarto piso p=[3;14;15;4;18;19]; %posición de la columna izquierda Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna izquierda p=[3;16;17;4;20;21]; %posición de la columna derecha Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna derecha %Para el quinto piso p=[22:25]'; %posicion de los GL de la viga Kg(p,p)=Kv+Kg(p,p); %aporte del quinto piso p=[4;18;19;5;22;23]; %posición de la columna izquierda



Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna izquierda p=[4;20;21;5;24;25]; %posición de la columna derecha Kg(p,p)=Kcp2+Kg(p,p); %aporte de la columna derecha



Matriz de rigidez global Kg queda como se indica la figura



Condensación estática Ahora bien, mediante condensación estáticamente se debe reducir la matriz a los grados de libertad asociados al desplazamiento horizontal del marco. En este caso se condensarán los asociados al desplazamiento lateral del marco.

 K aa K  ba

K ab  u   F     Kbb     0 

Con u=grados de libertad activos (desplazamientos laterales) y el resto de los grados de libertad. De esta manera,

 Kaa   u   K ab      F   Kba   u   Kbb      0 Finalmente la rigidez lateral queda expresada como

 K    Kaa    Kab    Kbb   Kba  1

Se condensará para la submatriz de dimensiones 25x25, considerando que sólo se agregaron estos términos para determinar un procedimientos de ensamble de la matriz de rigidez global eficiente. Truncado los últimos 5 GL, se tiene a=[1:5]; % GL activos b=[6:25]; % GL pasivos K=Kg(a,a)-Kg(a,b)*inv(Kg(b,b))*Kg(b,a) %condensacion a GL laterales

K= 14.7948 -9.8816 -9.8816 15.4632 2.6659 -10.0323 -0.4655 2.6500 0.0756 -0.3583

2.6659 -0.4655 -10.0323 2.6500 15.4597 -9.8125 -9.8125 14.1538 2.0987 -6.5904

0.0756 -0.3583 2.0987 -6.5904 4.7790

Una alternativa al procedimiento anterior consiste en generar los vectores de conectividad de manera individual. El problema radica en la lentitud en relación al procedimiento anterior.



% vigas Kv=kviga4x4(Es,Iv,L); % columnas EA=Es*Ac; EI=Es*Ic; Kcp1=kcol(EI,EA,hp1); Kcp2=kcol(EI,EA,hp2); % Matrices de compatibilidad para cada barra % vigas a1=zeros(4,25); a1(1,6)=1; a1(2,7)=1; a1(3,8)=1; a1(4,9)=1; a2=zeros(4,25); a2(1,10)=1; a2(2,11)=1; a2(3,12)=1; a2(4,13)=1; a3=zeros(4,25); a3(1,14)=1; a3(2,15)=1; a3(3,16)=1; a3(4,17)=1; a4=zeros(4,25); a4(1,18)=1; a4(2,19)=1; a4(3,20)=1; a4(4,21)=1; a5=zeros(4,25); a5(1,22)=1; a5(2,23)=1; a5(3,24)=1; a5(4,25)=1; % columnas a6=zeros(6,25); a6(4,1)=-1; a6(5,6)=1; a6(6,7)=1; a7=zeros(6,25); a7(4,1)=-1; a7(5,8)=1; a7(6,9)=1; a8=zeros(6,25); a8(1,1)=-1; a8(2,6)=1; a8(3,7)=1; a8(4,2)=-1; a8(5,10)=1; a8(6,11)=1; a9=zeros(6,25); a9(1,1)=-1; a9(2,8)=1; a9(3,9)=1; a9(4,2)=-1; a9(5,12)=1; a9(6,13)=1; a10=zeros(6,25); a10(1,2)=-1; a10(2,10)=1; a10(3,11)=1; a10(4,3)=-1; a10(5,14)=1; a10(6,15)=1; a11=zeros(6,25); a11(1,2)=-1; a11(2,12)=1; a11(3,13)=1; a11(4,3)=-1; a11(5,16)=1; a11(6,17)=1; a12=zeros(6,25); a12(1,3)=-1; a12(2,14)=1; a12(3,15)=1; a12(4,4)=-1; a12(5,18)=1; a12(6,19)=1; a13=zeros(6,25); a13(1,3)=-1; a13(2,16)=1; a13(3,17)=1; a13(4,4)=-1; a13(5,20)=1; a13(6,21)=1; a14=zeros(6,25); a14(1,4)=-1; a14(2,18)=1; a14(3,19)=1; a14(4,5)=-1; a14(5,22)=1; a14(6,23)=1; a15=zeros(6,25); a15(1,4)=-1; a15(2,20)=1; a15(3,21)=1; a15(4,5)=-1; a15(5,24)=1; a15(6,25)=1; % Matriz de rigidez Kvigas=a1'*Kv*a1 + a2'*Kv*a2 + a3'*Kv*a3 + a4'*Kv*a4 + a5'*Kv*a5; %aporte de las vigas Kcolumnas=a6'*Kcp1*a6 + a7'*Kcp1*a7 + a8'*Kcp2*a8 + a9'*Kcp2*a9 + a10'*Kcp2*a10 + a11'*Kcp2*a11 + a12'*Kcp2*a12 + a13'*Kcp2*a13 + a14'*Kcp2*a14 + a15'*Kcp2*a15; %aporte de las columnas Kt=Kvigas+Kcolumnas; %matriz de rigidez completa a=[1:5]; %grados de libertad horizontales b=[6:25]; %resto de los grados de libertad del sistema Kaa=Kt(a,a); Kab=Kt(a,b); Kba=Kt(b,a); Kbb=Kt(b,b); K=Kaa-Kab*inv(Kbb)*Kba; %matriz de rigidez condensada a grados de libertad laterales con masa asociada Kt=Kv+Kc; %matriz de rigidez completa



Dinámica estructural Subrutina ‘ode45’ ‘ode45’ corresponde a un método adaptativo; es una combinación de un método Runge-Kutta de ordenes cuatro y cinco de cinco y seis evaluaciones respectivamente. La combinación de estos métodos permite estimar el error en la aproximación numérica en cada paso y así la subrutina puede ajustar el largo de paso "h" en forma dinámica para mantener el error global menor de una tolerancia especificada por el usuario. Requiere del nombre de la función ‘f’, los tiempos inicial y final, la condición inicial, la tolerancia para cada iteración, y un indicador si la rutina imprime o no los resultados calculados. Normalmente este indicador se toma como cero que es su valor por omisión. A continuación se detalla un ejemplo que permite visualizar su uso. En específico, se resolverá la ecuación general de movimiento para un caso particular. Sea la ecuación del movimiento:

mx(t )  cx (t )  kx(t )  F0 sin  t , 0  t  10 Considerar el siguiente caso: m = 10; c = 0,05; k = 50; Fo = 5 Las condiciones iniciales serán las siguientes:

x(0)  0

x (0)  0 Para usar ode45 con ecuaciones de 2do orden, como es el caso de la ecuación general de movimiento:

mx  cx  kx  F (t )  x(t 0 )  x0   x (t 0 )  x1  se debe convertir la ecuación en un sistema equivalente de primer orden, haciendo y1 = x, y2 = x’:

y '1  y 2    k  y1  c  y 2  F (t ) y  2  m 



Con

y1 (t 0 )  x0 y 2 (t 0 )  x1 Luego se debe crear y guardar un archivo .m que retornará una función evaluada en cada t. Para el caso del ejemplo, se procedería de la siguiente forma: function dy = dfun(t,y) Fo=5; A=Fo*sin(pi*t); c=0.05; k=50; m=10; dy=zeros(2,1); dy(1)=y(2); dy(2)=(A-c*y(2)-k*y(1))/m; end

Notar que el sistema está descrito en las dos últimas filas antes de end y lo anterior es para definir los coeficientes del sistema. Este archivo se guardará como “dfun.m”. Finalmente, para obtener los resultados, se tiene que ingresar lo siguiente: [t,x]=ode45(@”nombre del archivo.m”,[to tf],[x(0) x’(0)]);

t0 y tf son el tiempo inicial y final del análisis y x(0) y x’(0) las condiciones iniciales. Reemplazando en el caso ejemplo, queda: [t,x]=ode45(@dfun,[0 10],[0 1]);

t será el vector de los instantes de tiempo en que se efectuó el cálculo según el criterio de error de este proceso, x será una matriz de dos columnas, de las cuales, la primera irá definiendo la posición y la segunda la velocidad según el vector t. Con lo cual, si se desea mostrar el resultado en un gráfico, se implementará el siguiente código: hold on plot(t,w(:,1),'b') plot(t,w(:,2),'r') xlabel('tiempo') ylabel('Posición y Velocidad') title('Solución EDO') legend('Posición','Velocidad',3) grid on hold off



Y el resultado por pantalla será el siguiente: Solución EDO 2 1.5

Posición y Velocidad

1 0.5 0 -0.5 -1 -1.5 -2

Posición Velocidad 0

1

2

3

4

5 tiempo

6

7

8

9

10



Factores de Amplificación Dinámica (FAD) En el análisis de sistemas dinámicos, un punto que importa conocer es la amplificación que sufre la respuesta del sistema ante solicitaciones con respecto a la respuesta estática. De esta manera, se define el Factor de amplificación dinámico de un sistema, para una cierta razón de amortiguamiento, como el cociente entre la amplitud de la respuesta (para una cierta excitación) y el desplazamiento (velocidad, aceleración) estático. Para poder definir este parámetro, es necesario conocer la frecuencia que posee la excitación a la cual se esta sometiendo el sistema. Conocida esta, es posible analizar lo que ocurre con la respuesta del sistema para diferentes excitaciones, y a su vez, analizar el comportamiento para un determinado amortiguamiento. Para un mejor entendimiento, veamos un ejemplo: Sea la EDM:

mu(t )  cu (t )  ku(t )  F0 sin p  t  Rescrita en forma canónica:

u(t )  2      u (t )   2 u (t ) 

F0 sin  p  t  m

La solución a esta ecuación diferencial de segundo orden, esta dada por una parte homogénea y otra particular. Resolviendo, se obtiene la siguiente solución para un sistema amortiguado de un grado de libertad:

u (t ) 

F0 m 2

1 2

2   p 2   1       2    p       w  

sin  p  t   

   2     p        tan 1  2   1   p         Donde:

 : Razón de amortiguamiento crítico. F0 : Amplitud de la carga.



p : Frecuencia angular de la carga.

u st : Desplazamiento estático del sistema. Luego, el factor de amplificación dinámica para desplazamiento se define numéricamente de la siguiente manera:

FAD u 

1 2

2   p 2   1       2    p       w  

Es posible conocer también, los factores de amplificación dinámica para la velocidad y la aceleración: Diferenciando la respuesta del sistema, se tiene que la velocidad es de la forma:

u (t ) 

 p    

F0 cos  p  t    m  2 2 2  1   p     2    p       w  

De esta manera, el FAD para velocidad es:

FAD v 

 p     2

2   p 2   1       2    p       w  

La aceleración:

u(t ) 

F0 m

 p     2

2

2   p 2   1       2    p       w  

sin  p  t   



De esta manera, el FAD para aceleración es:

FAD a 

 p    

2

2

2   p 2   p  1       2          w  

N   m 

Por ejemplo, sea: m=100 [kg], k  5.2  10 3 

Se realiza un barrido de la frecuencia de la carga desde 0 hasta 100 [rad/s] y un barrido de la razón de amortiguamiento crítico comenzando en 0.1 hasta 1.1. Implementando el siguiente código en MATLAB se obtiene la gráfica que ilustra el FAD para desplazamiento, velocidad y aceleración: %FACTORES DE RESPUESTA DINÁMICA (FAD) %Para una excitación conocida de la forma: F(t)=Fo*sin(p*t) tic clear all clc %Para un oscilador de 1 Grado de Libertad con distintos niveles de amortiguamiento... %Sea m=100;%[kg] k=5.2*10^3;%[N/m] w=(k/m)^(1/2);%[rad/s] p=(0:.5:100);%Frecuencias de la carga aplicada al sistema. beta=[0.1 0.2 0.4 0.5 0.7 0.9 1 1.1];%razones de amortiguamiento aux1=size(p,2); aux2=size(beta,2); FADu=zeros(aux1,aux2); FADv=zeros(aux1,aux2); FADa=zeros(aux1,aux2); for j=1:1:aux2 for i=1:1:aux1 aux3=((1-(p(i)/w)^2)^2+(2*beta(j)*(p(i)/w))^2)^(1/2); FADu(i,j)=1/aux3; FADv(i,j)=(p(i)/w)/aux3; FADa(i,j)=((p(i)/w)^2)/aux3; end end p_w=p/w; figure(1) hold on plot(p_w,FADu(:,1),'-k') plot(p_w,FADu(:,2),'-r')



plot(p_w,FADu(:,4),'-b') plot(p_w,FADu(:,5),'-g') plot(p_w,FADu(:,6),'-m') plot(p_w,FADu(:,7),'-y') legend('\beta=0.1','\beta=0.2','\beta=0.5','\beta=0.7','\beta=0.9','\beta =1') hold off xlabel('p/\omega') ylabel('FAD') title('FAD DESPLAZAMIENTO OSCILADOR SIMPLE AMORTIGUADO-CARGA SINUSOIDAL') grid on figure(2) hold on plot(p_w,FADv(:,1),'-k') plot(p_w,FADv(:,2),'-r') plot(p_w,FADv(:,4),'-b') plot(p_w,FADv(:,5),'-g') plot(p_w,FADv(:,6),'-m') plot(p_w,FADv(:,7),'-y') legend('\beta=0.1','\beta=0.2','\beta=0.5','\beta=0.7','\beta=0.9','\beta =1') hold off xlabel('p/\omega') ylabel('FAD') title('FAD VELOCIDAD OSCILADOR SIMPLE AMORTIGUADO-CARGA SINUSOIDAL') grid on figure(3) hold on plot(p_w,FADa(:,1),'-k') plot(p_w,FADa(:,2),'-r') plot(p_w,FADa(:,4),'-b') plot(p_w,FADa(:,5),'-g') plot(p_w,FADa(:,6),'-m') plot(p_w,FADa(:,7),'-y') legend('\beta=0.1','\beta=0.2','\beta=0.5','\beta=0.7','\beta=0.9','\beta =1') hold off xlabel('p/\omega') ylabel('FAD') title('FAD ACELERACION OSCILADOR SIMPLE AMORTIGUADO-CARGA SINUSOIDAL') grid on toc







Respuesta ante carga “Rampa-Escalón” Considerando que se tiene una excitación como

200  t kN  para t  0.5 s  F (t )    100 kN  para t  0.5 s  Considerando un sistema de 1 GL, se tiene que la EDM es

100  t para t  t d s   mu(t )  cu (t )  ku(t )   t d  100 para t  t d s  La respuesta por tanto debe separarse en 2 tramos. i)

Respuesta en la fase de carga lineal creciente t  t d La solución es

u  u (t )  st e  wt t d  ii)

  2   2   2   2  1    sin  d  t    cos d  t     t      d      

Respuesta en la fase de carga lineal creciente t  t d La solución es

   u (t )     u (td )  ust    sin d  t  td   ust u (t )  e  wt td  u (t d )  ust  cosd  t  td    d d     Donde u (t d ), u (t d ) corresponde a las condiciones iniciales provistas por la fase anterior.

De esta manera, se obtiene la respuesta, que de manera gráfica, permite apreciar que en la fase de carga creciente la respuesta del sistema crece linealmente hasta la posición de equilibrio estático, siendo igual para distintas razones de amortiguamiento. Luego, oscila en torno a ella hasta mantenerse en dicha posición.



%RESPUESTA OSCILACIONES FORZADAS %Para una excitación tipo rampa-escalon m=5000;%[kg] a=.3;%[m] I=(a^4)/12;%[m^4] E=2.94*10^(10);%[N/m^2] K=768*E*I/875;%[N/m] d=0.05; Fo=100000;%[N] i=0; tr=.5; H=@(t)(((Fo*t/tr)*(t>=0 & t<=tr))+(Fo*(t>tr & t<=inf))); ho=@(t)(H(t)/m); for T=0.1:0.01:20; wn=2*pi/T; xst=Fo/(m*wn^2); y=@(t,y)([y(2);ho(t)-2*d*wn*y(2)-wn^2*y(1)]); i=i+1; [t,w]=ode45(y,[0 20],[0 0]); TT(i)=T; xmax(i)=max(abs(w(:,1)))/xst; figure(1) plot(0.5./TT,xmax,'g') title('Espectro de Respuesta, \xi=0.05') xlabel('t_d/T') ylabel('R_d=u_m_á_x/u_s_t') legend('Rampa-Escalón') grid on



Respuestas a partir de un registro de aceleraciones Considerando que se tiene un sistema como el mostrado en la figura, donde u g (t ) representa la aceleración basal discretizada en pasos de t s  , se requiere resolver el problema utilizando el método de Newmark con aceleración lineal (es una alternativa). Esto permite determinar el la relación “tiempo-historia” de la estructura, tanto a nivel de desplazamientos, como de velocidad y aceleración.

La EDM es

mu(t )  u g (t )   cu (t )  ku(t )  0  mu(t )  cu (t )  ku(t )  mu g (t )  u(t )  2      u (t )   2 u (t )  u g (t ) Ahora bien, implementando “Newmark con aceleración lineal”, se tiene

3 6 t     6   K  t  c  t 2  m  u  F  c  3  u (t )  2  u(t )  m   t  u (t )  3  u (t ) Con condiciones iniciales

u (0)  0 u (0)  0 u(0) 

F (0)  c  u (0)  K  u (0) m



De esta forma, para un sistema de período T=0.1 (s) y razón de amortiguamiento de   0.02 , se obtiene el gráfico para el desplazamiento lateral u(t).

%Time-History para SDOF mediante Newmark con aceleración lineal %Se carga el registro de aceleraciones. dato=load('aceleraciones.dat'); dato1=dato*9.81; dt=0.02;%Incremento de tiempo del registro [s]. n=length(dato); tmax=(n-1)*dt; t=(0:dt:tmax)'; d=zeros(n,1); v=zeros(n,1); a=zeros(n,1); %Condiciones iniciales. d(1)=0;%Desplazamiento inicial nulo. v(1)=0;%Velocidad inicial nula. a(1)=(dato1(1)-c*v(1)-K*d(1))/m;%Aceleración inicial . K_eq=K+(3*c/dt)+(6*m/(dt*dt));%Masa equivalente. dd=d(1); vv=v(1); aa=a(1); for i=1:1:(n-1) Apuntes de MatLab para Dinámica de Estructuras – II Semestre 2012 – USM Campus Santiago – RevB // Ricardo Bustamante H.


delta_F=(dato1(i+1)-dato1(i))+c*(3*vv+(dt*aa/2))+m*((6*vv/dt)+3*aa); du=delta_F/K_eq; dv=(3*du/dt)-3*vv-(dt*aa/2); da=(6*du/(dt^2))-(6*vv/dt)-3*aa; d(i+1)=dd+du; v(i+1)=vv+dv; a(i+1)=aa+da; dd=d(i+1); vv=v(i+1); aa=a(i+1); end SD=max(abs(d*100)); SV=max(abs(v*100)); SA=max(abs(a))/9.81; figure(1) subplot(2,1,1); plot(t,dato,'r'); title('Acelerograma'); xlabel('Tiempo [s]'); ylabel('Aceleración [g]'); subplot(2,1,2); plot(t,d*100,'b'); legend('T= 0.1 [s], \xi=0.02') xlabel('Tiempo [s]'); ylabel('Deformación u(t) [cm]');

El espectro de respuesta para un registro se calcula para un barrido de períodos en un oscilador de un grado de libertad y para una razón de amortiguamiento fija. Los espectros se definen de la siguiente manera:

S d  u(t )

max

 2   Sv   S d  T 

 2   Sa    Sd  T  2

%Espectro de Respuesta Lineal para SDOF mediante Newmark con aceleración lineal %Se carga el registro de aceleraciones. dato=load('aceleraciones.dat'); dato1=dato*9.81; dt=0.02;%Incremento de tiempo del registro [s]. n=length(dato); tmax=(n-1)*dt; t=(0:dt:tmax)'; d=zeros(n,1); v=zeros(n,1); a=zeros(n,1);



T=(0.01:0.01:7)'; z=length(T); Sd=zeros(z,1); Sv=zeros(z,1); Sa=zeros(z,1); for j=1:1:z wn=2*pi/T(j); K=m*wn^2;%Rigidez del sistema en [N/m]. c=2*sqrt(m*K)*p;%Viscosidad en [N*s/m]. %Condiciones iniciales. d(1)=0;%Desplazamiento inicial nulo. v(1)=0;%Velocidad inicial nula. a(1)=(dato1(1)-c*v(1)-K*d(1))/m;%Aceleración inicial . K_eq=K+(3*c/dt)+(6*m/(dt*dt));%Masa equivalente. dd=d(1); vv=v(1); aa=a(1); for i=1:1:(n-1) delta_F=(dato1(i+1)-dato1(i))+c*(3*vv+(dt*aa/2))+m*((6*vv/dt)+3*aa); du=delta_F/K_eq; dv=(3*du/dt)-3*vv-(dt*aa/2); da=(6*du/(dt^2))-(6*vv/dt)-3*aa; d(i+1)=dd+du; v(i+1)=vv+dv; a(i+1)=aa+da; dd=d(i+1); vv=v(i+1); aa=a(i+1); end Sd(j,1)=max(abs(d)); Sv(j,1)=(2*pi/T(j))*Sd(j,1); Sa(j,1)=((2*pi/T(j))^2)*Sd(j,1); end figure(1) plot(T,Sd*100) title('Espectro de Respuesta SDOFS \beta=0.02') xlabel('Periodo T_n [s]'); ylabel('Sa [cm]');



Espectro Lineal de Aceleración Obtener el Espectro Lineal de Aceleración, utilizando un algoritmo de Newmark, para una estructura de un grado de libertad, entre los periodos de T = 0 a T = 3 *s+, del registro “n” del terremoto del 27 de Febrero de 2010, considerando los amortiguamientos ξ = *0,a,b+% (n = 1: Constitución E-O; n = 2: Concepción N-S; n = 3: Angol E-O; n=4: Curicó E-O) El “.txt” del registro, contiene un vector con las aceleraciones en *g+, separadas por un intervalo Δt. A continuación, se indican las especificaciones y condiciones iniciales de los registros: Registro

Δt [s]

N° datos

Desplazamiento Inicial Acelerógrafo [cm]

Velocidad Inicial Acelerógrafo [cm/s]

Constitución E-O

0,005

28.657

-0,011

-0,144

Concepción N-S

0,005

28.338

-0,035

1,641

Angol E-O

0,01

18.001

-0,003

0,099

Curicó E-O

0,01

18.001

0,011

0,036

Para obtener el Espectro Lineal de Aceleración se modela un sistema de un grado de libertad de la siguiente forma:

Con la ecuación diferencial de este problema como:

mt   ut   ct   ut   K t   ut   F t 

donde: m: Masa del sistema en el tiempo t en [kg] Apuntes de MatLab para Dinámica de Estructuras – II Semestre 2012 – USM Campus Santiago – RevB // Ricardo Bustamante H.


ü( t): Aceleración del sistema en el tiempo t en [m/s2] c: Viscosidad del sistema en el tiempo t en [N·s/m] u (t): Velocidad del sistema en el tiempo t en [m/s] K: Rigidez del sistema en el tiempo t en [N/m] u(t): Posición del sistema en el tiempo t en [m] F(t): Fuerza externa sobre el sistema en el tiempo t en [N] Si se evalúa el sistema en un tiempo t + Δt:

mt   ut  t   ct   ut  t   K t   ut  t   F t  t  Notar que se consideran las propiedades masa, viscosidad y rigidez, constantes durante el intervalo Δt, tomándose para el análisis, el valor de ellas en el tiempo t (principio del intervalo). Y luego se resta la expresión anterior con la ecuación para el tiempo t, donde:

ut   ut  t   ut  ut   ut  t   ut  ut   ut  t   ut 

F t   F t  t   F t  La ecuación de equilibrio incremental se obtiene para el intervalo de tiempo Δt:

mt   ut   ct   ut   K t   ut   F t  donde: m(t): Masa del sistema al inicio del intervalo en [kg] Δü (t): Incremento de la aceleración del sistema en el intervalo en [m/s2] c(t): Viscosidad del sistema al inicio del intervalo en [N·s/m] Δ u (t): Incremento de la velocidad del sistema en el intervalo en [m/s] K(t): Rigidez del sistema al inicio del intervalo en [N/m] Δu(t): Incremento de la posición del sistema en el intervalo en *m+ ΔF(t): Incremento de la fuerza externa sobre el sistema en el intervalo en *N+ La anterior ecuación será la que se resolverá para cada intervalo Δt, para ella se conocen los siguientes datos provenientes de la resolución del intervalo anterior: m(t), c(t), K(t), ü(t), u (t), u(t); además de ser dato ΔF(t), ya que son datos F(t) y F(t + Δt). Quedando como incógnitas Δü( t ), Δ u (t) y Δu(t), o sea ü(t + Δt), u (t + Δt), u(t + Δt).

En cada paso Δt la solución realizando una integración directa es:

u t  t   u t  

t  t

 u d t



u t  t   u t  

t  t

 u d t

Dado que no se conoce como cambia la aceleración ni la velocidad durante el intervalo, se supone un tipo de variación, usándose para este ejemplo, la de Aceleración Lineal con β = 1/6 del Método de Newmark:

Realizando los reemplazos respectivos:


Se obtiene:



Si de despejan Δ u (t) y Δü(t) en función de Δu(t), quedan las siguientes expresiones:

3 t  ut   3  u t    ut  t 2 6 6 ut   2  ut    u t   3  ut  t t u t  

Reemplazando en la ecuación de equilibrio incremental:

6 t  6  3  mt    2  u t    u t   3  ut   ct     ut   3  u t    ut   K t   ut   F t  t 2  t   t  Ordenando:

3 t  6  6     t 2  mt   t  ct   K t   u t   F t    t  u t   3  ut   mt   3  u t   2  ut   ct  Con:

K t  

6 3  mt    ct   K t  como la rigidez efectiva 2 t t



t 6    F t   F t     u t   3  ut   mt   3  u t    ut   ct  como el incremento de fuerza 2  t    efectiva Se obtiene el valor para Δu(t):

u t  

F t  K t 

Algoritmo Método β = 1/6 de Newmark: 1) Se conocen las condiciones iniciales del problema: m(0), c(0), K(0), u (0), u(0), ΔF(0), con ello se despeja ü(0) de la ecuación de diferencial del problema:

u0 

F 0  c0  u 0  K 0  u 0 m0

2) Una vez conocidas las condiciones iniciales del problema o los datos finales del intervalo anterior: m(t), c(t), K(t), ü(t), u (t), u(t), ΔF(t), se calculan K y F

K t  

6 3  mt    ct   K t  2 t t t 6    F t   F t     u t   3  ut   mt   3  u t    ut   ct  2  t    3) Calcular Δu( t ) y con ello se calculan Δ u (t) y Δü(t)

u t  

F t  K t 


4) Calcular ü( t + Δt ), u ( t + Δt ), u( t + Δt ):

ut  t   ut   ut  ut  t   ut   ut  Apuntes de MatLab para Dinámica de Estructuras – II Semestre 2012 – USM Campus Santiago – RevB // Ricardo Bustamante H.


ut  t   ut   ut  5) 5) Repetir desde el paso 2, ocupando ü(t + Δt), u (t + Δt), u(t + Δt) como los ü(t), u (t), u(t) del siguiente intervalo. El método de Newmark, tiene límites de convergencia y estabilidad, según el Δt y el periodo T de la estructura: Límite de convergencia: t  0,39  T Límite de estabilidad: t  0,55  T

t 0,39 t  T 0,55  T

Dado el Δt en el problema, debería considerarse el límite de convergencia como criterio de elección del periodo más pequeño a evaluar. En el caso del problema la fuerza externa se puede expresar en función de la aceleración del registro üg( t ), que viene con sus datos en [g], por lo cual, la fuerza externa se define como:

F t   m  ug t   9,81 Como el problema es un sistema lineal, las propiedades de éste son constantes en todo momento:

mt   m K t   k

ct   2    k  m , con ξ amortiguamiento del sistema Se debe realizar el Método β = 1/6 de Newmark para cada periodo. Luego se obtiene el mayor valor absoluto de aceleración para cada periodo y con estos datos se construye el gráfico Espectro v/s Periodo de la respuesta. Los espectros de respuesta normalmente entregan el espectro en unidades de [g], por lo cual, se debería dividir la máxima aceleración absoluta, por 9,81. Notar que el Método β = 1/6 de Newmark tiene criterio de convergencia, por lo cual, el mínimo periodo a evaluar, debiese ser T 

t , luego asignar un paso ojalá menor igual a 0,05 [s] en 0,39

incrementación del periodo, hasta evaluar T = 3 [s]. A continuación, se presentan los periodos mínimos a usar para cada registro y las condiciones iniciales del problema:



Registro

Δt [s]

Constitución E-O Concepción N-S Angol E-O Curicó E-O

0,005 0,005 0,01 0,01

T mínimo a usar[s] 0,05 > 0,0128 0,05 > 0,0128 0,05 > 0,0256 0,05 > 0,0256

Desplazamiento Inicial Acelerógrafo [m] -0,00011 -0,00035 -0,00003 0,00011

Velocidad Inicial Acelerógrafo [m/s] -0,00144 0,01641 0,00099 0,00036

Ya que el dato a evaluar es el periodo, el sistema se puede definir con una masa cualquiera y la rigidez del sistema quedará en función de m y T:

 2   k  m   T 

2

Todo lo anterior, hay que hacerlo para ξ = 0%; ξ = a% y ξ = b%. Alguno de los espectros de respuesta se adjuntan a continuación, con ΔT = 0,005: Espectro Registro Constitución Dirección E-O 20

=0  = 0.03  = 0.055

18 16 14

Sa (g)

12 10 8 6 4 2 0

0

0.5

1

1.5 Periodo (s)

2

2.5

3



Espectro Registro Concepción Dirección N-S 8

=0  = 0.03  = 0.055

7 6

Sa (g)

5 4 3 2 1 0

0

0.5

1

1.5 Periodo (s)

2

2.5

3

Espectro Registro Angol Dirección E-O 25

=0  = 0.03  = 0.055

20

Sa (g)

15

10

5

0

0

0.5

1

1.5 Periodo (s)

2

2.5

3



Espectro Registro Curicó Dirección E-O 14

=0  = 0.03  = 0.055

12

10

Sa (g)

8

6

4

2

0 0

0.5

1

1.5 Periodo (s)

2

2.5

3



%ESPECTRO LINEAL DE ACELERACION DF=11;%Dígito Verificador dT=0.005;%Intervalo de periodos %Constitución E-O %dt=0,005 Datos=28657 %di = -0,011 cm %vi = -0,144 cm/s %Angol E-O %dt=0,01 Datos=18001 %di = -0,003 cm %vi = 0,099 cm/s %Concepción N-S %dt=0,005 Datos=28338 %di = -0,035 cm %vi = 1,641 cm/s %Curicó E-O %dt=0,01 Datos=18001 %di = 0,011 cm %vi = 0,036 cm/s %fid=fopen('Constitución.txt'); %registro=fscanf(fid,'%f',[28657,1]); % % % % % %

fid=fopen('Concepción N-S.txt'); registro=fscanf(fid,'%f',[18001,1])/981; fclose(fid); fid=fopen('Concepción.txt','wt'); fprintf(fid,'%f\n',registro); fclose(fid);

registro=load('Constitución.txt'); x=length(registro); color=['b','r','k']; ddd=[0,(DF+5)/200,(DF+10)/200]; m=1000; df=m*(registro(2:x,:)-registro(1:x-1))*9.81; dt=0.005; T=0.05:dT:3; for j=1:length(ddd); k=m*(2*pi./T).^2; c=2*ddd(j)*(2*pi./T)*m; K=k+3/dt*c+6/dt^2*m*ones(1,length(T)); d=zeros(length(registro),length(T)); v=d; a=v; d(1,:)=-0.011/100; v(1,:)=-0.144/100; a(1,:)=(m*registro(1)*ones(1,length(T))*9.81-c.*v(1,:)-k.*d(1,:))/m; for i=1:x-1;



dF=df(i)*ones(1,length(T))+c.*(3*v(i,:)+dt/2*a(i,:))+m*(6/dt*v(i,:) +3*a(i,:)); dd=dF./K; dv=3/dt*dd-3*v(i,:)-dt/2*a(i,:); da=6/dt^2*dd-6/dt*v(i,:)-3*a(i,:); a(i+1,:)=a(i,:)+da; v(i+1,:)=v(i,:)+dv; d(i+1,:)=d(i,:)+dd; end Sa=max(abs(a)); figure(1) hold on plot(T,Sa/9.81,color(j)) end ylabel('Sa (g)') xlabel('Periodo (s)') title('Espectro Registro Constitución Dirección E-O') legend('\xi = 0',char(['\xi = ',num2str(ddd(2))]),char(['\xi = ',num2str(ddd(3))]),1)



Espectro No Lineal de Aceleración Se modela el mismo sistema del ejemplo anterior, consistente de un grado de libertad de la siguiente forma:

La ecuación diferencial de este problema es la siguiente:

mt   ut   ct   ut   f s t   F t 

Donde: m: Masa del sistema en el tiempo t en [kg] ü( t ): Aceleración del sistema en el tiempo t en [m/s2] c: Viscosidad del sistema en el tiempo t en [N·s/m] u ( t ): Velocidad del sistema en el tiempo t en [m/s] fs( t ): Fuerza interna del sistema en el tiempo t en [N] F( t ): Fuerza externa sobre el sistema en el tiempo t en [N] La fuerza interna se define:

u t   k  f s t    f y f y 

si u t   k  f y si u t   k  f y si u t   k   f y

donde: k: Rigidez del sistema en [N/m] fy: Tensión de fluencia en [N] Si se evalúa el sistema en un tiempo t + Δt:

mt   ut  t   ct   ut  t   f s t  t   F t  t 

* Se consideran las propiedades masa, viscosidad y rigidez, constantes durante el intervalo Δt, tomándose para el análisis, el valor de ellas en el tiempo t (principio del intervalo). Y luego se resta la expresión anterior con la ecuación para el tiempo t:

donde:



ut   ut  t   ut  ut   ut  t   ut 

K t   ut   f s t  t   f s t 

F t   F t  t   F t 

Se obtiene la ecuación de equilibrio incremental para el intervalo de tiempo Δt:

mt   ut   ct   ut   K t   ut   F t 

con: m(t): Masa del sistema al inicio del intervalo en [kg] Δü(t): Incremento de la aceleración del sistema en el intervalo en *m/s2] c(t): Viscosidad del sistema al inicio del intervalo en [N·s/m] Δ u (t): Incremento de la velocidad del sistema en el intervalo en [m/s] K(t): Rigidez según ley elastoplástica en el instante t en [N/m] Δu(t): Incremento de la posición del sistema en el intervalo en *m+ ΔF(t): Incremento de la fuerza externa sobre el sistema en el intervalo en [N] La anterior ecuación será la que se resolverá para cada intervalo Δt, para ella se conocen los siguientes datos provenientes de la resolución del intervalo anterior: m(t), c(t), K(t), ü(t), u (t), u(t); además de ser dato ΔF(t), ya que son datos F(t) y F(t + Δt). Quedando como incógnitas Δü( t ), Δ u (t) y Δu(t), o sea ü(t + Δt), u (t + Δt), u(t + Δt). En cada paso Δt, la solución realizando una integración directa es:

u t  t   u t  

t  t

 u d t

u t  t   u t  

t  t

 u d t

Dado que no se conoce como cambia la aceleración ni la velocidad durante el intervalo, se supone un tipo de variación, usándose para este ejemplo, la de Aceleración Lineal con β = 1/6 del Método de Newmark:



Realizando los reemplazos respectivos:


Se obtiene:

Si de despejan Δ u (t) y Δü t) en función de Δu (t), quedan las siguientes expresiones:




Reemplazando en la ecuación de equilibrio incremental:

6 t  6  3  mt    2  u t    u t   3  ut   ct     ut   3  u t    ut   K t   ut   F t  t 2  t   t  Ordenando:

3 t  6  6     t 2  mt   t  ct   K t   u t   F t    t  u t   3  ut   mt   3  u t   2  ut   ct 

Con:

K t  

6 3  mt    ct   K t  ; Rigidez efectiva 2 t t t 6    F t   F t     u t   3  ut   mt   3  u t    ut   ct  ; Incremento de fuerza efectivo 2  t    Se obtiene el valor para Δu( t ):

u t  

F t  K t 

Algoritmo Método β = 1/6 de Newmark: 1) Se conocen las condiciones iniciales del problema: m( 0 ), c( 0 ), K( 0 ), u ( 0 ), u( 0 ), ΔF( 0 ), con ello se despeja ü( 0 ) de la ecuación de diferencial del problema, también con ello, se obtienen la fuerza interna fs:

u0 

F 0  c0  u 0  K 0  u 0 m0 f s 0  K 0  u0

2) Una vez conocidas las condiciones iniciales del problema o los datos finales del intervalo anterior: m( t ), c( t ), K( t ), ü( t ), u ( t ), u( t ), ΔF( t ), fs( t - Δt ), se calculan fs( t ), F , K en el siguiente orden:

t 6    F t   F t     u t   3  ut   mt   3  u t    ut   ct  2  t   

k si K t    0 si

f s t  t   u t  t   k  f y

f s t  t   u t  t   k  f y



 f t  t   u t  t   k  f s t    fy   fy  K t  

si si si

f s t  t   u t  t   k  f y f s t  t   u t  t   k  f y f s t  t   u t  t   k   f y

6 3  mt    ct   K t  2 t t

3) Calcular Δu( t ) y con ello se calculan Δ u ( t ) y Δü( t ).

u t  

F t  K t 


4) Calcular ü( t + Δt ), u ( t + Δt ), u( t + Δt ):

ut  t   ut   ut  ut  t   ut   ut  ut  t   ut   ut 

5) Repetir desde el paso 2, ocupando ü(t + Δt), u (t + Δt), u(t + Δt), fs(t) como los ü(t), u (t), u(t), fs(t - Δt) del siguiente intervalo. El método de Newmark, tiene límites de convergencia y estabilidad, según el Δt y el periodo T de la estructura: Límite de convergencia: t  0,39  T Límite de estabilidad: t  0,55  T

t 0,39 t  T 0,55  T

Dado el Δt en el problema, debería considerarse el límite de convergencia como criterio de elección del periodo más pequeño a evaluar. En el caso del ejemplo, la fuerza externa se puede expresar en función de la aceleración del registro üg (t), que viene con sus datos en [g], por lo cual, la fuerza externa se define como:

F t   m  ug t   9,81

Al ser este un sistema elastoplástico, las propiedades del sistema son las siguientes:

mt   m



k si K t    0 si

f s t  t   u t  t   k  f y

f s t  t   u t  t   k  f y ct   2    k  m

 f t  t   u t  t   k  f s t    fy   fy 

si si si

f s t  t   u t  t   k  f y f s t  t   u t  t   k  f y f s t  t   u t  t   k   f y

Con ξ Amortiguamiento del sistema y fy, la tensión de fluencia:

Dado que se pide el espectro no lineal de aceleración, se debe realizar el Método β = 1/6 de Newmark para cada periodo, luego se obtiene el mayor valor absoluto de aceleración para cada periodo, con estos datos, se construye el gráfico Espectro v/s Periodo de la respuesta. Los espectros de respuesta normalmente entregan el espectro en unidades de [g], por lo cual, se debería dividir la máxima aceleración absoluta, por 9,81. Notar que el Método β = 1/6 de Newmark tiene criterio de convergencia, por lo cual, el mínimo periodo a evaluar, debiese ser T 

t , luego asignar un paso ojalá menor igual a 0,05 [s] en 0,39

incrementación del periodo, hasta evaluar T = 3 [s]. A continuación, se presentan los periodos mínimos a usar para cada registro y las condiciones iniciales del problema:



Registro

Δt [s]

Constitución E-O Concepción N-S Angol E-O Curicó E-O

0,005 0,005 0,01 0,01

T mínimo a usar[s] 0,05 > 0,0128 0,05 > 0,0128 0,05 > 0,0256 0,05 > 0,0256

Desplazamiento Inicial Acelerógrafo [m] -0,00011 -0,00035 -0,00003 0,00011

Velocidad Inicial Acelerógrafo [m/s] -0,00144 0,01641 0,00099 0,00036

Ya que el dato a evaluar es el periodo, el sistema se puede definir con una masa cualquiera y la rigidez del sistema quedará en función de m y T:

 2   k  m   T 

2

La obtención del espectro para el ejemplo se realiza para cada tensión de fluencia desde fy = 0,2 * W, con un paso de B * W, hasta la tensión de fluencia con la que el sistema se comporte en forma lineal, es decir, hasta que en todos los puntos evaluados fs < fy. El resultado para el registro de Constitución se presenta en el gráfico adjunto. No se parte de períodos bajos puesto que la respuesta para intervalos de mucha rigidez distorsiona la visualización de la respuesta para el rango de análisis usual.



%ESPECTRO LINEAL DE ACELERACION %Resetea todas las variables DF=1; %Dígito Verificador dT=0.1;%Intervalo de periodos ddd=2*DF/100; %Amortiguamiento dfy=0.2; %Paso deltaFy = dfy*W %Se ingresa vector con distintos colores dados por las propiedades RGB, que no están definidos por defecto con una letra en MATLAB, con el fin de graficar colores=[0.5 0.5 0.5;0.5 0 0.5;0 0.5 0.5;1 0.5 0;0.5 0.25 0;0.5 0 0;0 0.5 0]; %Se ingresa vector con los colores definidos por una letra en MATLAB, para graficar color=['y';'b';'k';'m';'c';'r';'g']; m=1000; %Masa del sistema (Este puede tomar cualquier valor >= 0 ) T=0.4:dT:3; %Rango de periodos para evaluar el espectro k=m*(2*pi./T).^2; %Cálculo de la rigidez para cada periodo c=2*ddd*(2*pi./T)*m; %Cálculo de la viscosidad para cada periodo %Carga de los datos del acelerograma del registro de Constitución registro=load('Constitución.txt'); %El registro se encuentra en [g], se pasa a [m/s^2] registro=registro*9.81; %Se mide el largo de datos del registro x=length(registro); %Se obtiene el incremento de fuerza del registro, restando cada punto del registro con su anterior dato, y todo ello multiplicado por la masa df=m*(-registro(2:x,:)+registro(1:x-1)); %deltaF = - m * üg(t) dt=0.005; %Paso de tiempo del registro para Constitución %Cálculo de la rigidez efectiva, se crea un vector del mismo tamaño que la cantidad de periodos a evaluar K=k+3/dt*c+6/dt^2*m*ones(1,length(T)); %Valores iniciales para ingresar al while fs=0.2*m*9.8; %Valor dado para poder ingresar al while fy=0.2*m*9.8*ones(1,length(T)); %Fluencia inicial de 0,2*W j=0;%Iniciación del contador %Se crea un proceso repetitivo while, que seguirá funcionando mientras algún fs alcance el valor de fluencia en algún punto del registro de la respuesta



while sum(max(abs(fs))~=fy)~=length(T); %Solo se saldrá del while, si todos los valores de fs en valor absoluto, son menores a la fluencia evaluada j=j+1; %Contador %Para la primera iteración, la fluencia fue iniciada afuera del while if j~=1; fy=fy+dfy*m*9.8*ones(1,length(T)); %Incrementar Fy en 0,2*W end %Iniciación de matrices de posición, velocidad, aceleración y fs, en que cada columna representa un periodo y cada fila representa un instante de tiempo de la respuesta d=zeros(length(registro),length(T)); v=d; a=v; fs=a; %Condiciones iniciales del acelerógrafo d(1,:)=-0.011/100; v(1,:)=-0.144/100; %Cálculo de la aceleración inicial del acelerógrafo a(1,:)=(-m*registro(1)*ones(1,length(T))-c.*v(1,:)-k.*d(1,:))/m; %fs inicial fs(1,:)=0; %Primer fs con las condiciones iniciales fs(2,:)=k.*d(1,:); %Newmark inelástico for i=1:x-1;%Irá recorriendo cada instante delta t %Cálculo del incremento de fuerza dF=df(i)+c.*(3*v(i,:)+dt/2*a(i,:))+m*(6/dt*v(i,:)+3*a(i,:)); %El fs inicial es dado según las condiciones iniciales, fuera del for, los siguientes se calculan, según las condiciones del paso anterior if i~=1; %Función por tramos, si fs + k*du en valor absoluto, es menor a la fluencia, el nuevo fs será fs + k*du, sino, tomará el valor de fy si fs + k*du es positivo o -fy si es negativo fs(i+1,:)=(fs(i,:)+k.*dd).*(abs(fs(i,:)+k.*dd)<=fy(j))+fy(j)* (fs(i,:)+k.*dd>fy(j))-fy(j)*(fs(i,:)+k.*dd<-fy(j)); end %Cálculo del incremento de u, o sea, du dd=dF./(K-k.*(abs(fs(i+1,:))==fy(j))); %Cálculo del incremento de v, o sea, dv dv=3/dt*dd-3*v(i,:)-dt/2*a(i,:); %Cálculo del incremento de a, o sea, da da=6/dt^2*dd-6/dt*v(i,:)-3*a(i,:); %Cálculo de las condiciones para el siguiente paso a(i+1,:)=a(i,:)+da; v(i+1,:)=v(i,:)+dv; d(i+1,:)=d(i,:)+dd;



end %Se obtiene el máximo valor de posición en valor absoluto para cada periodo Sd=max(abs(d)); %Cálculo del Pseudo Espectro de Aceleración PSa=Sd.*((2*pi)./T).^2; %Gráfica PSa v/s T para cada fluencia, de tal manera que cada curva, quede con un %color distinto y/o grosor de línea distinto. Notar que los espectros %se dividen por g, para que queden en unidades de [g] figure(1) if j<=6; plot(T,PSa/9.81,color(j+1))% Curvas con líneas de colores definidos en MATLAB y con grosor de línea estándar = 1 elseif j>=14 && j<=20; plot(T,PSa/9.81,color(j-13),'LineWidth',2)% Curvas con líneas de colores definidos en MATLAB y con grosor de línea = 2 elseif j>=7 && j<=13; plot(T,PSa/9.81,'Color',colores(j-6,:))% Curvas con líneas de colores definidos por usuario y con grosor de línea estándar = 1 elseif j>=21 && j<=27; plot(T,PSa/9.81,'Color',colores(j-20,:),'LineWidth',2)% Curvas con líneas de colores definidos por usuario y con grosor de línea = 2 elseif j>=28 && j<=34; plot(T,PSa/9.81,color(j-27),'LineWidth',2.5)% Curvas con líneas de colores definidos en MATLAB y con grosor de línea = 2,5 elseif j>=35 && j<=41; plot(T,PSa/9.81,'Color',colores(j-34,:),'LineWidth',2.5)% Curvas con líneas de colores definidos por usuario y con grosor de línea = 2,5 elseif j>=42 && j<=48; plot(T,PSa/9.81,color(j-41),'LineWidth',3)% Curvas con líneas de colores definidos en MATLAB y con grosor de línea = 3 else plot(T,PSa/9.81,'Color',colores(j-48,:),'LineWidth',3)% Curvas con líneas de colores definidos por usuario y con grosor de línea = 3 end %Comando para que las gráficas se superpongan en el mismo gráfico y no se borren cada vez que se grafica encima hold on end % Se introduce reseña para el eje y ylabel('PSa (g)') % Se introduce reseña para el eje x xlabel('Periodo (s)') % Se introduce el título, acá se pueden pasar variables numéricas a texto: % Ej: char(['texto',num2str(variable numérica),'texto']) title(char(['Pseudoespectro Constitución Dirección E-O para distintas tensiones de fluencia, con \xi = ',num2str(ddd*100),'%'])) %Leyenda que determina formato de línea para cada fluencia legend(char(['Fy/W = ',num2str(0.2)]),char(['Fy/W = ',num2str(0.2+1*dfy)]),char(['Fy/W = ',num2str(0.2+2*dfy)]),char(['Fy/W = ',num2str(0.2+3*dfy)]),char(['Fy/W = ',num2str(0.2+4*dfy)]),char(['Fy/W = ',num2str(0.2+5*dfy)]),char(['Fy/W = ',num2str(0.2+6*dfy)]),char(['Fy/W =



',num2str(0.2+7*dfy)]),char(['Fy/W = ',num2str(0.2+8*dfy)]),char(['Fy/W = ',num2str(0.2+9*dfy)]),char(['Fy/W = ',num2str(0.2+10*dfy)]),char(['Fy/W = ',num2str(0.2+11*dfy)]),char(['Fy/W = ',num2str(0.2+12*dfy)]),char(['Fy/W = ',num2str(0.2+13*dfy)]),char(['Fy/W = ',num2str(0.2+14*dfy)]),char(['Fy/W = ',num2str(0.2+15*dfy)]),char(['Fy/W = ',num2str(0.2+16*dfy)]),char(['Fy/W = ',num2str(0.2+17*dfy)]),char(['Fy/W = ',num2str(0.2+18*dfy)]),char(['Fy/W = ',num2str(0.2+19*dfy)]),char(['Fy/W = ',num2str(0.2+20*dfy)]),char(['Fy/W = ',num2str(0.2+21*dfy)]),char(['Fy/W = ',num2str(0.2+22*dfy)]),char(['Fy/W = ',num2str(0.2+23*dfy)]),char(['Fy/W = ',num2str(0.2+24*dfy)]),char(['Fy/W = ',num2str(0.2+25*dfy)]),char(['Fy/W = ',num2str(0.2+26*dfy)]),char(['Fy/W = ',num2str(0.2+27*dfy)]),char(['Fy/W = ',num2str(0.2+28*dfy)]),char(['Fy/W = ',num2str(0.2+29*dfy)]),char(['Fy/W = ',num2str(0.2+30*dfy)]),char(['Fy/W = ',num2str(0.2+31*dfy)]),char(['Fy/W = ',num2str(0.2+32*dfy)]),char(['Fy/W = ',num2str(0.2+33*dfy)]),char(['Fy/W = ',num2str(0.2+34*dfy)]),char(['Fy/W = ',num2str(0.2+35*dfy)]),char(['Fy/W = ',num2str(0.2+36*dfy)]),char(['Fy/W = ',num2str(0.2+37*dfy)]),char(['Fy/W = ',num2str(0.2+38*dfy)]),char(['Fy/W = ',num2str(0.2+39*dfy)]),char(['Fy/W = ',num2str(0.2+40*dfy)]),char(['Fy/W = ',num2str(0.2+41*dfy)]),char(['Fy/W = ',num2str(0.2+42*dfy)]),char(['Fy/W = ',num2str(0.2+43*dfy)]),char(['Fy/W = ',num2str(0.2+44*dfy)]),char(['Fy/W = ',num2str(0.2+45*dfy)]),char(['Fy/W = ',num2str(0.2+46*dfy)]),char(['Fy/W = ',num2str(0.2+47*dfy)]),char(['Fy/W = ',num2str(0.2+48*dfy)]),char(['Fy/W = ',num2str(0.2+49*dfy)]),char(['Fy/W = ',num2str(0.2+50*dfy)]),char(['Fy/W = ',num2str(0.2+51*dfy)]),char(['Fy/W = ',num2str(0.2+52*dfy)]),char(['Fy/W = ',num2str(0.2+53*dfy)]),char(['Fy/W = ',num2str(0.2+54*dfy)]),1)



Superposición modal Plantear las ecuaciones del movimiento del marco de acero de un vano y cinco pisos analizado en clase, pero considerando las vigas flexuralmente deformables y columnas axial-flexuralmente deformables; suponer que todas las columnas tienen el mismo momento de inercia I=9230cm4 y área 62cm2, y las vigas un momento de inercia I=12250cm4 mayor, que el primer piso es de 4m de altura y los restantes de 3.5m de altura; y que el vano es de 7.5m de luz. Las masas de los pisos, excepto el último, son todas iguales, la del último es solo el 65% de una de ellas, y se supone que están concentradas a nivel de cada piso. Plantear la ecuación condensando los grados de libertad no asociados a masas (los grados de libertad inerciales son solamente los 5 desplazamientos laterales de los pisos). a) Determinar las formas modales, y calibrar la razón entre rigidez y masas de modo que el primer modo tenga un período 0.2s. Determinar la matriz de amortiguamiento de manera que este sea clásico, con 5% de amortiguamiento en todos los modos. b) Obtener la respuesta a la excitación de la componente ux del registro de Melipilla del terremoto de 1985 (obtenerla del archivo Melip.m) a) integrando directamente con ode23 b) integrando directamente con aceleración constante c) mediante superposición modal con 2, 3 y 5 modos (integrando cada modo independientemente en la forma que quiera) c) Determinar la estimación CQC del máximo momento flector en la base de columna del primer piso, el máximo del esfuerzo axial de la columna, y el máximo de una combinación lineal del momento y el esfuerzo axial (los mismos valores de los coeficientes de la combinación). d) Integrar la ecuación:

y  2 y   2 y  ug considerando u g como el ux del registro de Melipilla del terremoto de 1985 del archivo rec_melip.mat, y para los 5 valores de correspondientes a las frecuencias modales del marco de un vano y cinco pisos y obtener de esas respuestas los máximos que deberían ser valores de las correspondientes ordenadas espectrales. Aplicar superposición modal espectral CQC para estimar los máximos de: a) El desplazamiento del primer piso b) El desplazamiento del quinto piso c) El desplazamiento relativo de entrepiso del tercer piso d) El corte basal e) El corte del tercer entrepiso Comparar con los resultados de usar las coordenadas modales obtenidas para calcular los máximos “exactos”.



a) Matriz de Rigidez Se requiere obtener la matriz de rigidez de un edificio de 5 pisos compuesto por columnas y vigas. Las vigas se considerarán flexurablemente deformables y las columnas axial-flexurablemente deformables. De esta manera los grados axiales de las vigas no se consideraran, utilizando una matriz de rigidez reducida para vigas en posición horizontal que no incluye los términos axiales, mientras que la matriz de rigidez de cada columna contiene todos los términos. Para mayor facilidad en el desarrollo se consideraran primero los grados de libertad que representan el desplazamiento lateral de cada piso, para luego considerar el resto. De esta forma los primero cinco términos de la matriz de rigidez resultante serán los que se usarán para determinar la respuesta a nivel de piso debido al registro. La siguiente figura indica los grados de libertad globales considerados para el edificio de 5 pisos.

Vigas La matriz de rigidez local se obtiene no considera los grados de libertad axiales, por lo que el elemento queda como lo indica la figura:

Considerando la matriz elemental para una viga horizontal y tomando en cuenta sólo los 4 grados de libertad locales indicados, la matriz de rigidez local del elemento viga se expresa como



6 3L   6 3L 3L 2 L2 3L L2  EI  kv  2 c 3 v   6 3L  L  6 3L   2 3L 2 L2  3L L

Columnas En relación a la matriz de rigidez de la columna, se consideran todos los grados de libertad de la matriz de rigidez elemental. De esta forma, se tiene que el elemento columna se expresa en coordenadas locales como:

Considerando los grados de libertad locales de la columna, la matriz de rigidez local del elemento columna se expresa como  Ec I c  12 h 3 i   0    Ec I c  6 2 hi kc    Ec I c  12 3 hi   0    EI  6 c 2c hi 

0

6

Ec I c hi 2

Ec Ac hi



0

0



Ec Ac hi

Ec I c hi 2

0

Ec I c hi 2

12

Ec I c hi 3

0

0

6

0 2

0

6

4

Ec Ac hi

Ec I c hi 3

Ec I c hi

0

0

12

Ec I c hi

0 6

Ec I c hi 2

Ec Ac hi 0

Ec I c  hi 2    0   EI  2 c c  hi  Ec I c  6 2  hi   0   EI  4 c c  hi  6

Se debe hacer la diferencia sobre la altura del primer piso con la del resto de los pisos, para lo cual se generarán matrices de rigidez para el primer piso y para el resto. Todas las vigas tienen las mismas propiedades, por lo cual no se hará diferencia con las matrices de rigidez de éstas. Los valores a utilizar son



 tonf Ec  2038.902  2  cm

   29000  ksi 

Ac  62 cm2  I c  9230 cm4 

I v  12250 cm4 

Lv  750 cm

  400  cm , para el primer piso hi    350  cm , para el resto delos pisos

i  0.05 Una vez realizado lo anterior, es necesario establecer la relación entre los grados de libertad locales y globales. Para ello se enumeran los elementos, con la finalidad de establecer los vectores de conectividad que permitan ensamblar la matriz de rigidez del sistema. La figura que se indica a continuación establece la numeración de los elementos. Se requieren 5 vectores de conectividad para las 5 vigas y 10 vectores de conectividad para las 10 columnas.

El procedimiento descrito a continuación se puede estandarizar utilizando indexación, como aparece en el Capítulo de Análisis Estructural del tutorial, página 23. Resolución respuestas a la excitación de la componente ‘ux’ del registro de Melipilla



% Matrices de rigidez barras % vigas kv=2*Ec*Iv/L^3*[6 3*L -6 3*L; 3*L 2*L^2 -3*L L^2; -6 -3*L 6 -3*L; 3*L L^2 -3*L 2*L^2]; % columnas EA=Ec*Ac; EI=Ec*Ic; Kcp1 = [12*EI/hp1^3,0,-6*EI/hp1^2,-12*EI/hp1^3,0,-6*EI/hp1^2; 0,EA/hp1,0,0,-EA/hp1,0; -6*EI/hp1^2,0,4*EI/hp1,6*EI/hp1^2,0,2*EI/hp1; -12*EI/hp1^3,0,6*EI/hp1^2,12*EI/hp1^3,0,6*EI/hp1^2; 0,-EA/hp1,0,0,EA/hp1,0; -6*EI/hp1^2,0,2*EI/hp1,6*EI/hp1^2,0,4*EI/hp1]; Kcp2 = [12*EI/hp2^3,0,-6*EI/hp2^2,-12*EI/hp2^3,0,-6*EI/hp2^2; 0,EA/hp2,0,0,-EA/hp2,0; -6*EI/hp2^2,0,4*EI/hp2,6*EI/hp2^2,0,2*EI/hp2; -12*EI/hp2^3,0,6*EI/hp2^2,12*EI/hp2^3,0,6*EI/hp2^2; 0,-EA/hp2,0,0,EA/hp2,0; -6*EI/hp2^2,0,2*EI/hp2,6*EI/hp2^2,0,4*EI/hp2]; % Matrices de compatibilidad para cada barra % vigas a1=zeros(4,25); a1(1,6)=1; a1(2,7)=1; a1(3,8)=1; a1(4,9)=1; a2=zeros(4,25); a2(1,10)=1; a2(2,11)=1; a2(3,12)=1; a2(4,13)=1; a3=zeros(4,25); a3(1,14)=1; a3(2,15)=1; a3(3,16)=1; a3(4,17)=1; a4=zeros(4,25); a4(1,18)=1; a4(2,19)=1; a4(3,20)=1; a4(4,21)=1; a5=zeros(4,25); a5(1,22)=1; a5(2,23)=1; a5(3,24)=1; a5(4,25)=1; % columnas a6=zeros(6,25); a6(4,1)=-1; a6(5,6)=1; a6(6,7)=1; a7=zeros(6,25); a7(4,1)=-1; a7(5,8)=1; a7(6,9)=1; a8=zeros(6,25); a8(1,1)=-1; a8(2,6)=1; a8(3,7)=1; a8(4,2)=-1; a8(5,10)=1; a8(6,11)=1; a9=zeros(6,25); a9(1,1)=-1; a9(2,8)=1; a9(3,9)=1; a9(4,2)=-1; a9(5,12)=1; a9(6,13)=1; a10=zeros(6,25); a10(1,2)=-1; a10(2,10)=1; a10(3,11)=1; a10(4,3)=-1; a10(5,14)=1; a10(6,15)=1; a11=zeros(6,25); a11(1,2)=-1; a11(2,12)=1; a11(3,13)=1; a11(4,3)=-1; a11(5,16)=1; a11(6,17)=1; a12=zeros(6,25); a12(1,3)=-1; a12(2,14)=1; a12(3,15)=1; a12(4,4)=-1; a12(5,18)=1; a12(6,19)=1; a13=zeros(6,25); a13(1,3)=-1; a13(2,16)=1; a13(3,17)=1; a13(4,4)=-1; a13(5,20)=1; a13(6,21)=1; a14=zeros(6,25); a14(1,4)=-1; a14(2,18)=1; a14(3,19)=1; a14(4,5)=-1; a14(5,22)=1; a14(6,23)=1; a15=zeros(6,25); a15(1,4)=-1; a15(2,20)=1; a15(3,21)=1; a15(4,5)=-1; a15(5,24)=1; a15(6,25)=1; % Matriz de rigidez Kv=a1'*kv*a1 + a2'*kv*a2 + a3'*kv*a3 + a4'*kv*a4 + a5'*kv*a5; %aporte de las vigas Kc=a6'*Kcp1*a6 + a7'*Kcp1*a7 + a8'*Kcp2*a8 + a9'*Kcp2*a9 + a10'*Kcp2*a10 + a11'*Kcp2*a11 + a12'*Kcp2*a12 + a13'*Kcp2*a13 + a14'*Kcp2*a14 + a15'*Kcp2*a15; %aporte de las columnas Kt=Kv+Kc; %matriz de rigidez completa



Ahora bien, se debe condensar estáticamente la matriz para reducirla a los grados de libertad asociados al desplazamiento horizontal del modelo. La idea es que los grados de libertad dentro de la matriz de rigidez se pueden separar entre sí. En este caso se ‘aislarán’ los asociados al desplazamiento lateral del modelo.


K ab  u   F     Kbb     0 

 Kaa   u   K ab      F   Kba   u   Kbb      0 La rigidez lateral queda expresada como

 K    Kaa    Kab    Kbb   Kba  1

K= 14.7948 -9.8816 -9.8816 15.4632 2.6659 -10.0323 -0.4655 2.6500 0.0756 -0.3583

2.6659 -0.4655 -10.0323 2.6500 15.4597 -9.8125 -9.8125 14.1538 2.0987 -6.5904

0.0756 -0.3583 2.0987 -6.5904 4.7790

a=[1:5]; %grados de libertad horizontales b=[6:15]; %resto de los grados de libertad del sistema Kaa=Kt(a,a); Kab=Kt(a,b); Kba=Kt(b,a); Kbb=Kt(b,b); K=Kaa-Kab*inv(Kbb)*Kba; %matriz de rigidez condensada a grados de libertad laterales con masa asociada

Matriz de Masa Se requiere determinar la matriz de masa tal que el período fundamental sea igual a T  2  seg  . Las masas están concentradas a nivel de piso y son iguales a ‘m’, salvo en el techo, en la cual es igual a ‘0.65m’. De esta forma la matriz de masa es

1 0   M   m 0  0 0

0 0 0 1 0 0 0 1 0 0 0 1 0 0 0

0  0  0   0  0.65



La ecuación de movimiento en coordenadas generalizadas permite determinar las frecuencias del edificio. Resolviendo la ecuación, se tiene ..

 M  q   K  q  0 , con

q  f (t )  v

..

 M  f (t ) v   K  f (t )  v  0 ..

f (t )  M   v   K   v  0 f (t ) m 2 G  f (t )  y  f (t )  K  y  0 m 2 y  G 

1

K  y

Los valores propios aparecen multiplicados por

m al resolver la ecuación m 2 y  G 

1

K  y .

Se debe imponer que el periodo sea T  2  seg  , lo que en términos de la frecuencia equivale a

w  10 . Es decir, al resolver para G 

1

 K  , el valor de la primera frecuencia es w1  10

m.

Como se busca determinar la masa, se debe resolver para m=1 y reemplazar w1 resultante en la expresión determinada con anterioridad. 2

 w1    m  10  % Matriz M G=eye(5); G(5,5)=0.65; [V,D]=eig(G\K); w=sqrt(diag(D)); T=2*pi./w; [T,p]=sortT(T); V=V(:,p); w=w(p); m=(0.1*w(1)/pi)^2; M=m*G; G(5,5)=1/0.65; invM=G/m;



-4 Lo anterior permite determinar la masa, con m  2.6711·10  tonf



s2   cm 



M= 1.0e-003 * 0.2671 0 0 0 0 0 0.2671 0 0 0 0 0 0.2671 0 0 0 0 0 0.2671 0 0 0 0 0 0.1736

Matriz de Amortiguamiento

Se sabe que

C *  V T  C V I  V T  M V

Lo anterior permite determinar la matriz de amortiguamiento

C  M V  C* V T  M Con

C*  2  i  i

Resolución de frecuencias y modos de vibración Teniendo todo lo anterior se procede a resolver para las frecuencias y modos de vibración del edificio. De esta forma, se tiene: [V,D]=eig(K,M); w=sqrt(diag(D)); T=2*pi./w; [T,p] = sortT(T); V=V(:,p); w=w(p); Mm=V'*M*V; for i=1:5 norm=(1/Mm(i,i))^0.5; V(:,i)=norm*V(:,i); end Km=V'*K*V; Cm=diag(2*bi*w); C=M*V*Cm*V'*M; e1=ones(5,1);v Masapormodo=(V'*M*e1).*(V'*M*e1)/(e1'*M*e1);



T= 0.2000 0.0625 0.0342 0.0227 0.0174 w= 31.4159 100.5183 183.8185 276.7465 360.8302 V= 9.2713 20.2458 29.7331 36.4124 39.9604

-25.8557 -37.3297 -19.7283 14.2420 40.9438

35.7330 13.1592 -32.5955 -16.9718 38.0930

-34.9753 25.0694 11.7072 -34.9833 28.5840

22.1162 -33.7443 35.6501 -26.5169 14.7720

Mm = 1.0e-003 * 0.6262 -0.0000 -0.0000 0.0000 0.0000

-0.0000 0.5965 0.0000 0.0000 -0.0000

-0.0000 0.0000 0.6891 0.0000 -0.0000

0.0000 0.0000 0.0000 0.8171 -0.0000

0.0000 -0.0000 -0.0000 -0.0000 0.7868

Km = 1.0e+005 * 0.0099 -0.0000 0. 0000 0.0000 0.0000

-0.0000 0.1010 0.0000 0.0000 -0.0000

-0.0000 0.0000 0.3379 0.0000 -0.0000

0.0000 0.0000 0.0000 0.7659 -0.0000

0.0000 -0.0000 -0.0000 -0.0000 1.3020

Masapormodo = 0.8499 0.1016 0.0333 0.0122 0.0029



b) Respuesta del registro a) Integrando directamente con ‘ode23’

Para la resolución mediante la función ‘ode23’ se utiliza el procedimiento visto en clases, consistente en reducir el grado de la ecuación diferencial para luego, de manera matricial, resolver el problema.

 x1 (t )   q(t )     q(t )  x ( t )   2  

 x(t )   tal que

x1 (t )    M   C   x1 (t )   M    K   x2 (t ) 1

1

Se re-escribe la ecuación como

x(t )   A x  b   M 1  C  M 1  K  ug  e      x   I 0  0   

 x(t )  

A=[-invM*C -invM*K; eye(5) zeros(5)]; B=[-ones(5,1) ; zeros(5,1)]; xo=zeros(10,1); t1=0; t2=160; %cond. iniciales y tiempo integración registro tsp=[t1; t2]; [t_ode,x_ode]=ode23('shearbuilding',tsp,xo);

La funcion Shearbuilding corresponde a: function xdot = shearbuilding(t,x) global A B S DT N d = t/DT; i = floor(d)+1; if i < N s = S(i)*(i-d)+S(i+1)*(1+d-i); else s = 0; end xdot = A*x + B*s;



b) Integrando directamente con aceleración constante El procedimiento de aceleración constante consiste en expresar la ecuación diferencial de movimiento en pasos asociados al registro. Así, se tiene que la ecuación

 M  qi  C  qi   K  qi  uig  M e Se puede re-escribir como

 K  q   Q  Con

2 4  K    K   C   2  M  t t

 4   Q   ug  M e  2  M  qi    M   2 C   q i  t 

qi 1  qi  q qi 1  qi 1 

2 q  qi t

4 2 q  qi  qi 2 t t

[t_ac,x_ac,xd_ac,xdd_ac] = consacc(K,M,C,e1,S,DT,t2); La función Consacc corresponde a: function [t,q,qd,qdd] = consacc(K,M,C,e1,ug,dt,tt) n = size(K,1); m = length(ug); K = K + 2*C/dt + 4*M/dt^2; R = chol(K); Cm = 2*C + 4*M/dt; Me = M*e1; if tt>m*dt m = tt/dt; ug(m) = 0; end t = zeros(1,m); q = zeros(n,m); qd = zeros(n,m); qdd = zeros(n,m); for i=1:m-1 dQ =-(ug(i+1)-ug(i))*Me + 2*M*qdd(:,i) + Cm*qd(:,i); z = R'\dQ; dq = R\z;



t(i+1) = dt + t(i); q(:,i+1) = dq + q(:,i); qd(:,i+1) = 2*dq/dt - qd(:,i); qdd(:,i+1) = 4*dq/dt^2 -4*qd(:,i)/dt - qdd(:,i); end

c) mediante superposición modal con 2, 3 y 5 modos El procedimiento de superposición modal supone encontrar soluciones ligadas a las formas de vibrar, por lo que se puede encontrar para cada modo ‘i’ lo siguiente:

Vi : vector propio del modo ‘i’

i : frecuencia de vibración del modo ‘i’ i : amortiguamiento del modo ‘i’ M i*  viT  M  vi  1 K i*  viT  K  vi  i 2 Qi*  u g Li  u g  viT  M  e Ci*  viT  M  vi  2  i  i Lo anterior se reemplaza en la ecuación diferencial de movimiento para el modo ‘i’ como:

viT  M  vi  y  viT  C  vi  y  viT  K  vi  y  ug  viT  M  e y  2  i  i  y  i 2  y  ug  Li Para obtener la respuesta completa se requiere la combinación modal sumando todos los modos ‘i’. De esta forma, n

q(t )   y (t ) Vi i 1

for i=1:5 [t_sm,y_sm(i,:),yd_sm(i,:),ydd_sm(i,:)] =consacc(Km(i,i),1,Cm(i,i),fi(:,i)'*M*e1,S,DT,t2); end x_sm_2modo=fi(:,1:2)*y_sm(1:2,:); % con 2 modos x_sm_3modo=fi(:,1:3)*y_sm(1:3,:); % con 3 modos x_sm_5modo=fi(:,1:5)*y_sm(1:5,:); % con 5 modos



Resolución gráfica de respuestas a la excitación de la componente ‘ux’ del registro de Melipilla Para el contraste de las soluciones se determinará el desplazamiento del último piso para el registro de Melipilla, para luego contrastarlos entre sí utilizando colores para diferenciar los desplazamientos obtenidos por los distintos métodos. Los gráficos indican el desplazamiento del último piso (piso 5), en la ordenada del gráfico, en el tiempo, en la abscisa del gráfico. Integrando directamente con ‘ode23’ figure; plot(t_ode,x_ode(:,10), 'g'); title('ode23');

Integrando directamente con aceleración constante figure(2); plot(t_ac,x_ac(5,:), 'c'); title('aceleracion constante');



Mediante superposición modal con 2 modos figure(3); plot(t_sm,x_sm_2modo(5,:), 'k'); title('sp 2 modos');

Mediante superposición modal con 3 modos figure(4); plot(t_sm,x_sm_3modo(5,:), 'b'); title('sp 3 modos');





Mediante superposición modal con 5 modos figure(5); plot(t_sm,x_sm_5modo(5,:), 'r'); title('sp 5 modos');

Contraste entre los métodos En el gráfico que se indica a continuación se bosquejan todas las soluciones enumeradas con anterioridad, superpuestas en un solo gráfico.

Por inspección se aprecia que no hay diferencia notable en la comparación de los métodos, por lo que se puede aseverar que existen diferencias numéricas muy menores entre los métodos.



c) Valores por CQC

Se requiere determinar los desplazamientos modales que son los valores de entrada para realizar el CQC. Los desplazamientos modales se definen como la multiplicación entre la forma modal ‘i’, expresada como vector, por el valor del máximo desplazamiento producto de un registro dado para el período (o frecuencia) asociada a la forma modal, expresada como un escalar. De esta forma se tienen que determinar los valores del espectro de desplazamiento asociados a las frecuencias (períodos) de las formas modales. El espectro de desplazamiento se obtiene a partir de un oscilador de un grado de libertad, del cual se requerirá el máximo desplazamiento para una frecuencia dada del oscilador. Así, se tiene que n

n

 x (a,1)  V (:, i)  S i 1

i

i 1

d

(i) 

L(i) abs( L(i ))

Donde

x(a,1) : desplazamiento modal máximo asociado a los grados de libertad activos ‘a’ del modo ‘i’ V (:,1) : vector de vectores propios Sd (i) : desplazamiento obtenido del espectro de desplazamiento asociado al modo ‘i’ L(i ) : inclusión del signo de las masas modales consideradas en la integración de la respuesta abs( L(i ))

Ahora bien, se hace alusión a los grados activos. Esto se refiere a que se debe condensar estáticamente la matriz para reducirla a los grados de libertad asociados de los desplazamientos horizontales del modelo.


K ab  a   F     Kbb  b   0 

 Kaa   a   Kab   b  F   Kba   a   Kbb   b  0 La rigidez asociada a los términos no activos queda expresada como

 Kb    Kbb   Kba  1

n

n

i 1

i 1

 xi (b,1)   Kbb 

1

 Kba   xi (a,1)

De esta manera se pueden determinar los esfuerzos requeridos al multiplicar los términos de la matriz de rigidez local del elemento por los desplazamientos modales correspondientes. Basta con asociar los desplazamientos de los grados de libertad al vector de desplazamientos modales.



a=[1:5]; b=[6:25]; for i=1:5 x=zeros(25,1); L(i)=V(:,i)'*M*e1; Sd(i)= max(abs(y_sm(i,:)))*L(i)/abs(L(i)); x(a,1)=-V(:,i)*Sd(i); %desplaz. laterales maximos de cada forma modal x(b,1)=inv(Kbb)*Kba*x(a,1); %desplaz. laterales maximos de cada forma modal vi=[1;6;7]; Fzas_izq(:,i)=Kcp1*[zeros(3,1);x(vi,1)]; vd=[1;8;9]; Fzas_der(:,i)=Kcp1*[zeros(3,1);x(vd,1)]; end Notar que hay un signo menos en los desplazamientos laterales máximos de cada forma modal puesto que la definición del elemento columna se definió en sentido contrario al supuesto para el desplazamiento, por lo que los desplazamientos son en sentido contrario y requieren de un signo menos. Con lo anterior se debe aplicar la combinación CQC. Esta se hace para un esfuerzo en particular, es decir, si se requiere conocer el momento en la base de la columna izquierda se obtienen los momentos de en la base de la columna izquierda asociados a cada modo ‘i’ para luego aplicar el operador ij , definido el coeficiente de acoplamiento modal, cuya es expresión corresponde a

ij 



i

8 2 i   j  i 3 j 3

2

  j 2   4 2i j i   j  2

2

El valor requerido se desprende de la ecuación: n

n

r 2   ij  ri  rj i 1 j 1

Donde ri : corresponde al valor del esfuerzo mediante CQC. Notar que se debe aplicar la raíz cuadrada.

ri : corresponde al valor del esfuerzo para el modo i, a combinar con el coeficiente de acoplamiento modal. Mizq_cqc=0; Mder_cqc=0; Nizq_cqc=0; Nder_cqc=0; N_Mizq_cqc=0; N_Mder_cqc=0; for i=1:5 for j=1:5



if i==j rho=1; else rho=(8*bi^2*(w(i)+w(j))*(w(i)^3*w(j)^3)^0.5)/((w(i)^2w(j)^2)^2+4*bi^2*w(i)*w(j)*(w(i)+w(j))^2); end Mizq_cqc=Mizq_cqc+rho*Fzas_izq(3,i)*Fzas_izq(3,j);% Momento en la base columna izquierda primer piso Mder_cqc=Mder_cqc+rho*Fzas_der(3,i)*Fzas_der(3,j);% Momento en la base columna derecha primer piso Nizq_cqc=Nizq_cqc+rho*Fzas_izq(2,i)*Fzas_izq(2,j);% Carga axial en la base columna izquierda primer piso Nder_cqc=Nder_cqc+rho*Fzas_der(2,i)*Fzas_der(2,j);% Carga axial en la base columna derecha primer piso N_Mizq_cqc=N_Mizq_cqc+rho*Fzas_izq(3,i)*Fzas_izq(2,j);% Interaccion N-M en la base columna izquierda primer piso N_Mder_cqc=N_Mder_cqc+rho*Fzas_der(3,i)*Fzas_der(2,j);% Interaccion N-M en la base columna izquierda primer piso end end Mizq_cqc=Mizq_cqc^0.5; Mder_cqc=Mder_cqc^0.5; Nizq_cqc=Nizq_cqc^0.5; Nder_cqc=Nder_cqc^0.5;

De esta manera, los esfuerzos solicitados son Columna izquierda

M  247.3253Tf ·cm N  2.3730 Tf 

N  M  580.8136 Tf 2·cm Columna derecha

M  247.3253Tf ·cm N  2.3730 Tf  N  M  580.8136 Tf 2·cm



d) Superposición modal para la respuesta Mediante superposición modal se puede obtener la respuesta para cada modo. Esto induce a pensar que también se puede obtener la respuesta completa sumando la respuesta de cada modo. El procedimiento de superposición modal supone encontrar soluciones ligadas a las formas de vibrar, por lo que se puede encontrar para cada modo ‘i’ lo siguiente:

Vi : vector propio del modo ‘i’

i : frecuencia de vibración del modo ‘i’ i : amortiguamiento del modo ‘i’ M i*  viT  M  vi  1 K i*  viT  K  vi  i 2 Qi*  u g Li  u g  viT  M  e Ci*  viT  M  vi  2  i  i Lo anterior se reemplaza en la ecuación diferencial de movimiento para el modo ‘i’ como:

viT  M  vi  y  viT  C  vi  y  viT  K  vi  y  ug  viT  M  e

y  2  i  i  y  i 2  y  ug  Li Para obtener la respuesta completa se requiere la combinación modal sumando todos los modos ‘i’. De esta forma, n

x(t )   y (t ) Vi i 1

ti=0; tf=160; % tiempo integración registro for i=1:5 [t_sm,y_sm(i,:),yd_sm(i,:),ydd_sm(i,:)]=consacc(Km(i,i),1,Cm(i,i),1,S,DT,tf); L(i)=V(:,i)'*M*e1; end x_modo1=V(:,1)*y_sm(1,:); % respuesta modo 1 x_modo2=V(:,2)*y_sm(2,:); % respuesta modo 2 x_modo3=V(:,3)*y_sm(3,:); % respuesta modo 3 x_modo4=V(:,4)*y_sm(4,:); % respuesta modo 4 x_modo5=V(:,5)*y_sm(5,:); % respuesta modo 5



x_sm=x_modo1+x_modo2+x_modo3+x_modo4+x_modo5; % superposicion d1_th=0; d1_max=0; d5_th=0; d5_max=0; drift3_th=0; drift3_max=0; Vb_th=0; Vb_max=0; V3_th=0; V3_max=0; a=[1:5]; b=[6:25]; fin_reg=size(x_sm,2); for i=1:fin_reg x_th=zeros(25,1); x_th(a,1)=x_sm(:,i); %desplazamientos laterales maximos de cada forma modal x_th(b,1)=-inv(Kbb)*Kba*x_th(a,1); %desplazamientos pasivos maximos de cada forma modal Ftot=Kt*x_th; d1_th=x_th(1,1); if (abs(d1_max)


Notar que como no se consideraron los grados de libertad en la base de las columnas por saber que son nulos, se tuvo que confeccionar el vector de desplazamientos de las columnas del primer piso para luego multiplicarlo por la matriz de rigidez local de cada columna. Al realizar esto se obtienen las fuerzas y momentos sobre el elemento. Ahora bien, para estimar por CQC el desplazamiento del primer y quinto piso, el drift del piso 3 y los cortes basales y del tercer piso se requiere hacer una estimación mediante una combinación. Para esto se requiere determinar los desplazamientos modales que son los valores de entrada para realizar el CQC. Los desplazamientos modales se definen como la multiplicación entre la forma modal ‘i’, expresada como vector, por el valor del máximo desplazamiento producto de un registro dado para el período (o frecuencia) asociada a la forma modal, expresada como un escalar. De esta forma se tienen que determinar los valores del espectro de desplazamiento asociados a las frecuencias (períodos) de las formas modales. El espectro de desplazamiento se obtiene a partir de un oscilador de un grado de libertad, del cual se requerirá el máximo desplazamiento para una frecuencia dada del oscilador. Así, se tiene que n

n

i 1

i 1

 xi (a,1)  V (:, i)  Sd (i)  L(i) Donde

x(a,1) : desplazamiento modal máximo asociado a los grados de libertad activos ‘a’ del modo ‘i’ V (:,1) : vector de vectores propios Sd (i) : desplazamiento obtenido del espectro de desplazamiento asociado al modo ‘i’ L(i) : la masa modal considerada en la integración de la respuesta

Ahora bien, se hace alusión a los grados activos. Esto se refiere a que se debe condensar estáticamente la matriz para reducirla a los grados de libertad asociados de los desplazamientos horizontales del modelo.


K ab  a   F     Kbb  b   0 

 Kaa   a   Kab   b  F   Kba   a   Kbb   b  0 La rigidez asociada a los términos no activos queda expresada como

 Kb    Kbb   Kba  1

n

n

 x (b,1)   K   K   x (a,1) i 1

i

i 1

1

bb

ba

i



De esta manera se pueden determinar los esfuerzos requeridos al multiplicar los términos de la matriz de rigidez local del elemento por los desplazamientos modales correspondientes. Basta con asociar los desplazamientos de los grados de libertad al vector de desplazamientos modales. a=[1:5]; b=[6:25]; for i=1:5 x_cqc=zeros(25,1); Sd(i)= max(abs(y_sm(i,:))); x_cqc(a,1)=-V(:,i)*Sd(i)*L(i); %desp laterales maximos de cada forma modal x_cqc(b,1)=inv(Kbb)*Kba*x_cqc(a,1); %desp pasivos maximos de cada forma modal d1(i)=x_cqc(1,1); d5(i)=x_cqc(5,1); drift3(i)=x_cqc(3,1)-x_cqc(2,1); Vb(i)=Sd(i)*(L(i)*w(i))^2; Fcp1=Kcp2*[x_cqc(2,1) x_cqc(10,1) x_cqc(11,1) x_cqc(3,1) x_cqc(14,1) x_cqc(15,1)]'; Fcp2=Kcp2*[x_cqc(2,1) x_cqc(12,1) x_cqc(13,1) x_cqc(3,1) x_cqc(16,1) x_cqc(17,1)]'; V3(i)=Fcp1(1,1)+Fcp2(1,1); end d1_cqc=0; d5_cqc=0; drift3_cqc=0; Vb_cqc=0; V3_cqc=0; for i=1:5 for j=1:5 if i==j rho=1; else rho=(8*bi^2*(w(i)+w(j))*(w(i)^3*w(j)^3)^0.5)/((w(i)^2w(j)^2)^2+4*bi^2*w(i)*w(j)*(w(i)+w(j))^2); end d1_cqc=d1_cqc+rho*d1(i)*d1(j);% desplazamiento primer piso d5_cqc=d5_cqc+rho*d5(i)*d5(j);% desplazamiento quinto piso drift3_cqc=drift3_cqc+rho*drift3(i)*drift3(j);% drift tercer piso Vb_cqc=Vb_cqc+rho*Vb(i)*Vb(j);% Carga axial en la base columna derecha primer piso V3_cqc=V3_cqc+rho*V3(i)*V3(j);% Interaccion N-M en la base columna izquierda primer piso end end d1_cqc=d1_cqc^0.5; d5_cqc=d5_cqc^0.5; drift3_cqc=drift3_cqc^0.5;



Vb_cqc=Vb_cqc^0.5; V3_cqc=V3_cqc^0.5;

Resultados Para comparar los resultados, éstos se presentan como pares ordenados Desplazamiento primer piso desplaz_primer_piso=[d1_cqc,d1_max] 0.5304 0.5532 Desplazamiento quinto piso desplaz_quinto_piso=[d5_cqc,d5_max] 2.2773 2.2636 Desplazamiento relativo del tercer piso drift_tercer_piso=[drift3_cqc,drift3_max] 0.5412 0.5409 Corte basal corte_basal=[Vb_cqc,Vb_max] 1.8387 1.9415 Corte tercer piso corte_tercer_piso=[V3_cqc,V3_max] 1.3868 1.3893

Para comparar la calidad de los resultados obtenidos se presentan los errores para cada par ordenado Error desplazamiento primer piso error_desplaz_primer_piso=(d1_max-d1_cqc)/d1_max*100 4.1262 Error desplazamiento quinto piso error_desplaz_quinto_piso=(d5_max-d5_cqc)/d5_max*100 -0.6074 Error desplazamiento relativo del tercer piso error_drift_tercer_piso=(drift3_max-drift3_cqc)/drift3_max*100 -0.0660



Error corte basal error_corte_basal=(Vb_max-Vb_cqc)/Vb_max*100 5.2934 Error corte tercer piso error_corte_tercer_piso=(V3_max-V3_cqc)/V3_max*100 0.1842


Matlab Para Dinamica de Estructuras RevB

Recommend Documents