Input file for Analysis of Highway Bridge Structure.

You will need Aladdin 2.0 to run this problem.


/*
 *  =====================================================================
 *  Analysis of 4 span isolated bridge, 3 piers, 2 abutments, 5 isolators
 *
 *  Modeling Assumptions:
 *
 *     --  Bottoms of the abuntments' isolators and piers are fixed.
 *     --  Deck is modeled as FRAME_3D elements.                    
 *     --  Pier columns are modeled as FRAME_3D elements.          
 *     --  Isolators are modeled as FIBER_3D elements.             
 *     --  Perfect elastic-plastic for isolators.                   
 *     --  Peak earthquake = 0.5*g. (rec.6)                         
 *     --  5% of structural damping is added, C=A0*M+A1+K.          
 *
 *  Written By : Wane-Jang Lin                               October 1997
 *  =====================================================================
 */

NDimension         = 3;
NDofPerNode        = 6;
MaxNodesPerElement = 2;
GaussIntegPts      = 2;

StartMesh();

H2    = 14 m;
H3    =  7 m;
H4    = 21 m;
iso_h = 20 cm;

/* Assemble mesh of nodes and elements */

   total_node = 0;
   x = 0 m;  y = 0 m;  z = 0 m;
   for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) {
       AddNode( node_no, [x,y,z] );
       x = x + 50 m;
   }

   total_node = node_no-1;
   x = 0 m;  y = -iso_h;  z = 0 m;
   for( node_no=total_node+1 ; x<=200m ; node_no=node_no+1 ) {
       AddNode( node_no, [x,y,z] );
       x = x + 50 m;
   }

   total_node = node_no-1;
   x = 50 m;  y = -iso_h-3.5m;  z = 0 m;
   for( node_no=total_node+1 ; y>=-(iso_h+H2) ; node_no=node_no+1 ) {
       AddNode( node_no, [x,y,z] );
       y = y - 3.5 m;
   }

   total_node = node_no-1;
   fix2       = total_node;
   x = 100 m;  y = -iso_h-3.5m;  z = 0 m;
   for( node_no=total_node+1 ; y>=-(iso_h+H3) ; node_no=node_no+1 ) {
       AddNode( node_no, [x,y,z] );
       y = y - 3.5 m;
   }

   total_node = node_no-1;
   fix3       = total_node;
   x = 150 m;  y = -iso_h-3.5m;  z = 0 m;
   for( node_no=total_node+1 ; y>=-(iso_h+H4) ; node_no=node_no+1 ) {
       AddNode( node_no, [x,y,z] );
       y = y - 3.5 m;
   }

   total_node = node_no-1;
   fix4       = total_node;

   total_elmt = 0;
   for( elmt_no=total_elmt+1 ; elmt_no<=4 ; elmt_no=elmt_no+1 ) {
       AddElmt( elmt_no, [elmt_no, elmt_no+1], "deck_attr" );
   }

   AddElmt(  5, [1,6],  "isolator_attr1" );
   AddElmt(  6, [2,7],  "isolator_attr2" );
   AddElmt(  7, [3,8],  "isolator_attr3" );
   AddElmt(  8, [4,9],  "isolator_attr4" );
   AddElmt(  9, [5,10], "isolator_attr1" );

   total_elmt = 9;
   node_no = fix2 - H2/(3.5m) + 1;
   AddElmt( total_elmt+1, [7, node_no], "column_attr2" );

   total_elmt = total_elmt+1;
   for( elmt_no=total_elmt+1 ; node_no < fix2 ; elmt_no=elmt_no+1 ) {
       AddElmt( elmt_no, [node_no, node_no+1], "column_attr2" );
       node_no = node_no+1;
   }

   total_elmt = elmt_no-1;
   bot_elmt2  = total_elmt;
   node_no    = fix3 - H3/(3.5m) + 1;
   AddElmt( total_elmt+1, [8, node_no], "column_attr3" );
   total_elmt = total_elmt+1;
   for( elmt_no=total_elmt+1 ; node_no < fix3 ; elmt_no=elmt_no+1 ) {
       AddElmt( elmt_no, [node_no, node_no+1], "column_attr3" );
       node_no = node_no+1;
   }

   total_elmt = elmt_no-1;
   bot_elmt3  = total_elmt;
   node_no    = fix4 - H4/(3.5m) + 1;
   AddElmt( total_elmt+1, [9, node_no], "column_attr4" );
   total_elmt = total_elmt+1;
   for( elmt_no=total_elmt+1 ; node_no < fix4 ; elmt_no=elmt_no+1 ) {
       AddElmt( elmt_no, [node_no, node_no+1], "column_attr4" );
       node_no = node_no+1;
   }

   total_elmt = elmt_no-1;
   bot_elmt4  = total_elmt;

   ElementAttr("deck_attr") { type     = "FRAME_3D";
                              section  = "deck_sec";
                              material = "deck_mat";
                            }
   SectionAttr("deck_sec")  { Iyy = 103.77 m^4; Izz = 5.26 m^4;
                              J       = 16.13 m^4;
                              area    = 6.88 m^2;
                              unit_weight = 200 kN/m;
                            }
   MaterialAttr("deck_mat") { poisson = 0.25;
                              E       = 2.3e7 kN/m^2;
                            }

   b  = 2.2    m;  h = 4.0    m;
   Ix = 7.39 m^4; Iy = 0.67 m^4;
   K2 =       30700 kN/m; K3 =      124400 kN/m; K4 = 13650 kN/m;
   E2 = K2*H2*H2*H2/3/Ix; E3 = K3*H3*H3*H3/3/Ix; E4 = K4*H4*H4*H4/3/Ix;

   ElementAttr("column_attr2") { type     = "FRAME_3D";
                                 section  = "col_sec";
                                 material = "col_mat2";
                               }

   ElementAttr("column_attr3") { type     = "FRAME_3D";
                                 section  = "col_sec";
                                 material = "col_mat3";
                               }

   ElementAttr("column_attr4") { type     = "FRAME_3D";
                                 section  = "col_sec";
                                 material = "col_mat4";
                               }

   SectionAttr("col_sec") { area  = 4.16 m^2;
                            width = h;  depth = b;
                            Iyy   = Ix; Izz   = Iy;
                            unit_weight  = 104 kN/m;
                            shear_factor = 1.0;
                          }

   MaterialAttr("col_mat2") { poisson = 0.25;  E = E2; }
   MaterialAttr("col_mat3") { poisson = 0.25;  E = E3; }
   MaterialAttr("col_mat4") { poisson = 0.25;  E = E4; }

   h = 50 cm;
   b = 50 cm;
   E_lead = 356200 kN/m^2;
   G1  = 4240 kN/m^2;  G2 = 12880 kN/m^2;   G3 = 9200 kN/m^2;   G4 = 37600 kN/m^2;
   fvy = 8400 kN/m^2;

   ElementAttr("isolator_attr1") { type     = "FIBER_3D";
                                   section  = "iso_sec";
                                   material = "iso_mat1";
                                   fiber    = "iso_fib";
                                 }

   ElementAttr("isolator_attr2") { type     = "FIBER_3D";
                                   section  = "iso_sec";
                                   material = "iso_mat2";
                                   fiber    = "iso_fib";
                                 }

   ElementAttr("isolator_attr3") { type     = "FIBER_3D";
                                   section  = "iso_sec";
                                   material = "iso_mat3";
                                   fiber    = "iso_fib";
                                 }

   ElementAttr("isolator_attr4") { type     = "FIBER_3D";
                                   section  = "iso_sec";
                                   material = "iso_mat4";
                                   fiber    = "iso_fib";
                                 }

   SectionAttr("iso_sec") { area  = b*h;
                            width = b;
                            depth = h;
                            unit_weight  = 0.0001 kN/m;
                            shear_factor = 1.0;
                          }

   MaterialAttr("iso_mat1") { poisson = 0.25; G = G1; Gt = 0.001*G1; shear_yield = fvy/2; }
   MaterialAttr("iso_mat2") { poisson = 0.25; G = G2; Gt = 0.001*G2; shear_yield = fvy; }
   MaterialAttr("iso_mat3") { poisson = 0.25; G = G3; Gt = 0.001*G3; shear_yield = fvy; }
   MaterialAttr("iso_mat4") { poisson = 0.25; G = G4; Gt = 0.001*G4; shear_yield = fvy; }

   no_fiber = 4;
   iso_fcoord = Matrix([ 2, no_fiber ]);
   iso_farea  = Matrix([ 1, no_fiber ]);
   iso_fmap   = Matrix([ 1, no_fiber ]);

   iso_fcoord[1][1] =  b/4;
   iso_fcoord[2][1] =  h/4;
   iso_fcoord[1][2] = -b/4;
   iso_fcoord[2][2] =  h/4;
   iso_fcoord[1][3] =  b/4;
   iso_fcoord[2][3] = -h/4;
   iso_fcoord[1][4] = -b/4;
   iso_fcoord[2][4] = -h/4;

   for( i=1 ; i <= no_fiber ; i=i+1 ) {
      iso_farea[1][i]  = b*h/no_fiber;
      iso_fmap[1][i]   = 1;
   }

   iso_fattr = Matrix([ 3, 1 ]);
   iso_fattr[1][1] = E_lead;
   iso_fattr[2][1] = E_lead;
   iso_fattr[3][1] = fvy*1e6;

   FiberAttr( no_fiber, "iso_fib" ) { FiberMaterialAttr = iso_fattr;
                                      FiberCoordinate   = iso_fcoord;
                                      FiberArea         = iso_farea;
                                      FiberMaterialMap  = iso_fmap;
                                    }

