input-biaxial

有限元程序

    -40.02;   -43.52;   -27.74;   -11.32;     9.04;    22.97;    39.00;    41.79;
     36.32;    33.21;    24.78;    19.35;    25.37;    41.24;    29.07;    32.97;
     22.07;    24.15;   -43.88;   -37.85;   -43.29;   -47.30;   -41.99;   -41.04;
    -61.55;   -50.26;   -45.22;   -35.65;   -25.37;   -11.40;    -0.30;    16.00;
     26.08;    26.63;    25.93;    25.26;    23.29;    30.02;    23.56;    19.43;
     16.99;    19.74;    23.29;    19.11;     6.87;    -3.88;    -6.44;     0.69;
     19.86;    38.21;    11.48;    -8.37;     3.37;    16.91;    26.95;    40.53;
     46.63;    42.07;    19.55;    -0.57;     0.57;    -6.75;    -5.45;   -14.55;
    -31.95;   -40.61;   -40.22;   -37.89;   -28.13;   -18.60;    -9.55;    -0.57;
      6.56;    21.52;    19.67;     2.07;    -4.15;   -17.11;   -22.22;   -20.85;
    -13.09;    -5.41;   -18.72;   -32.42;   -19.78;   -17.54;    -6.67;     2.62;
     11.48;    -0.85;   -12.54;   -17.18;   -10.49;    -7.81;     1.40;     8.09;
     14.70;    12.11;     4.04;    -3.13;   -11.91;   -18.52;    -6.16;     3.41;
     17.34;    23.84;     7.11;    -3.76;   -17.78;   -27.26;   -23.21;   -20.93;
    -15.14;    -8.05;    -1.44;     5.10;     3.96;     4.27;     3.37;    10.96;
     10.10;     4.82;     2.11;    -4.23;    -8.37;   -15.65;   -12.70;    -1.75;
      9.55;     8.41;     3.05;     2.38;    -5.30;   -10.69;   -17.26;   -21.95;
    -20.69;   -29.55;   -21.83;   -11.48;     3.92;    26.04;    33.13;    23.01;
     15.26;    -4.82;   -22.85;   -37.81;   -34.82;   -24.00;   -16.83;   -27.15;
    -29.51;   -14.67;    -2.93;     9.07;    21.08;    30.30;    36.59;    32.46;
     24.78;    18.37;     9.94;     3.21;    -8.13;   -19.98;   -30.30;   -21.59;
    -15.77;    -5.69;   -11.04;   -13.72;   -19.78;   -16.55;   -10.37;    -2.30;
      7.66;    18.17;    23.64;    16.87;    12.66;     6.59;     1.24;    -6.16;
    -13.01;   -16.52;   -13.84;    -8.33;   -10.26;   -15.57;   -19.43;   -24.00;
    -28.41;   -31.12;   -24.19;   -18.33;   -15.69;   -26.12;   -27.30;   -21.91;
    -15.14;   -10.14;     0.77;    -8.64;   -18.05;    -7.11;    -6.08;     1.87;
      8.25;    -1.16;   -10.81;   -20.02;   -30.85;   -40.65;   -51.55;   -39.04;
    -12.18;    15.26;    37.70;    29.86;    43.56;    46.91;    42.74;    14.86;
     24.47;    19.70;    32.11;    -1.20;   -10.45;    -6.52;   -14.51;   -14.27;
     -8.88;     5.06;    17.42;    22.46;    11.16;     6.83;     1.52;    14.07;
     31.67;    37.97;    40.96;    41.99;    29.74;    26.83;    27.93;    28.05;
     29.35;    31.87;    32.30;    34.51;     7.46;     0.81;     4.23;     1.63;
      5.93;     7.58;    11.16;     9.98;     1.28;    -5.49;   -15.41;   -19.39;
    -14.74;   -14.74;   -12.18;   -11.99;   -26.44;   -38.05;   -30.49;   -20.33;
     -3.21;    12.58;    16.87;     2.26;   -14.31;   -26.32;   -41.48;   -29.90;
    -19.90;     3.45;    11.59;     8.72;     4.43;     5.10;     9.07;    13.60;
     10.77;    -2.11;   -13.33;   -28.80;   -33.92;   -38.17;   -17.62;   -22.54;
    -39.82;   -17.58;     2.50;    22.62;    36.63;    37.66;    16.52;    -2.07;
    -15.96;    -4.43;     2.22;    13.05;    28.13;    17.19;   -16.40;    -2.62;
    -10.85;    -2.62;   -10.61;   -18.05;   -12.30;   -10.81;    -6.67;    -2.07;
      0.57;    -4.74;   -10.02;   -17.70;   -26.55;   -19.78;    -5.53;    11.36;
     22.66;     3.68;     5.61;   -12.26;     0.33;    10.45;    18.64;    15.45;
      2.42;    -9.98;   -15.10;    -3.68;     9.27;    29.90;    26.67;    19.19;
     11.28;     1.59;    -6.83;   -24.23;   -40.93;   -55.73;   -53.44;   -47.50;
    -29.07;   -26.04;   -30.85;   -23.84;   -13.17;    -3.84;    13.92;    31.91;
     47.19;    55.53;    65.77;    73.52;    81.48;    76.87;    70.65;    63.80;
     55.73;    47.19;    44.82;    46.08;    26.87;    13.92;    -4.59;   -20.89;
    -37.93;   -54.11;   -66.44;     0.
    ], [in/sec/sec] );

