Analisis modal espectral

Realizar el análisis dinámico tridimensional de la estructura siguiente que está configurada con pórticos no arriostrados de concreto reforzado y con placas aligeradas. Cuenta con un modulo de elasticidad, Ec igual a 20000MPa, y tiene las siguientes propiedades geométricas y de masa:
Dimensión de vigas: 0.40x0.45m
Dimensión de columnas: 0.5x0.45m
Masa en los pisos: 1000 kg/m2
Masa en lacubierta: 700 kg/m2

Vamos a empezar hallando la matriz de masas de la estructura teniendo en cuenta las masas rotacionales:
mr=mA*Jo
mt=m*A
A=5*4+12*5=80m2
Jo=1288.3334m4

El momento polar de inercia de la figura fue calculado en Autocad
Omitiendo más cálculos esta es la matriz construida:
Matriz de Masas 9x9 | | | | | | | |
56000,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00| 0,00 |
0,00 | 56000,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 11272,92 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 80000,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 80000,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 16104,17 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 |80000,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 80000,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 16104,17 |

Otro cálculo es el centroide de la figura donde se van a condensar los grados de libertad:
Xc=6m
Zc= -1.25m
Para el análisis tridimensional dividimos la estructura en pórticos planos y le calculamos la rigidez G.

División de laestructura en pórticos planos:
Pórticos A, B

Código de la programación en Matlab para la rigidez de un elemento y proceso hasta ensamblar la matriz de estos tipos de pórticos.
Parte del código en Matlab
%rigidez de vigas de 4m de longitud::::::::::::
E=2*(10^10);
b = 0.40;
h = 0.45;
I=b*h^3/12;
A=b*h;
L1 = 4;
L=L1;
ThetaR1=0;
Theta1=ThetaR1*pi/180;

K1=[A*E/L 0 0 -A*E/L 0 0;0 12*E*I/(L^3) 6*E*I/(L^2) 0 -12*E*I/(L^3) 6*E*I/(L^2);
0 6*E*I/(L^2) 4*E*I/(L) 0 -6*E*I/(L^2) 2*E*I/L;
-A*E/L 0 0 A*E/L 0 0;
0 -12*E*I/(L^3) -6*E*I/(L^2) 0 12*E*I/(L^3) -6*E*I/(L^2);
0 6*E*I/(L^2) 2*E*I/L 0 -6*E*I/(L^2) 4*E*I/L];

T1=[cos(Theta1) sin(Theta1) 0 0 0 0;
-sin(Theta1) cos(Theta1) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 cos(Theta1) sin(Theta1) 0;
0 0 0-sin(Theta1) cos(Theta1) 0;
0 0 0 0 0 1];

K1 MPa-m | | | | | |
9,00E+15 | 0,00 | 0,00 | -9,00E+15 | 0,00 | 0,00 |
0,00 | 1,14E+16 | 2,28E+16 | 0,00 | -1,14E+16 | 2,28E+16 |
0,00 | 2,28E+16 | 6,08E+15 | 0,00 | -2,28E+16 | 3,04E+16 |
-9,00E+15 | 0,00 | 0,00 | 9,00E+15 | 0,00 | 0,00 |
0,00 | -1,14E+16 | -2,28E+16 | 0,00 | 1,14E+16 | -2,28E+16 |
0,00 | 2,28E+16 | 3,04E+16 |0,00 | -2,28E+16 | 6,08E+15 |

Matriz de transformación T1
T1 | | | | | |
1,00 | 0,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 1,00 | 0,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 1,00 | 0,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 1,00 | 0,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 1,00 | 0,00 |
0,00 | 0,00 | 0,00 | 0,00 | 0,00 | 1,00 |

Matriz global del elemento 1
K1Ga= T1'*K1*T1;*código matlab

K1Ga | | | | | |
9,00E+15 | 0,00 | 0,00 | -9,00E+15 | 0,00 | 0,00 |
0,00 | 1,14E+16 | 2,28E+16 | 0,00 | -1,14E+16 | 2,28E+16 |
0,00 | 2,28E+16 | 6,08E+15 | 0,00 | -2,28E+16 | 3,04E+16 |
-9,00E+15 | 0,00 | 0,00 | 9,00E+15 | 0,00 | 0,00 |
0,00 | -1,14E+16 | -2,28E+16 | 0,00 | 1,14E+16 | -2,28E+16 |
0,00 | 2,28E+16 | 3,04E+16 | 0,00 | -2,28E+16 | 6,08E+15 |

Para lossiguientes miembros tipo vigas horizontales del pórtico tenemos que:
K2Ga=K1Ga K3Ga=K1Ga K4Ga=K1Ga
K5Ga=K1Ga K6Ga=K1Ga K7Ga=K1Ga
K8Ga=K1Ga K9Ga=K1Ga
Para los elementos verticales tipo columnas y que también componen el pórtico tenemos una rigidez calculada de la siguiente forma:
Parte del código en Matlab
% rigidez de columnas de 3m de longitud::::::::::::

E=2*(10^10);...