input-springs3

有限元程序

/*
 *  =====================================================================
 *  Plastic Analysis for a System of Two 1-DOF Spring Elements
 *
 *     -- general matrix format for bi-linear material
 *     -- energy calculation for balance checking
 *
 *  Dynamic analysis using average acceleration method.
 *
 *  Written By : Wane-Jang Lin                              October, 1997
 *  =====================================================================
 */

/* plastic spring material, Bi-linear */

   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;

   m1  = 1 kg;        m2  = 0.5 kg;
   c1  = 2 N*sec/m;   c2  = 1 N*sec/m;  /* damping coefficient */

   Ks = [ks1; ks2];
   Kt = [kt1; kt2];
   Fy = [fy1; fy2];
   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];
   q_saved  = [0 cm; 0 cm];

/* Initial flags and stresses and strains, initialed before 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 Ks at each step */

   tangent = [ks1 ; ks2];

/* 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 = [ ks1+ks2, -ks2;  -ks2, ks2 ];

/* assemble mass matrix and damping matrix */

   BigM = [    m1, 0 kg;  0 kg, m2 ];
   BigC = [ c1+c2,  -c2;   -c2, c2 ];

   M_inv = Inverse(BigM);

/* initial time = 0 sec conditions */

   total_step = 400;
   dt         = 0.01 sec;
   time       = 0 sec;
   print "\n\ndt =",dt,"\ttotal time=",total_step*dt,"\n";

/* Setup initial velocity and acceleration */

   vel   = p/(1 sec);
   accel = vel/(1 sec);

   new_p     = p;
   new_vel   = vel;
   new_accel = accel;

/* Setup initial internal force and damping force */

   Fs = P;
   Fd = P;

   P_old  = P;
   Fs_old = P;
   Fd_old = P;

/* Allocate the matrices for storing output history        */
/* column[1] : external applied load at end node (2)       */
/* column[2] : total elogation at end node (2) = [3]+[4]   */
/* column[3] : element elogation 1, node 1                 */
/* column[4] : element elogation 2, node 2                 */
/* column[5] : element force 1, node 1                     */
/* column[6] : element force 2, node 2                     */
/* column[7] : dynamic natrual period 1 (T1) for each step */
/* column[8] : dynamic natrual period 2 (T2) for each step */

   result = ColumnUnits( Matrix([total_step+1,8]), [N,cm,cm,cm,N,N,sec,sec] );
   result[1][1] = P[2][1];
   result[1][2] = p[2][1];
   result[1][3] = p[1][1];
   result[1][4] = p[2][1] - p[1][1];
   result[1][5] = P[2][1];
   result[1][6] = P[2][1];

/* Compute eigen problem */

   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] );

/* column[1] : internal strain energy for element 1 */
/* column[2] : elastic strain energy for element 1  */
/* column[3] : internal strain energy for element 2 */
/* column[4] : elastic strain energy for element 2  */

   element_energy = ColumnUnits( Matrix([total_step+1,4]), [Jou] );
   element_energy[1][1] = 0 Jou;
   element_energy[1][2] = 0 Jou;
   element_energy[1][3] = 0 Jou;
   element_energy[1][4] = 0 Jou;