/* Setup initial external load, internal load and error pseudo-load */

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

/* Setup initial displacement, velocity and acceleration */

   displ = Solve( stiff, P_ext-P_ext );
   velocity = displ/(1 sec);
   mass_inv = Inverse(mass);
   accel = mass_inv*(P_ext-damp*velocity-stiff*displ);

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

/* Setup initial internal force and damping force */

   Fs = InternalLoad( new_displ );
   Fd = damp*velocity;

/* Setup the influence vector in both dir-X and dir-Y */

   rx = displ/(1 m);
   ry = displ/(1 m);
   accel_dir_x = 1;
   accel_dir_y = 2;
   for( i=1 ; i<=total_node ; i=i+1 ) {
      dof = GetDof([i]);
      if( dof[1][accel_dir_x] > 0 ) { rx[ dof[1][accel_dir_x] ][1] = 1; }
      if( dof[1][accel_dir_y] > 0 ) { ry[ dof[1][accel_dir_y] ][1] = 1; }
   }

/* Setup time interval and analysis time */

   divid_no = 1;
   dt   = 0.02 sec;
   time   = 0.0  sec;
   stepno = 0;

   dimen = Dimension(ground_accel_x);
   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,"\ttime =",total_time,"\n";
   print "quake step =",dimen[1][1],"\tquake time =",quake_time,"\n";

/* Allocate response matrices */

   deck_displ_x  = ColumnUnits(Zero([total_stepno,2]),[m]);
   deck_force_x  = ColumnUnits(Zero([total_stepno,2]),[N]);
   deck_displ_y  = ColumnUnits(Zero([total_stepno,2]),[m]);
   deck_force_y  = ColumnUnits(Zero([total_stepno,2]),[N]);

/* Using average acceleration method */

   /* Time-History Analysis */

   while(stepno < total_stepno) {

      /* Update time and step no */

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

      /* Compute effective incremental loading */

      if( stepno <= dimen[1][1] ) then {
         P_ext = -mass*(rx*ground_accel_x[stepno][1]+ry*ground_accel_y[stepno][1]);
      }
      else { P_ext = -mass*(rx*(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 */
      tol = 0.001;
      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, ", P_err =" ,x, ", err =" ,err, "\n";
         ii = ii+1;
         if( ii > 10 ) { flag=1; err=tol; }
      }

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

      UpdateResponse();
      P_old    = P_ext;
      displ    = new_displ;
      velocity = new_velocity;
      accel    = new_accel;

      i = 1;
      for( node_no=2 ; node_no<=3 ; node_no=node_no+1 ) {
         dof = GetDof([node_no]);
         deck_displ_x[stepno][i] = displ[ dof[1][accel_dir_x] ][1];
         deck_force_x[stepno][i] = Fs[ dof[1][accel_dir_x] ][1];
         deck_displ_y[stepno][i] = displ[ dof[1][accel_dir_y] ][1];
         deck_force_y[stepno][i] = Fs[ dof[1][accel_dir_y] ][1];
         i = i+1;
      }
      if( flag==1 ) { stepno = total_stepno; }
   }

PrintMatrix(deck_displ_x,deck_force_x);
PrintMatrix(deck_displ_y,deck_force_y);
quit;