Force-Displacement Response of Mass-Spring System

[ Problem Statement ] [ Load-Displacement Curve ] [ State Determination ]
[ Force-Displacement Calculation ] [ Strain Energy Calculation ]
[ Time-History Analysis ] [ Input/Output Files ]

You will need Aladdin 2.0 to run this problem.


PROBLEM STATEMENT

Consider the mass-spring system shown in Figure 1.

Figure 1 : Nonlinear Mass-Spring System

In this example we use Aladdin's matrix language to calculate the load-displacement response of a nonlinear mass-spring system subject to a well-defined external loading. Our analysis will be divided into two parts:

  1. First we assume that the externally applied loads are static, and the underlying equations of equilibrium are not influenced by damping or inertial effects.

  2. The second part of our analysis will account for dynamic effects -- that is, our equations of equilibrium will include inertia and damping force effects.

By using the matrix language for the complete analysis, we are really creating a high-level software prototype suitable for testing algorithms and new ideas. This particular example assumes bi-linear hysteretic behavior of the springs. (A reasonable variation on this problem would be to examine the behavior of a system having a Romberg-Osgood stress-strain curve.)

Once everything works, the follow-up implementation would involve a finite element implementation of the modeling mesh and element-level behavior. The second implementation will have improved performance because sections of interpreted code in the first version will now be executed as low-level compiled code in the finite element version.

Section and Material Properties

The lower-half of Figure 1 defines the mass-spring system properties.

Springs 1 and 2 both have a bi-linear force-displacement relationships which follow the kinematic hardening rule . The basic shape of the force-displacement constituitive relationship is defined by the Aladdin variables:

   /* Define bi-linear spring behavior */

   ks1 = 2.0 N/cm;   ks2 = 1.5 N/cm;
   kt1 = 0.8 N/cm;   kt2 = 0.5 N/cm;
   fy1 =  18 N;      fy2 =  15 N;

   Ks = [ks1; ks2];
   Kt = [kt1; kt2];
   Fy = [fy1; fy2];

The undeformed springs start out with an initial stiffness "ks" that lasts until the load exceeds the yield force "fy". The tangent stiffness then changes to "kt". When the loading is reversed, the tangent stiffness switches back to the initial stiffness "ks" until the yielding reappears. The elastic force-displacement range is 2.fy.

External Loads

At each step of the analysis the end-load is incremented by either 1 N or -1 N, depending on whether the load is increasing or decreasing.

The relevant section of Aladdin code is:

   for ( k = 1; k <= total_step ; k = k+1 ) {

      /* define force increment for each step */

      if( k<=20 )           { d_P =  [0 N; 1 N]; }
      if( k>20  && k<=60  ) { d_P = -[0 N; 1 N]; }
      if( k>60  && k<=105 ) { d_P =  [0 N; 1 N]; }
      if( k>105 && k<=155 ) { d_P = -[0 N; 1 N]; }
      if( k>155 && k<=190 ) { d_P =  [0 N; 1 N]; }
      if( k>190 && k<=200 ) { d_P = -[0 N; 1 N]; }
      if( k>200 )           { d_P =  [0 N; 0 N]; }

      P = P + d_P;

      .... details of analysis removed .....
   }

Here d_P is the load increment and P is the external load. The resulting graph of external load versus step number is shown in Figure 2.

Figure 2 : Applied Load "P" versus Step Number


LOAD-DISPLACEMENT CURVE

Figure 3 shows the load-deflection curve for springs 1 and 2. Similar force-displacement curves can be drawn for the individual springs of course.

Figure 3 : Load-Deflection Curve for Springs 1 and 2

During the first 20 steps of the analysis the external load is incrementally increased from 0 N to 20 N. The system response will be elastic for the first 15 steps -- that is until the external load equals 15 N. The effective stiffness of the mass-spring system prior to first yield is:

    Effective stiffness =  1 / ( 1/ks1 + 1/ks2 )
                        =  6/7 N/cm.

Hence we expect that when the external force equals 15 N, the system displacement will be:

    Displacement at first yield = 15 N / 6/7 N/cm
                                = 17.5 cm.

For external loading between 15 and 18 N, the effective stiffness will be:

    Intermediate effective stiffness =  1/ ( 1/ks1 + 1/kt2 )
                                     =  2/5 N/cm.

Both springs will yield when the external load reaches 18 N. The effective stiffness drops to 4/13 N/cm for externalloads between 18 and 20 N.