/* Apply boundary conditions */

   FixNode(   1, [1,0,0,1,0,0]);
   FixNode(   5, [1,0,0,1,0,0]);
   FixNode(   6, [1,1,1,1,1,1]);
   FixNode(  10, [1,1,1,1,1,1]);
   FixNode(fix2, [1,1,1,1,1,1]);
   FixNode(fix3, [1,1,1,1,1,1]);
   FixNode(fix4, [1,1,1,1,1,1]);

/* Compile and print f.e. mesh    */

   EndMesh();
   PrintMesh();

/* Compute initial stiffness and mass matrices */

   stiff = Stiff();
   mass  = Mass([1]);
   mass_inv = Inverse( mass );

/* Compute and print eigenvalues */

   no_eigen   = 2;
   eigen      = Eigen(stiff, mass, [no_eigen]);
   eigenvalue = Eigenvalue(eigen);
   for(i=1 ; i <= no_eigen ; i=i+1) {
      w = sqrt( eigenvalue[i][1] );
      print "Mode",i,",\tw=",w,",\tperiod=",2*PI/w,"\n";
   }

   eigenvector = Eigenvector(eigen);
   PrintMatrix(eigenvector);

   wa = sqrt( eigenvalue[1][1] );
   wb = sqrt( eigenvalue[2][1] );