/* increase external load, structure determination */

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

   print "*** start step no = ", k, "\n";

   time = time + dt;

   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;

   /* Compute effective incremental loading */

   dPeff = d_P + ((4/dt)*BigM + 2*BigC)*vel + 2*BigM*accel;

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

   /* Compute effective stiffness from tangent stiffness */

      Keff = BigK + (2/dt)*BigC + (4/dt/dt)*BigM;

   /* Solve for d_displacement, d_velocity */

      d_p = Solve( Keff, dPeff );
      d_v = (2/dt)*d_p - 2*vel;

   /* Compute displacement, velocity */

      new_p   = p + d_p;
      new_vel = vel + d_v;

   /* State determination for each element -- check material yielding */
   /* and compute new stiffness                                       */

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

         if( ele==1 ) { d_pe = [ 0 m; d_p[1][1] ]; }
         if( ele==2 ) { d_pe = [ d_p[1][1]; d_p[2][1] ]; }
         Q  = Q_saved[ele][1];   /* retrieve values from (j-1) */
         q  = q_saved[ele][1];
         Dx = bx*Q;
         dx = bx*q;
         K  = tangent[ele][1];
         kx = tangent[ele][1];
         fx = 1/kx;

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

         d_q = QuanCast(L*d_pe);  /* element deformation increment */
         q   = q + d_q;

         /* iterate for convervence of element displacements */

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

            /* get new section tangent flexibility f(x) from new d(x)     */
            /* and section resisting force DR(x) from material properties */

            /* set flags and stresses and strains for each i-th iteration */

            if( index[ele][1] == 1 ) {
               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]; }
               index[ele][1] = index[ele][1] + 1;
            }

            yielding[ele][1]  = flag1[ele][1];
            pre_range[ele][1] = flag2[ele][1];
            pre_load[ele][1]  = flag3[ele][1];

            sr[ele][1] = sr_saved[ele][1];
            er[ele][1] = er_saved[ele][1];
            s0[ele][1] = s0_saved[ele][1];
            e0[ele][1] = e0_saved[ele][1];

            /* material is still in elastic range, does not have any plastic residual */

            if( yielding[ele][1] == 0 ) then {
            if( abs(dx) <= ey[ele][1] ) then {
                kx = Ks[ele][1];
                sx = 0 N;
                ex = 0 cm;
                yielding[ele][1] = 0;
                pre_range[ele][1] = 0;
            } else {
                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 {   /* plastic residual occurs , yielding[ele][1] = 1 */

            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 {
                /* line 2 */
                if( loading[ele][1]*dx <= loading[ele][1]*er[ele][1] ) then {
                    /* line 2 */
                    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 {
                    /* line 4 */
                    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 */
           }
           }
           }

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

        }  /* j, while loop to check element convergence */

        /* 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]);

        if( abs(Q) > 2*Fy[ele][1] ) then {
        if( kx == Ks[ele][1] ) then {
            delta_Q = abs(sr[ele][1]) - 2*Fy[ele][1];
        } else {
            delta_Q = abs(Q) - 2*Fy[ele][1];
        }
        delta_x = delta_Q*(1/Kt[ele][1]-1/Ks[ele][1]);
        element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1] + 0.5*delta_x*delta_Q;
        } else {
        if( abs(sr[ele][1]) > 2*Fy[ele][1] ) then {
        if( kx == Ks[ele][1] ) then {
            delta_Q = abs(sr[ele][1]) - 2*Fy[ele][1];
            delta_x = delta_Q*(1/Kt[ele][1]-1/Ks[ele][1]);
            element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1] + 0.5*delta_x*delta_Q;
        } else {
            if(((Q>0N)&&(sr[ele][1]>0N))||((Q<0N)&&(sr[ele][1]<0N))) then {
                element_energy[k+1][2*ele] = 0.5*Q*Q/Kt[ele][1];
            } else {
                element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1];
            }
        }
        } else {
            element_energy[k+1][2*ele] = 0.5*Q*Q/Ks[ele][1];
        }
        }

        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] */

     }  /* ele */

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

     /* assemble new internal load Fs, damping load Fd, and acceleration */

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

     Fd = BigC*new_vel;
     new_accel = M_inv*( P-Fs-Fd );

     /* Store analysis results */

     result[k+1][1] = P[2][1];
     result[k+1][2] = p[2][1];
     result[k+1][3] = p[1][1];
     result[k+1][4] = p[2][1] - p[1][1];
     result[k+1][5] = Fs[1][1] + Fs[2][1];
     result[k+1][6] = Fs[2][1];

     /* Compute eigen problem */

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

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

     /* Updating history: save flags and reversal points for each load step k */

     for( ele=1 ; ele<=2 ; ele=ele+1 ) {
          pre_load[ele][1] = loading[ele][1];
          flag1[ele][1] = yielding[ele][1];
          flag2[ele][1] = pre_range[ele][1];
          flag3[ele][1] = pre_load[ele][1];

          sr_saved[ele][1] = sr[ele][1];
          er_saved[ele][1] = er[ele][1];
          s0_saved[ele][1] = s0[ele][1];
          e0_saved[ele][1] = e0[ele][1];
          sx_saved[ele][1] = Q_saved[ele][1];
          ex_saved[ele][1] = q_saved[ele][1];
     }

     /* Update information for this step */

     P_old  = P;
     Fs_old = Fs;
     Fd_old = Fd;
     p      = new_p;
     vel    = new_vel;
     accel  = new_accel;
}  /* k in for() loop, increase to next time step */

PrintMatrix(result);
PrintMatrix(element_energy);

quit;