From this point on the load-deflection is as shown in Figure 3, but defined by the loading history shown in Figure 2. When the analysis ceases at Step 200, the system will have a -12 cm residual displacement.


STATE DETERMINATION FOR NONLINEAR SPRINGS

Figure 4 shows the notation that we will use for the nonlinear spring state determination.

Figure 4 : Nomenclature for State Determination

We will use the variables "s" and "e" for the force and displacement at the integration section.

Variable(s)
Description
DRx
dx
Current force and displacement
s0
e0
Force and displacement at the point where the two tangent stiffness lines meet, like point A in Figure 3.
sr
er
Force and displacement at the point where the last displacement reversal with force, like point B in Figure 3.
sx
ex
Temporary force and displacement variables used in evaluation of section forces.
sx_saved
ex_saved
sr_saved
er_saved
Working arrays that store response from previous analysis step.

Table 1 : Arrays used in the State Determination

For the bi-linear force-displacement curve:

    DRx = sx + kx*( dx - ex ) 

The state determination can be viewed as a sequence of transitions from lines 0 through 4.

  1. Elastic Range : When the element deformations are still within the elastic range, and the material does not have any plastic residual, array element yielding[ele][1] = 0.
    if( yielding[ele][1] == 0 ) then {
    
        if( abs(dx) <= ey[ele][1] ) then {  /* line 0 on Figure 3         */
            kx = Ks[ele][1];
            sx = 0 N;
            ex = 0 cm;
            yielding[ele][1]  = 0;
            pre_range[ele][1] = 0;
        } else {                            /* move from line 0 to line 1 */
            kx = Kt[ele][1];
            s0[ele][1] = Fy[ele][1]*loading[ele][1];
            e0[ele][1] = ey[ele][1]*loading[ele][1];
            sx = s0[ele][1];
            ex = e0[ele][1];
            yielding[ele][1]  = 1;
            pre_range[ele][1] = 1;
        }
    
    } else { 
    
        .... source code for plastic residual .....
    
    }
    

    The two cases to consider are: (1) the response remains on line 0, and (2) the response moves from line 0 to either of lines 1 or 3.

    The array loading[ele][1] is given by:

         if( d_dx >  0 cm ) { loading[ele][1] =  1; }
         if( d_dx <  0 cm ) { loading[ele][1] = -1; }
         if( d_dx == 0 cm ) { loading[ele][1] = pre_load[ele][1]; }
    

    thereby allowing for yielding in both tension and compression.

  2. Plastic Range : The element has yielded (i.e., yielding[ele][1] = 1). A residual deformation will exist.
    if( yielding[ele][1] == 0 ) then {
    
        .... source code for elastic response .....
    
    } else { 
    
        if( pre_load[ele][1] != loading[ele][1] ) then {
    
            if( pre_range[ele][1] == 1 ) then {
                sr[ele][1] = sx_saved[ele][1];
                er[ele][1] = ex_saved[ele][1];
                kx = Ks[ele][1];
                sx = sr[ele][1];
                ex = er[ele][1];
                pre_range[ele][1] = 0;
            } else {   /* pre_range[ele][1] = 0 */
                if( loading[ele][1]*dx <= loading[ele][1]*er[ele][1] ) then {
                    kx = Ks[ele][1];
                    sx = sr[ele][1];
                    ex = er[ele][1];
                    pre_range[ele][1] = 0;
                } else {
                    if( loading[ele][1]*ex_saved[ele][1] >= loading[ele][1]*er[ele][1] ) then {
                        kx = Ks[ele][1];
                        sx = sr[ele][1];
                        ex = er[ele][1];
                        pre_range[ele][1] = 0;
                    } else {
                        kx = Kt[ele][1];
                        sx = sr[ele][1];
                        ex = er[ele][1];
                        pre_range[ele][1] = 1;
                    }
                }
                }
    
        } else {   /* pre_load[ele][1] == loading[ele][1] */
    
            if( pre_range[ele][1] == 1 ) then {
                kx = Kt[ele][1];
                sx = s0[ele][1];
                ex = e0[ele][1];
                pre_range[ele][1] = 1;
            } else {  /* pre_range[ele][1] = 0 */
    
                if( loading[ele][1]*ex_saved[ele][1] <= loading[ele][1]*er[ele][1] ) then {
                    if( loading[ele][1]*dx <= loading[ele][1]*er[ele][1] ) then {
                        kx = Ks[ele][1];
                        sx = sr[ele][1];
                        ex = er[ele][1];
                        pre_range[ele][1] = 0;
                    } else { /* line 2 -> line 1 */
                        kx = Kt[ele][1];
                        sx = sr[ele][1];
                        ex = er[ele][1];
                        pre_range[ele][1] = 1;
                    }
                } else { /* line 4 */
                    if( abs(dx-er[ele][1]) <= 2*ey[ele][1] ) then {
                        kx = Ks[ele][1];
                        sx = sr[ele][1];
                        ex = er[ele][1];
                        pre_range[ele][1] = 0;
                    } else { /* line 4 -> line 1 */
                        s0[ele][1] = sr[ele][1] + loading[ele][1]*2*Fy[ele][1];
                        e0[ele][1] = er[ele][1] + loading[ele][1]*2*ey[ele][1];
                        kx = Kt[ele][1];
                        sx = s0[ele][1];
                        ex = e0[ele][1];
                        pre_range[ele][1] = 1;
                    }
                } /* end of line 4 */
            }
        }
    }
    

