[ Problem Statement ]
[ Finite Element Mesh ]
[ Dynamic Analysis ]
[ Displacements ]
[ Energy Balance Calculations ]
[ Input/Output Files ]
You will need Aladdin 2.0 to run this problem.
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 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.
The finite element mesh is shown in Figure 2.
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.
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.
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.
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; }
The first and second floor displacements in the x- and y- directions are plotted in Figures 5 and 6, respectively.
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.
Points to note:
References