/* Setup Rayleigh damping and damping matrix */

   rdamping = 0.20;
   A0   = 2*rdamping*wa*wb/(wa+wb);
   A1   = 2*rdamping/(wa+wb);
   damp = A0*mass + A1*stiff;

/* El Centro Earthquake record */

   Elcentro1 = ColumnUnits( [
      4.56;    13.77;     6.13;     3.73;     1.32;    -6.81;   -16.22;   -22.41;
    -26.67;   -24.10;   -21.86;   -14.02;    -1.62;     1.67;    18.12;    39.26;
     28.40;    16.72;    16.06;    12.54;     0.15;    -8.60;   -14.00;   -12.01;
     -3.34;     5.14;    14.04;    16.36;    29.90;    43.09;    32.13;    17.28;
     18.44;    22.31;    13.24;    -0.91;    -6.81;   -15.20;   -14.88;    -9.54;
    -13.32;   -10.60;    11.80;    10.20;   -12.90;   -22.88;     1.42;    21.73;
      7.92;     6.47;    11.39;    -8.29;   -33.44;   -26.25;   -23.99;   -31.05;
    -25.42;   -26.97;   -38.48;   -35.01;   -27.78;   -33.35;   -29.56;   -13.60;
    -13.97;   -11.63;    -1.31;    13.86;    35.71;    27.62;    25.33;    46.45;
     44.12;    28.02;    19.13;    12.57;    11.60;    12.29;     9.73;    10.38;
     18.93;    29.43;    30.25;    26.23;    23.15;    22.29;    21.11;    14.84;
     11.02;    14.74;    21.65;    26.66;    28.45;    30.34;    30.21;    28.82;
     31.31;    26.66;    15.29;     9.10;     3.58;    -6.65;   -15.33;   -18.05;
    -17.46;   -11.68;    -3.48;    -5.25;   -15.43;   -25.38;   -32.32;   -29.67;
    -22.01;   -17.97;   -22.34;   -30.51;   -38.89;   -45.29;   -55.26;   -67.84;
    -76.74;   -84.58;   -86.83;   -86.63;   -80.87;   -69.83;   -59.81;   -54.51;
    -48.94;   -39.94;   -38.67;   -40.13;   -38.04;   -34.78;   -31.30;   -21.29;
     -8.85;     3.58;    13.97;    30.88;    35.09;     8.84;   -10.41;   -16.30;
    -20.10;    -3.69;    33.71;    51.69;    43.52;    27.04;    34.70;    38.30;
     36.60;    29.85;    -9.22;   -23.98;   -17.28;   -14.20;    11.58;    47.99;
     64.26;    70.36;    68.93;    62.67;    45.30;    28.43;    23.73;    28.51;
     34.46;    43.12;    43.39;    26.64;    21.23;    31.23;    40.30;    16.76;
    -10.80;   -22.88;    -0.51;    12.51;   -10.01;   -12.52;    -6.05;   -24.18;
    -43.17;     7.01;    50.06;    54.19;    68.83;    69.26;    57.53;     8.73;
    -38.98;   -43.59;   -51.12;   -41.86;   -19.51;   -15.78;   -11.54;    -2.18;
     -2.19;    20.84;    23.40;    -0.91;   -11.45;   -22.94;   -40.84;   -53.17;
    -44.02;   -21.96;    -3.06;    -0.49;    -9.59;   -11.83;    -9.86;    -7.02;
      0.88;    24.00;    50.57;    56.01;    51.75;    52.30;    47.79;    16.77;
     -6.09;   -20.13;   -20.42;    -4.35;    -0.33;    -6.33;     1.74;    19.90;
     37.97;    53.36;    57.68;    53.32;    28.37;     7.54;     6.90;     4.73;
     -4.97;   -21.97;   -28.16;   -23.25;   -28.43;   -18.57;    -3.93;   -16.60;
    -31.92;   -35.08;   -34.77;   -30.85;   -30.37;    -6.83;    24.57;    27.39;
     23.72;    29.30;     5.79;   -35.35;   -54.31;   -57.00;   -59.38;   -19.00;
     32.14;     7.01;   -21.44;   -21.26;   -26.01;   -35.66;   -33.68;   -24.09;
    -21.39;   -26.29;   -27.73;   -23.91;   -14.27;    -6.41;    -0.59;     4.21;
      1.06;    -4.59;   -15.16;   -18.12;    -7.09;     2.41;     6.08;     1.20;
    -10.53;    -8.69;    -4.70;   -13.05;    -9.51;    -0.67;    -7.10;    -7.19;
      5.36;     9.80;     7.87;    12.63;    20.28;    22.44;    24.86;    24.31;
     15.25;     6.54;     6.86;    11.08;    18.12;    23.64;    22.47;    23.76;
     21.57;    19.80;    25.47;    27.94;    14.81;    -6.94;   -13.57;    -7.80;
     -2.65;    -3.97;    -1.01;     0.20;     5.14;     2.31;   -16.89;   -18.18;
    -11.55;   -12.10;    -3.07;     6.16;     8.87;    19.72;    21.31;    16.45;
     12.70;     4.15;    -1.59;     0.57;    11.04;    21.39;    27.42;    15.82;
     -9.25;   -23.67;   -17.34;     0.67;     8.69;    11.01;    12.01;    11.12;
     13.62;    11.31;     9.71;    17.68;    15.83;     5.59;     0.60;    -2.53;
     -4.46;    -4.09;    -7.81;   -15.06;   -19.43;   -15.49;    -8.93;    -5.44;
     -1.33;    -1.35;     7.19;    12.52;     0.21;    -4.89;     0.22;    11.51;
     20.17;    18.62;    16.72;    13.74;    10.62;    11.14;    12.74;    13.34;
     14.16;    14.59;    13.81;    10.94;     5.91;     5.74;    14.53;    19.56;
      9.27;    -6.75;   -14.69;   -14.25;   -11.48;    -3.96;     6.54;    10.59;
      2.26;   -12.97;   -24.01;   -25.73;   -18.75;    -6.99;    -2.22;     0.17;
      3.63;    -1.70;    -4.73;     6.13;    20.26;    22.36;    19.58;     9.14;
     -0.09;    -0.23;    -1.00;     4.82;    15.39;    19.75;     7.45;    -0.99;
      0.77;    -1.64;    -9.34;   -14.34;   -13.86;    -6.61;     5.49;     6.30;
      2.08;    -2.88;    -2.25;     2.99;     3.30;    -1.57;    -6.86;    -8.44;
    -11.03;   -16.68;   -17.43;   -12.39;   -12.20;   -15.43;   -16.28;   -17.99;
    -20.21;   -21.58;   -23.05;   -23.79;   -23.64;   -21.60;   -15.32;    -5.72;
     11.14;    28.30;    33.22;    30.27;    25.49;    18.45;     7.62;    -1.66;
     -7.86;   -10.04;    -5.71;    -5.16;   -12.56;   -19.39;   -22.16;   -18.18;
    -12.00;    -3.60;     3.99;     3.89;     0.85;     4.01;     6.05;     0.68;
     -0.14;     1.50;     3.17;     4.96;    -0.03;    -7.20;    -4.85;    -0.60;
     -3.05;    -6.01;    -6.08;    -6.40;    -6.70;    -5.70;    -6.14;    -8.41;
    -10.05;   -12.35;   -15.72;     0.00
     ], [cm/sec/sec] );

   dimen = Dimension(Elcentro1);
   x     = dimen[1][1];