Points to note:

  1. The array
        pre_load[][]
    

    stores the values of array loading from the previous analysis step.

  2. The array
        pre_range [][]
    

    indicates the elastic or plastic status of the previous step. The array element is set to 0 when the previous step in the stress-strain relationship is in elastic range. Conversely, its value is set to 1 if it was in the plastic range.

    So why do we need array pre_range[][]?

    If the previous step results in plastic range, and in the current step we have a reversal loading case, you know that the tangent stiffness for current step is elastic stiffness. But if it's elastic in previous step, even you have a reversal loading, you can't guarantee that it's elastic or plastic stiffness for this step, because it depends on the size of the displacement increment.


FORCE-DISPLACEMENT CALCULATION

The heart of the force-displacement calculation is contained in the block of code:

/* Define bi-linear spring behavior */

   .... details shown above .....

/* Compute and store displacements for incipient yielding of springs */

   ey = Matrix([2,1]);
   for( ele=1; ele<=2; ele=ele+1 ) {
      ey[ele][1] = Fy[ele][1]/Ks[ele][1];
   }

/* Initial condition, unstressed */

   P   = [0 N;   0 N];       /* structure force           */
   p   = [0 cm; 0 cm];       /* structure displacement    */
   PR  = [0 N;   0 N];       /* structure resisting force */
   PRe = [0 N;   0 N];       /* element resisting force   */

   Q_saved  = [0  N; 0  N];  /* matrix of element forces        */
   q_saved  = [0 cm; 0 cm];  /* matrix of element displacements */

/* initialize flags, stresses and strains before applying any load step */

   index = [0;0];  /* index for loading flag   */
   flag1 = [0;0];  /* yielding at load step k  */
   flag2 = [0;0];  /* pre_range at load step k */
   flag3 = [0;0];  /* pre_load at load step k  */

   yielding  = [0;0];
   pre_range = [0;0];
   pre_load  = [0;0];
   loading   = [0;0];
   sr = [0 N;   0 N];
   er = [0 cm; 0 cm];
   s0 = [0 N;   0 N];
   e0 = [0 cm; 0 cm];

   sr_saved = [0 N;   0 N];
   er_saved = [0 cm; 0 cm];
   s0_saved = [0 N;   0 N];
   e0_saved = [0 cm; 0 cm];
   sx_saved = [0 N;   0 N];
   ex_saved = [0 cm; 0 cm];

/* transformation matrix Lele, d_q = Lele*d_p[ele] for each element    */
/* transform structural displacements d_p to element deformations d_q  */

   Rigid = [-1,1];
   Transform = [1,0;0,1];
   L = Rigid*Transform;

/* temporary matrix to store element tangent stiffness Ks at each step */

   tangent = [Ks[1][1] ; Ks[2][1]];

/* force interpolation matrices b(x), D(x) = b(x)*Q */
/* relate section force D(x) with element force Q */

   bx = 1;

/* assemble initial structure tangent stiffness matrix BigK */

   BigK = [ Ks[1][1]+Ks[2][1], -Ks[2][1];
                    -Ks[2][1],  Ks[2][1] ];

/* Increase external load, structure determination */

