Input file for Response Spectra Analysis of 4 Story Shear Building


/*
 *  ================================================================
 *  Response Spectra Analysis of 4 Story Shear Building
 *
 *  Written By : Mark Austin                               June 1996
 *  ================================================================
 */

/* 
 *  [a] : Setup Matrix for Piece-wise linear approximation to earthquake
 *        acceleration Spectra
 *
 *        Column 1 : Natural Period (sec)
 *        Column 2 : Spectral Acceleration (at 2% damping)
 */ 

   print "*** ACCELERATION SPECTRA FOR MODAL ANALYSIS \n";
   print "*** ======================================= \n";

   npoints = 18;
   spectra = Matrix( [ npoints , 2] );
   spectra = ColumnUnits( spectra, [sec],      [1]);
   spectra = ColumnUnits( spectra, [cm/sec^2], [2]);

   spectra [ 1][1] = 0.0 sec;  spectra [ 1][2] = 981.0*0.15 cm/sec/sec;
   spectra [ 2][1] = 0.1 sec;  spectra [ 2][2] = 981.0*0.18 cm/sec/sec;
   spectra [ 3][1] = 0.2 sec;  spectra [ 3][2] = 981.0*0.25 cm/sec/sec;
   spectra [ 4][1] = 0.3 sec;  spectra [ 4][2] = 981.0*0.38 cm/sec/sec;
   spectra [ 5][1] = 0.4 sec;  spectra [ 5][2] = 981.0*0.50 cm/sec/sec;
   spectra [ 6][1] = 0.5 sec;  spectra [ 6][2] = 981.0*0.50 cm/sec/sec;
   spectra [ 7][1] = 0.6 sec;  spectra [ 7][2] = 981.0*0.40 cm/sec/sec;
   spectra [ 8][1] = 0.8 sec;  spectra [ 8][2] = 981.0*0.32 cm/sec/sec;
   spectra [ 9][1] = 1.0 sec;  spectra [ 9][2] = 981.0*0.25 cm/sec/sec;
   spectra [10][1] = 1.2 sec;  spectra [10][2] = 981.0*0.19 cm/sec/sec;
   spectra [11][1] = 1.4 sec;  spectra [11][2] = 981.0*0.14 cm/sec/sec;
   spectra [12][1] = 1.6 sec;  spectra [12][2] = 981.0*0.10 cm/sec/sec;
   spectra [13][1] = 1.8 sec;  spectra [13][2] = 981.0*0.07 cm/sec/sec;
   spectra [14][1] = 2.0 sec;  spectra [14][2] = 981.0*0.05 cm/sec/sec;
   spectra [15][1] = 2.4 sec;  spectra [15][2] = 981.0*0.04 cm/sec/sec;
   spectra [16][1] = 2.8 sec;  spectra [16][2] = 981.0*0.03 cm/sec/sec;
   spectra [17][1] = 3.2 sec;  spectra [17][2] = 981.0*0.02 cm/sec/sec;
   spectra [18][1] = 3.4 sec;  spectra [18][2] = 981.0*0.02 cm/sec/sec;

   PrintMatrix( spectra );

/* [b] : Setup mass and stiffness matrices */ 

   mass = ColumnUnits( 1500*[ 1, 0, 0, 0;
                              0, 2, 0, 0;
                              0, 0, 2, 0;
                              0, 0, 0, 3], [kg] );

   stiff = ColumnUnits( 800*[ 1, -1,  0,  0;
                             -1,  3, -2,  0;
                              0, -2,  5, -3;
                              0,  0, -3,  7], [kN/m] );

   PrintMatrix(mass, stiff);

/* [c] : Calculate natural periods of vibration and mode shapes */

   no_eigen = 3;
   eigen       = Eigen( stiff, mass, [ no_eigen ] );
   eigenvalue  = Eigenvalue( eigen );
   eigenvector = Eigenvector ( eigen );

   period = ColumnUnits( Matrix( [ no_eigen,1 ] ), [ sec ] );
   for( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
        period [ii][1] = 2*PI/sqrt(eigenvalue[ii][1]);
   }

   print "\n";
   for( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
        print "Mode", ii, " : w^2 = ", eigenvalue[ii][1];
        print " T = ", period [ii][1], "\n";
   }

   PrintMatrix( period );
   PrintMatrix( eigenvector );

/* [d] : Find Spectral Accelerations at Modal Periods */

   SpectralAccn = ColumnUnits( Matrix( [ no_eigen,1 ] ), [ m/sec^2 ] );
   for( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
        for( ij = 1; ij < npoints; ij = ij + 1 ) {
             period1 = spectra [  ij][1];
             period2 = spectra [ij+1][1];
             if(period [ii][1] >= period1 && period [ii][1] < period2 ) {
                dAccn   = spectra [ij+1][2] - spectra [ij][2];
                dPeriod = (period [ii][1] - period1)/(period2 - period1);
                SpectralAccn[ii][1] = spectra [ij][2] + dPeriod*dAccn;
             }
        }
   }

   PrintMatrix( SpectralAccn );

/* [e] : Generalised mass, stiffness, and loading matrices */

   eigenTrans = Trans (eigenvector);
   gmass      = eigenTrans*mass*eigenvector;
   gstiff     = eigenTrans*stiff*eigenvector;

   gload = eigenTrans*mass*[ 1; 1; 1; 1 ];

   PrintMatrix( gmass, gstiff, gload );