/* Interpolate accelerations */

   divid_no = 16;
   Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] );
   Elcentro[1][1] = Elcentro1[1][1]/16;
   Elcentro[2][1] = Elcentro1[1][1]/8;
   Elcentro[3][1] = Elcentro1[1][1]*3/16;
   Elcentro[4][1] = Elcentro1[1][1]/4;
   Elcentro[5][1] = Elcentro1[1][1]*5/16;
   Elcentro[6][1] = Elcentro1[1][1]*3/8;
   Elcentro[7][1] = Elcentro1[1][1]*7/16;
   Elcentro[8][1] = Elcentro1[1][1]/2;
   Elcentro[9][1] = Elcentro1[1][1]*9/16;
   Elcentro[10][1]= Elcentro1[1][1]*5/8;
   Elcentro[11][1]= Elcentro1[1][1]*11/16;
   Elcentro[12][1]= Elcentro1[1][1]*3/4;
   Elcentro[13][1]= Elcentro1[1][1]*13/16;
   Elcentro[14][1]= Elcentro1[1][1]*7/8;
   Elcentro[15][1]= Elcentro1[1][1]*15/16;
   Elcentro[divid_no*x][1] = Elcentro1[x][1];
   for( i=1 ; i < x; i=i+1 ) {
      Elcentro[divid_no*i][1] = Elcentro1[i][1];
      Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*15+Elcentro1[i+1][1] )/16;
      Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1] )/8;
      Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*13+Elcentro1[i+1][1]*3 )/16;
      Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1] )/4;
      Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*11+Elcentro1[i+1][1]*5 )/16;
      Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*3 )/8;
      Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1]*9 +Elcentro1[i+1][1]*7 )/16;
      Elcentro[divid_no*i+8][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1] )/2;
      Elcentro[divid_no*i+9][1] = ( Elcentro1[i][1]*7 +Elcentro1[i+1][1]*9 )/16;
      Elcentro[divid_no*i+10][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*5 )/8;
      Elcentro[divid_no*i+11][1] = ( Elcentro1[i][1]*5 +Elcentro1[i+1][1]*11 )/16;
      Elcentro[divid_no*i+12][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*3 )/4;
      Elcentro[divid_no*i+13][1] = ( Elcentro1[i][1]*3 +Elcentro1[i+1][1]*13 )/16;
      Elcentro[divid_no*i+14][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*7 )/8;
      Elcentro[divid_no*i+15][1] = ( Elcentro1[i][1]   +Elcentro1[i+1][1]*15 )/16;
   }