tol = 0.00001; total_step = 400;
for ( k=1; k<=total_step ; k=k+1 ) {

   /* Define incremental loading */

      .... details shown above ....

   err = tol+1;
   while( err > tol )  { /* i-th Newton-Raphson iteration */

      /* compute displacement increment */

      d_p = Solve( BigK, d_P);
      p   = p + d_p;

      /* Element level state determination */

      for( ele = 1; ele <= 2; ele = ele + 1 ) {

         /* Retrieve element level displacements   */

         if( ele==1 ) {
             d_pe = [       0 m; d_p[1][1] ];
         }
         if( ele==2 ) {
             d_pe = [ d_p[1][1]; d_p[2][1] ];
         }

         /* Retrieve values from previous iteration */

         Q  = Q_saved[ele][1];      /* element-level forces           */
         q  = q_saved[ele][1];      /* element-level displacements    */
         Dx = bx*Q;                 /* section forces                 */
         dx = bx*q;                 /* section displacements          */
         K  = tangent[ele][1];      /* tangent stiffness of element   */
         kx = tangent[ele][1];      /* tangent stiffness of element   */
         fx = 1/kx;                 /* tangent flexibility of element */

         rx  = 0 cm;
         DUx = 1E+7 N;

         /* Increment element deformation */

         d_q = QuanCast(L*d_pe);
         q   = q + d_q;

         /* iterative loop for convergence of element forces */

         while( abs(DUx) > 0.00001 N ) {

            d_Q = K*d_q;          /* element force increment */
            Q   = Q + d_Q;        /* update element force    */

            /* determine the section force increments           */
            /* repeat for all integration points of the element */

            d_Dx = bx*d_Q;        /* section force increment           */
            d_dx = rx + fx*d_Dx;  /* section deformation increment     */
            Dx = Dx + d_Dx;       /* update section force              */
            dx = dx + d_dx;       /* update section deformation vector */

            /* Section level state determination */

               ... compute "kx", "sx", and "ex" ....

            /* Calculate stress and flexibility  */

            fx  = 1/kx;             /*                              */
            DRx = sx + kx*(dx-ex);  /*                              */

            DUx = Dx - DRx;         /* section unbalanced force     */
            rx  = fx*DUx;           /* section residual deformation */

            /* finish for all integration points of the element      */
            /* update the element flexibility and stiffness matrices */

            F = bx*fx*bx;
            K = 1/F;

            /* check for element convergence */

            s = bx*rx;   /* element residual deformation */
            d_q = -s;
         }  

         /* Save element level forces and displacements */

         Q_saved[ele][1] = Q;
         q_saved[ele][1] = q;

         tangent[ele][1] = K; /* element stiffness Ke =LT*K*L=[K,-K;-K,K]    */
         PRe[ele][1] = Q;     /* element resistant force PRe = LT*Q = [-Q;Q] */
      } 

      /* assemble structure resistant force */

      PR[1][1] = PRe[1][1] - PRe[2][1];
      PR[2][1] = PRe[2][1];

      /* assemble new structure stiffness */

      BigK[1][1] =  tangent[1][1] + tangent[2][1];
      BigK[1][2] = -tangent[2][1];
      BigK[2][1] = -tangent[2][1];
      BigK[2][2] =  tangent[2][1];

      d_P = P - PR;
      err = L2Norm(d_P);

   }  /* i-th iteration in Newton-Raphson while loop */
} 


STRAIN ENERGY CALCULATION

Figure 5 shows the elastic and inelastic strain energy in a typical nonlinear system.

Figure 5 : Elastic and Inelastic Strain Energy

During loading, the work done is equal to the area under the curve (i.e., area OABCDO). When the load is removed the load-deflection diagram follows line BD. A permanent elongation OD remains when point B is beyond the elastic limit. Thus, the elastic strain energy recovered during unloading is represented by the shaded triangle BCD. Area OABDO represents energy that is lost in the process of permanently deforming the bar. This energy is known as the inelastic (or plastic) strain energy.

Summary of Results

Figure 6 shows components of strain-enegy for springs 1 and 2. Similar force-displacement curves can be drawn for the individual springs.

Figure 6 : Strain Energy Plot for Springs 1 and 2

Solution Procedure

In the static analysis, a Newton-Raphson procedure is used to calculate the load-deflection response of this spring system. We also calculate the elastic and plastic components of energy of the spring elements. Except for the spring properties and strain energy calculation, the input file is almost the same as for the composite bar example presented in the preceding chapter.

