Three-Dimensional Dynamic Analysis of Pier Structure

[ Problem Statement ] [ Finite Element Mesh ] [ Dynamic Analysis ]
[ Displacements ] [ Energy Balance Calculations ] [ Input/Output Files ]

You will need Aladdin 2.0 to run this problem.


PROBLEM STATEMENT

In this example we compute the nonlinear time-history response of a two story lumped mass pier subject to earthquake ground motions simultaneously applied in the x- and y- directions.

Figure 1 : Three-Dimensional Two-Story Lumped Mass Pier

Figure 1 is a three-dimensional view of the two-story lumped mass pier. The total mass of the first floor is 1500 kg. The second floor has mass 1000 kg.

We assume for the purposes of this analysis that each floor has the dimensions shown in the top right-hand corner of Figure 1. This effect includes the rotational inertias of both floors, namely:

    J_x  = L^2 / 12               J_y = B^2 / 12
    M_jx = M.J_x                 M_jy = M.J_y

The lumped mass matrix of each floor in local coordinates is:

    Mass = [ M_i  0.0  0.0  0.0;
             0.0  M_i  0.0  0.0;
             0.0  0.0  M_jx 0.0;
             0.0  0.0  0.0  M_jy ];

The Wilson-theta algorithm is used for numerical integration of the underlying equations of equilibrium.


FINITE ELEMENT MODEL

The finite element mesh is shown in Figure 2.

Figure 2 : Element Mesh for Two-Story Lumped Mass Pier

The finite element model has 3 nodes and 2 fiber elements. The boundary conditions are full fixity at the base. Axial deformations (u_z) and torsional rotations (theta_z) in each element are likewise assumed to be zero.

The pier is modeled with two "FIBER_3D" elements with 30 by 40 = 1200 fibers and 10 Gauss-Lobatto integration points for each element.

Figure 3 : Global DOF for Two-Story Lumped Pier

After the boundary conditions are applied, the model has only 8 degrees of freedom (d.o.f.). And this is in spite of the 2400 fibers used to model the pier deformations. See Figure 3.

Earthquake Accelerograms

The time-history response of the two-story lumped mass pier structure is computed for ground motion accelerograms simultaneously applied in the x and y directions. The ground motion accelerograms are ten second samples extracted from the 1940 El Centro record by Balling et al. [1]. See Figure 4.

Figure 4 : 1940 El Centro Earthquake Record

The standard way of solving this problem is to note that the total displacement "v^t" may be expressed as the sum of the relative displacement "v" and the pseudo-static displacements "v^s" that would result from a static-support displacement. That is

    v^t = v + v^s.

The pseudo-static displacements may be conveniently described by an influence vector "r" representing the displacements resulting from a unit support displacement. For multi-component time-history analysis. Using "r_x" to represent the influence coefficient vector "r" in the x direction, and "r_y" in the y direction, the effective-force vector generated by the earthquake ground motion components is

    P(t) = -[M][r][v_g(t)]
         = -[M]([r_x][v_{gx}(t)] + [r_y][v_{gy}(t)]).

If the structural degrees of freedom are as shown in Figure XXY, then a unit displacement in the x direction affects only d.o.f. 1 and 5, and a unit displacement in the y direction affects only d.o.f. 2 and 6. Hence

        [ r_x ] = [ 1 0 0 0 1 0 0 0 ]^T

    and [ r_y ] = [ 0 1 0 0 0 1 0 0 ]^T.


DYNAMIC ANALYSIS