/*
   divid_no = 8;
   Elcentro = ColumnUnits( Matrix([ divid_no*x, 1 ]),[in/sec/sec] );
   Elcentro[1][1] = Elcentro1[1][1]/8;
   Elcentro[2][1] = Elcentro1[1][1]/4;
   Elcentro[3][1] = Elcentro1[1][1]*3/8;
   Elcentro[4][1] = Elcentro1[1][1]/2;
   Elcentro[5][1] = Elcentro1[1][1]*5/8;
   Elcentro[6][1] = Elcentro1[1][1]*3/4;
   Elcentro[7][1] = Elcentro1[1][1]*7/8;
   Elcentro[divid_no*x][1] = Elcentro1[x][1];
   for( i=1 ; i < x; i=i+1 ) {
      Elcentro[divid_no*i][1] = Elcentro1[i][1];
      Elcentro[divid_no*i+1][1] = ( Elcentro1[i][1]*7+Elcentro1[i+1][1] )/8;
      Elcentro[divid_no*i+2][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1] )/4;
      Elcentro[divid_no*i+3][1] = ( Elcentro1[i][1]*5+Elcentro1[i+1][1]*3 )/8;
      Elcentro[divid_no*i+4][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1] )/2;
      Elcentro[divid_no*i+5][1] = ( Elcentro1[i][1]*3+Elcentro1[i+1][1]*5 )/8;
      Elcentro[divid_no*i+6][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1]*3 )/4;
      Elcentro[divid_no*i+7][1] = ( Elcentro1[i][1]  +Elcentro1[i+1][1]*7 )/8;
   }
*/

   x = Max( [ Max(Elcentro), abs(Min(Elcentro)) ] );
   Elcentro = (0.5*9.81/x)*Elcentro;

   print "\n*** Time History Analysis :";
   print "\n*** Use Average Acceleration Method :\n\n";

