input-springs2

有限元程序

/* 
 *  =========================================================================
 *  Force-displacement calculation for mass-spring system having springs 
 *  with bi-linear force-displacement characteristics.
 *
 *  We will compute : force displacement curve.
 *                    system-level energy versus time.
 *                    element-level energy versus time.
 *
 *  Written By: Wane-Jang Lin and Mark Austin       October 1997 - March 1998
 *  =========================================================================
 */ 

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

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

   total_step = 400;

/* Allocate storage matrix for force-displacement 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                */

   result = ColumnUnits( Matrix([total_step+1,4]), [N,cm,cm,cm] );
   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];

/* Allocate storage matrix system/element energy calculations   */

/* column[1] : total external work for whole system             */
/* column[2] : total internal energy = [3]+[4] for whole system */
/* column[3] : elastic strain energy for whole system           */
/* column[4] : plastic strain energy for whole system           */

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

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

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

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

   /* Define incremental loading */

   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;

   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;

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

            /* 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, and therefore 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 {
                  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 */
              }
              }
            }

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

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

         dElastic = 0.5*Q*Q/Ks[ele][1];
         element_energy[k+1][2*ele]   = dElastic;

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

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

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

   /* reassemble the system energy */

   dWork = 0.5*(result[k+1][1]+result[k][1])*(result[k+1][2]-result[k][2]);

   system_energy[k+1][1] =  system_energy[k][1]   + dWork;
   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];

} 

PrintMatrix(result);
PrintMatrix(element_energy);
PrintMatrix(system_energy);

quit;