The following abbreviated input file shows the essential details of computing the multi-component time-history analysis with the Wilson-theta method.

   /* Compute initial stiffness and mass matrices */

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

   /* Manually add lumped mass and rotational inertial to mass matrix. */

   dof = GetDof([2]);
   mass[dof[1][1]][dof[1][1]] = M1;
   mass[dof[1][2]][dof[1][2]] = M1;
   mass[dof[1][4]][dof[1][4]] = Mjx1;
   mass[dof[1][5]][dof[1][5]] = Mjy1;
   dof = GetDof([3]);
   mass[dof[1][1]][dof[1][1]] = M2;
   mass[dof[1][2]][dof[1][2]] = M2;
   mass[dof[1][4]][dof[1][4]] = Mjx2;
   mass[dof[1][5]][dof[1][5]] = Mjy2;

   mass_inv = Inverse(mass);

   /* Setup Rayleigh damping and damping matrix */

   rdamping = 0.05;
   A0 = 2*rdamping*w1*w2/(w1+w2);
   A1 = 2*rdamping/(w1+w2);
   damp = A0*mass + A1*stiff;

   /* Setup initial displacement, velocity and acceleration */

   NodeLoad( 1, [ 0 kN, 0 kN, 0 kN, 0 kN*m, 0 kN*m, 0 kN*m] );
   P_ext = ExternalLoad();
   displ = Solve( stiff, P_ext );
   velocity = displ/(1 sec);
   accel = velocity/(1 sec);

   /* 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; }
   }

   /* Define theta value for Wilson-theta method */

   theta = 1.420815;
   ds   = theta*dt;

   /* Wilson-theta time-history analysis */

   for( stepno=1 ; stepno <= total_stepno ; stepno=stepno+1 ) {

      /* [1] : Compute effective incremental loading */

      time = time + dt;
      if( time <= quake_time ) 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 = theta*(P_ext-P_old) + mass*((6/ds)*velocity+3*accel)
            + damp*(3*velocity+(ds/2)*accel);

      /* [2] : Compute effective stiffness from tangent stiffness */

      Keff = stiff + (3/ds)*damp + (6/ds/ds)*mass;

      /* [3] : Solve for estimated delta_displacement */

      dps = Solve( Keff, dPeff );

      /* [4] : Compute estimated displacement, velocity and acceleration */

      ds_accel = (6/ds/ds)*dps - (6/ds)*velocity - 3*accel;
      new_velocity = velocity + dt*accel + (dt/2/theta)*ds_accel;
      new_displ = displ + dt*velocity + (dt*dt/2)*accel
                  + (dt*dt/6/theta)*ds_accel;

      /* [5] : Check material yielding and compute new stiffness */

      dp = new_displ - displ;
      ElmtStateDet( dp );
      stiff = Stiff();

      /* [6] : Compute new internal load, damping force, acceleration */

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

      ...... details of energy balance calculation explained later ......

      /* [7] : Update histories for this time step */

      UpdateResponse();

      P_old    = P_ext;
      displ    = new_displ;
      velocity = new_velocity;
      accel    = new_accel;
   }


DISPLACEMENTS

The first and second floor displacements in the x- and y- directions are plotted in Figures 5 and 6, respectively.

Figure 5 : Displacement in x-direction (m) versus time (sec)

Figure 6 : Displacement in y-direction (m) versus time (sec)


ENERGY BALANCE CALCULATIONS

We now extend the previous two-story lumped mass pier example with the energy balance calculations. At the highest level of abstraction we require:

        W_internal + T = W_external

Here

        W_internal = work done on the internal system.

        T          = kinetic energy of the system.

        W_external = work done by the external loads.

In other words, the work done by external loads is converted to either kinetic energy or work done on the internal system. The latter can be stored either elastically or dissipated by plastic deformations.

This analysis uses a disscrete approximation to the energy balance equation, which includes the work done by the nodals loads that are developed from the straining of materials, damping forces, and so forth. The complete derivation can be found in Chapter 5 of Wane-Jang Lin's Ph.D. dissertation [2].

The following abbreviated input file is positioned after the end of each load step calculation, but before the updating of the element history response.

   /* Calculate the energy balance, Wint(n+1) + T(n+1) = Wext(n+1) */

   trans_vel     = Trans(velocity);
   trans_vel_new = Trans(new_velocity);
   trans_dis_new = Trans(new_displ);

   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*Ps_int + trans_vel_new*Fs);
       Wint_fd = dt/2*(trans_vel*Pd_int + 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];
   }

The matrix "energy" stores the results for each time step. The first column contains the internal energy done by straining force, the second column contains the internal energy done by damping force, the third column stores the kinetic energy, and the fourth column stores the external energy.

The overall results of our energy balance calculation are plotted in Figures 7 and 8.

Figure 7 : Internal Energy (Joules) versus Time (sec)

Figure 8 : External/Internal Energy (Joules) versus Time (sec)

Points to note:

  1. Figure 7 shows a great match of the internal energy and the external energy, therefore proving that the analysis of this nonlinear dynamic problem has stable and accurate computation.

  2. Note from Figure 8 that after the earthquake ground motions have ceased, there is no more external energy input into the system. Hence, the curve W_ext versus time is constant over the response interval t equals [10, 15] seconds.

References

  1. Balling R.J., Pister K.S., Polak E., DELIGHT.STRUCT : A Computer-Aided Design Environment for Structural Engineering , Computer Methods in Applied Mechanics and Engineering, January 1983, pp. 237-251.
  2. Lin W.J., Modern Computational Environments for Seismic Analysis of Highway Bridge Structures, Ph.D. Dissertation, University of Maryland, College Park, MD 20742, December 1997, (p. 197).


INPUT/OUTPUT FILES


Developed in July 1997 by Wane-Jang Lin and Mark Austin
Last Modified March 18, 1998
Copyright © 1998, Mark Austin, Department of Civil Engineering, University of Maryland