/* [f] : Compute and print floor level displacements */

   Y = ColumnUnits( Matrix([no_eigen,no_eigen]), [ m ] );
   for( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
        Y [ii][ii] = (gload[ii][1]/gmass[ii][ii]) *
                     (SpectralAccn[ii][1]/eigenvalue[ii][1]);
   }

   modaldispl = ColumnUnits ( eigenvector*Y, [cm]);
   PrintMatrix( Y , modaldispl );

   print "\n";
   print "Maximum Likely Floor Displacements (using SRSS)     \n";
   print "===============================================   \n\n";
   print "     Floor          Mode         Modes         Modes\n";
   print "        No             1         1 & 2    1, 2 and 3\n";

   for( ii = 1; ii <= 4; ii = ii + 1 ) {
        print (5-ii);
        for( ij = 1; ij <= no_eigen; ij = ij + 1 ) {
             sum = 0.0 m^2;
             for( ik = 1; ik <= ij; ik = ik + 1 ) {
                  sum = sum + modaldispl [ii][ik] * modaldispl [ii][ik];
             }
             print sqrt(sum) (cm); 
        }
        print "\n";
   }
 
   print "\n";
   print "Maximum Possible Floor Displacements (absolute values) \n";
   print "====================================================== \n\n";
   print "     Floor          Mode         Modes         Modes\n";
   print "        No             1         1 & 2    1, 2 and 3\n";

   for( ii = 1; ii <= 4; ii = ii + 1 ) {
        print (5-ii);
        for( ij = 1; ij <= no_eigen; ij = ij + 1 ) {
             sum = 0.0 cm;
             for( ik = 1; ik <= ij; ik = ik + 1 ) {
                  sum = sum + abs( modaldispl [ii][ik] );
             }
             print sum (cm); 
        }
        print "\n";
   }

/* [g] : Compute and print story drifts */

   print "\n";
   print "Maximum Likely Story Drifts (using SRSS) \n";
   print "======================================== \n\n";
   print "     Story          Mode         Modes         Modes\n";
   print "        No             1         1 & 2    1, 2 and 3\n";

   for( ii = 1; ii <= 4; ii = ii + 1 ) {
        print (5-ii);
        for( ij = 1; ij <= no_eigen; ij = ij + 1 ) {
             sum = 0.0 cm^2;
             for( ik = 1; ik <= ij; ik = ik + 1 ) {
                  if (ii == 4) then {
                      sum = sum + modaldispl[4][ik]^2;
                  } else {
                      sum = sum + (modaldispl [ii][ik] - modaldispl [ii+1][ik])^2;
                  } 
             }
             print sqrt(sum) (cm); 
        }
        print "\n";
   }

   print "\n";
   print "Maximum Possible Story Drifts (absolute values) \n";
   print "=============================================== \n\n";
   print "     Story          Mode         Modes         Modes\n";
   print "        No             1         1 & 2    1, 2 and 3\n";

   for( ii = 1; ii <= 4; ii = ii + 1 ) {
        print (5-ii);
        for( ij = 1; ij <= no_eigen; ij = ij + 1 ) {
             sum = 0.0 cm;
             for( ik = 1; ik <= ij; ik = ik + 1 ) {
                  if (ii == 4) then {
                      sum = sum + abs (modaldispl[4][ik]);
                  } else {
                      sum = sum + abs((modaldispl [ii][ik] - modaldispl [ii+1][ik]));
                  } 
             }
             print sum (cm); 
        }
        print "\n";
   }

/* [h] : Compute and print equivalent d.o.f. forces in each mode */

   print "\n";
   print "Inertia Forces for each mode \n";
   print "============================ \n\n";

   inertia_forces = stiff*modaldispl;
   PrintMatrix( inertia_forces );

/* [i] : Compute and print base shear force */

   print "\n";
   print "Shear Forces (at base of the structure) \n";
   print "======================================= \n\n";

   base_shear_forces = [1,1,1,1] * inertia_forces;
   PrintMatrix( base_shear_forces );

   shear1 = 0.0 N^2; shear2 = 0.0 N;
   for ( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
         shear1 = shear1 + base_shear_forces[ 1 ][ ii ]^2;
         shear2 = shear2 + abs( base_shear_forces[ 1 ][ ii ] );
   }

   print "\n";
   print "Base Shear Force : Maximum Likely   = ", sqrt(shear1) ,"\n";
   print "                 : Maximum Possible = ", shear2       ,"\n";

/* [j] : Compute and print overturning moments */

   print "\n";
   print "Overturning Moments (at base of the structure) \n";
   print "============================================== \n\n";

   floor_heights            = [12 m, 9 m, 6 m, 3 m];
   base_overturning_moments = floor_heights * inertia_forces;
   PrintMatrix( base_overturning_moments );

   mom1 = 0.0 (N*m)^2; mom2 = 0.0 N*m;
   for ( ii = 1; ii <= no_eigen; ii = ii + 1 ) {
         mom1 = mom1 + base_overturning_moments[ 1 ][ ii ]^2;
         mom2 = mom2 + abs( base_overturning_moments[ 1 ][ ii ] );
   }

   print "\n";
   print "Overturning Moments : Maximum Likely   = ", sqrt(mom1) ,"\n";
   print "                    : Maximum Possible = ", mom2       ,"\n";

   quit;


Developed in June 1996 by Mark Austin
Last Modified June 27, 1996
Copyright © 1996, Mark Austin, Department of Civil Engineering, University of Maryland