/* ----------------------------------------------------- */
/* Setup Time History Analysis  (and initial conditions) */
/* ----------------------------------------------------- */

/* Setup initial external load, internal load to zeros */

   NodeLoad( 1, [ 0 kN, 0 kN, 0 kN, 0 kN*m, 0 kN*m, 0 kN*m] );
   P_ext = ExternalLoad();
   P_old = P_ext;

/* Setup initial displacement, velocity and acceleration to zeros */

   displ    = Solve( stiff, P_ext );
   velocity = displ/(1 sec);
   accel    = velocity/(1 sec);

   new_displ    = displ;
   new_velocity = velocity;
   new_accel    = accel;
   displ_i      = displ;

/* Setup initial internal force and damping force */

   Fs_old = InternalLoad( displ );
   Fd_old = damp*velocity;
   Fs_i   = Fs_old;

/* Setup the influence vector */

   r = displ/(1 m);
   accel_dir = 3;
   for( i=1 ; i<=total_node ; i=i+1 ) {
      dof = GetDof([i]);
      if( dof[1][accel_dir] > 0 ) {
         r[ dof[1][accel_dir] ][1] = 1;
      }
   }

/* Setup time interval and analysis time */

   dt     = (0.02 sec)/divid_no;
   time   = 0.0  sec;
   tol    = 0.0001;
   stepno = 0;

   dimen = Dimension(Elcentro);
   total_stepno = dimen[1][1] + 250*divid_no;
   total_time   = total_stepno * dt;
   quake_time   = dimen[1][1] * dt;
   print "        dt =",dt,"\n";
   print "total step =",total_stepno," time       =",total_time,"\n";
   print "quake step =", dimen[1][1]," quake time =",quake_time,"\n";