The abbreviated input file is:

/* allocate the matrices for storing energy calculations */

   system_energy  = ColumnUnits( Matrix([total_step+1,4]), [Jou] );
   element_energy = ColumnUnits( Matrix([total_step+1,4]), [Jou] );

/* Increase External Load */

   for ( k = 1; k <= total_step ; k = k+1 ) {

   /* define force increment for each step */

      .... details of load increment found above ....

      element_energy[k+1][1] = element_energy[k][1];
      element_energy[k+1][3] = element_energy[k][3];

   /* i-th Newton-Raphson Iteration */

      index[1][1] = 1;   index[2][1] = 1;
      err = tol+1;
      while( err > tol ) {

         d_p = Solve(BigK,d_P);
         p = p + d_p;

         /* state determination for each element */

         for( ele=1;ele<=2;ele=ele+1 ) {

         ...... details about retrieving data from (j-1) ......

         q = q + d_q;

         /* element converge, j */

            while( abs(DUx) > 0.00001 N ) {

              ...... details of checking element convergence  ...... }
              ...... details of updating data at loop j ......

             /* energy calculations for each element ele */

            element_energy[k+1][2*ele-1] = element_energy[k+1][2*ele-1]
                           + 0.5*(Q_saved[ele][1]+Q)*(q-q_saved[ele][1]);
            element_energy[k+1][2*ele]   = 0.5*Q*Q/Ks[ele][1];
         } 

      /* assemble structure resistant force */

         ...... details of assembling structure resistant force ......

      /* assemble new structure stiffness */

         ...... details of assembling new structure stiffness ......

      d_P = P - PR;
      err = L2Norm(d_P);

   }  /* i-th iteration in Newton-Raphson while loop */

      ...... details of storing response removed ......

      ...... details of updating element history ......

   /* Reassemble the System Energy */

      system_energy[k+1][1] = system_energy[k][1] +
                              0.5*(result[k+1][1]+result[k][1])*(result[k+1][2]-result[k][2]);
      system_energy[k+1][2] = element_energy[k+1][1]+element_energy[k+1][3];
      system_energy[k+1][3] = element_energy[k+1][2]+element_energy[k+1][4];
      system_energy[k+1][4] = system_energy[k+1][2] - system_energy[k+1][3];
   }  

Points to note:

  1. In the input file, the matrix "element_energy" stores the element energy for all loading steps. The first column stores the internal work of element 1. The second column stores the element elastic strain energy of element 1. The third column stores the internal work of element 2. The fourth column stores the element elastic strain energy of element 2.

  2. A second matrix "system_energy" stores the system energy for all of the loading steps. The first column stores the total external work done by the system. The second column stores the total internal work of the system (it is the sum of internal work in the elements). The third column stores the total elastic strain energy of the system (it is the sum of elastic strain energy in the elements). The fourth column stores the total plastic strain energy of the system (again, this quantity is the sum of plastic strain energies in the elements).

  3. In the energy plots, you should observe that the plastic energy dissipation remains constant during periods of load until another yield point is reached.


TIME-HISTORY ANALYSIS

We now repeat the system analysis assuming that the load is incremented in time steps of only 0.01 seconds. The total loading time is 2 seconds. The system response is computed for a total of 4 seconds using the Newmark Integration algorithm with average acceleration across the time steps.

Figure 7 : Static/Dynamic Displacement-Responses

Figure 7 is a comparison of the static and dynamic displacement responses versus step no/time.

First Mode of Vibration

Figure 8 shows the change in the first-mode period versus time during the dynamic analysis.

Figure 8 : First mode dynamic period (sec) versus time (sec)

The time-variation in the first and second modes of vibration is computed by simply inserting:

   eigen = Eigen(BigK, BigM, [2]);
   eigenvalue  = Eigenvalue(eigen);

   result[1][7] = 2*PI/sqrt( eigenvalue[1][1] );
   result[1][8] = 2*PI/sqrt( eigenvalue[2][1] );

into the main loop of the time-history calculation. The 7th and 8th columns of "response" store the first and second mode periods versus step no.

System Energy

We also calculate the strain energy, elastic and plastic energies of the spring elements, and the natural periods of the system throughout the dynamic analysis. See Figure 9.

Figure 9 : Components of System Energy (Joules) versus Time (sec)


INPUT/OUTPUT FILES

Static Analysis

Dynamic Analysis


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