/* Allocate response matrices */

   deck_displ     = ColumnUnits(Zero([total_stepno,5]),[m]);
   pier_displ     = ColumnUnits(Zero([total_stepno,3]),[m]);
   bottom_shear   = ColumnUnits(Zero([total_stepno,3]),[N]);
   isolator_shear = ColumnUnits(Zero([total_stepno,5]),[N]);

   energy  = ColumnUnits(Zero([total_stepno,4]),[Jou]);

   /* Time-History Analysis */

   while(stepno < total_stepno) {

      /* Update time and step no */

      time   = time + dt;
      stepno = stepno + 1;

      print "*** Start at step  ", stepno, " : TIME = ", time, "\n";

      /* Compute effective incremental loading */

      if( stepno <= dimen[1][1] ) then {
          P_ext = -mass*r*Elcentro[stepno][1];
      } else {
          P_ext = -mass*r*(0.0 m/sec/sec);
      }

      dPeff = P_ext - P_old + ((4/dt)*mass + 2*damp)*velocity + 2*mass*accel;

      /* while loop to check converge, Keff*U = Peff */

      dp  = displ - displ;
      err = 1 + tol;

      ii = 1;
      while( err > tol ) {

         /* Compute effective stiffness from tangent stiffness */

         Keff = stiff + (2/dt)*damp + (4/dt/dt)*mass;

         /* Solve for d_displacement, d_velocity */

         dp_i = Solve( Keff, dPeff);
         dp = dp + dp_i;
         dv = (2/dt)*dp - 2*velocity;

         /* Compute displacement, velocity */

         new_displ    = displ + dp;
         new_velocity = velocity + dv;

         /* Compute incremental displacement and internal load using old stiffness */

         dFs = stiff*dp_i;

         Fs_i = Fs_i + dFs;
         if( ii==1 ) {
             x = L2Norm(dFs);
             if( x == 0 ) { x = 1; }
         }

         /* Check material yielding and compute new stiffness */

         ElmtStateDet( dp_i );
         stiff = Stiff();

         /* Compute new internal load, damping force, and acceleration */

         Fs = InternalLoad( new_displ );
         Fd = damp*new_velocity;
         new_accel = mass_inv*( P_ext-Fs-Fd );

         /* Calculate the unbalance force, and error percentage */

         P_err = Fs_i - Fs;
         y     = L2Norm(P_err);
         err   = y/x;

         /* Assign new effective incremental load */

         dPeff   = P_err;
         displ_i = new_displ;
         Fs_i    = Fs;

         print "*** In While Loop" ,ii, ", dFs = ", L2Norm(dFs), ", P_err =" , y , ", err =" ,err, "\n";
         ii = ii+1;
         if( ii > 50 ) { flag = 1;
                          err = tol;
                       }
      }

      print "*** FINISHED CONVERGENCE AT STEP", stepno, "\n";

      /* Calculate the energy balance, Wint(n+1) + T(n+1) = Wext(n+1) */
      /*                                                              */
      /* [1] column 1 is internal work by stiffness                   */
      /* [2] column 2 is internal work by damping                     */
      /* [3] column 3 is kinetic energy                               */
      /* [4] column 4 is external work                                */

      trans_vel     = Trans(velocity);         /* V(i-1) */
      trans_vel_new = Trans(new_velocity);     /* V(i)   */
      trans_dis_new = Trans(new_displ);        /* D(i)   */

      if( stepno == 1 ) then {
          Wint_fs = dt/2*(trans_vel_new*Fs);
          Wint_fd = dt/2*(trans_vel_new*Fd);
          T    = (trans_vel_new*mass*new_velocity)/2;
          Wext = dt/2*(trans_vel_new*P_ext);
          energy[1][1] = Wint_fs[1][1];
          energy[1][2] = Wint_fd[1][1];
          energy[1][3] = T[1][1];
          energy[1][4] = Wext[1][1];
      } else {
          Wint_fs = dt/2*(trans_vel*Fs_old + trans_vel_new*Fs);
          Wint_fd = dt/2*(trans_vel*Fd_old + trans_vel_new*Fd);
          T    = (trans_vel_new*mass*new_velocity)/2;
          Wext = dt/2*(trans_vel*P_old + trans_vel_new*P_ext);
          energy[stepno][1] = energy[stepno-1][1] + Wint_fs[1][1];
          energy[stepno][2] = energy[stepno-1][2] + Wint_fd[1][1];
          energy[stepno][3] = T[1][1];
          energy[stepno][4] = energy[stepno-1][4] + Wext[1][1];
      }

      print "*** flag 1\n";

      /* Tolerance is satisfied, update histories for this time step */

      UpdateResponse();
      P_old    = P_ext;
      Fs_old   = Fs;
      Fd_old   = Fd;
      displ    = new_displ;
      velocity = new_velocity;
      accel    = new_accel;

      print "*** flag 2\n";

      /* Store node displacement for this time step */

      i = 1;
      for( node_no=1 ; node_no<=5 ; node_no=node_no+1 ) {
         dof = GetDof([node_no]);
         deck_displ[stepno][i] = displ[ dof[1][accel_dir] ][1];
         i = i+1;
      }
      i = 1;
      for( node_no=7 ; node_no<=9 ; node_no=node_no+1 ) {
         dof = GetDof([node_no]);
         pier_displ[stepno][i] = displ[ dof[1][accel_dir] ][1];
         i = i+1;
      }

      print "*** flag 3\n";

      /* Store element shear force for this time step */

      /* In the meantime, comment this section out .....
      i = 1;
      for( elmt_no=5 ; elmt_no<=9 ; elmt_no=elmt_no+1 ) {
           stress = GetStress( [elmt_no], displ );
           isolator_shear[stepno][i] = stress[1][accel_dir];
           i=i+1;
      }
      */

      print "*** flag 4\n";

      /*
      stress = GetStress([bot_elmt2],displ);
      bottom_shear[stepno][1] = -stress[2][accel_dir];
      stress = GetStress([bot_elmt3],displ);
      bottom_shear[stepno][2] = -stress[2][accel_dir];
      stress = GetStress([bot_elmt4],displ);
      bottom_shear[stepno][3] = -stress[2][accel_dir];
      */

      print "*** flag 5\n";

      if( flag == 1 ) {
          stepno = total_stepno;
      }

      print "deck_displ : ";
      print  deck_displ[stepno][1],deck_displ[stepno][2],deck_displ[stepno][3],
             deck_displ[stepno][4],deck_displ[stepno][5],"\n";
      print "pier_displ : ";
      print  pier_displ[stepno][1], pier_displ[stepno][2],
             pier_displ[stepno][3],"\n";
      print "isol_shear : ";
      print  isolator_shear[stepno][1], isolator_shear[stepno][2],
             isolator_shear[stepno][3], isolator_shear[stepno][4],
             isolator_shear[stepno][5],"\n";
      print "bottom_shear : ";
      print  bottom_shear[stepno][1], bottom_shear[stepno][2],
             bottom_shear[stepno][3], "\n";
      print "energy       : ";
      print  energy[stepno][1], energy[stepno][2],
             energy[stepno][3], energy[stepno][4],"\n";
   }

   PrintMatrix(deck_displ,pier_displ);
   PrintMatrix(isolator_shear,bottom_shear);
   PrintMatrix(energy);

   quit;


Developed in March 1998 by Mark Austin
Copyright © 1998-1999, Department of Civil Engineering, University of Maryland