Input file for Cantilever with Material Nonlinearity

You will need Aladdin 2.0 to run this problem.

/*
 *  ============================================================================
 *  Load-displacement analysis of cantilever beam subject to incrementally
 *  increasing tip load
 *
 *     -- The cantilever is modeled with FIBER_2D elements having a bi-linear
 *        stiffness. Et = 0.1*E.
 *
 *  Written By : Wane-Jang Lin                                       August 1997
 *  ============================================================================
 */

/* Define problem specific parameters */ 

   NDimension         = 2;
   NDofPerNode        = 3;
   MaxNodesPerElement = 2;
   GaussIntegPts      = 5;

   SetUnitsType("US");
   StartMesh();

   total_node = 11;
   total_elmt = 10;

/* Generate grid of nodes for f.e. model */

   x = 0 in; y = 0 in;
   for( node=1 ; x<=50in ; node=node+1 ) {
      AddNode(node, [x, y]);
      x = x + 5 in;
   }

/* Attach elements to grid of nodes */

   for(elmt_no=1 ; elmt_no<=total_elmt ; elmt_no=elmt_no+1) {
      AddElmt(elmt_no,[elmt_no,elmt_no+1],"name_of_elmt_attr");
   }

/* Define element, section and material properties */

   ElementAttr("name_of_elmt_attr") { type     = "FIBER_2D";
                                      section  = "mysection";
                                      material = "mymaterial";
                                      fiber    = "fibername"; }
   h = 4.0 in;  b = 1.0 in;
   ks = 20000 psi;  kt = 0.1*ks;  fy = 500 psi;

   SectionAttr("mysection") { area  = b*h; 
                              width = b;
                              depth = h;
                              shear_factor  = 1.2; }
   MaterialAttr("mymaterial") { density = 2 lbf*sec^2/in/in^3;
                                poisson = 0.25;
                                E = ks;   Et = kt;  yield = fy;
                                shear_yield = 10000*fy; }
   no_fiber = 40;
   no_fiber_type = 1;
   dh = h/no_fiber;
   fcoord = Matrix([ 1, no_fiber ]);
   farea  = Matrix([ 1, no_fiber ]);
   fmap   = Matrix([ 1, no_fiber ]);
   for( i=1 ; i <= no_fiber ; i=i+1 ) {
      fcoord[1][i] = h/2 - dh/2 - (i-1)*dh;
      farea[1][i]  = b*h/no_fiber;
      fmap[1][i]   = 1;
   }
   fattr  = Matrix([ 3, no_fiber_type ]);
   for( j=1 ; j <= no_fiber_type ; j=j+1 ) {
      fattr[1][j] = ks;
      fattr[2][j] = kt;
      fattr[3][j] = fy;
   }

   FiberAttr(no_fiber,"fibername"){FiberMaterialAttr = fattr;
                                   FiberCoordinate   = fcoord;
                                   FiberArea         = farea;
                                   FiberMaterialMap  = fmap; }

/* Setup Boundary Conditions and End Mesh */

   FixNode(1, [1,1,1]);
   EndMesh();

/* Get the initial force and displacement. */

   Ks = Stiff();
   NodeLoad( total_node, [ 0 lbf, 0 lbf, 0 lbf*in ] ); 
   P_old = ExternalLoad();
   displ = Solve( Ks, P_old );

/* Setup the response matrices */

   total_step   = 60;
   tip_response = ColumnUnits( Zero([total_step+1,3]), [lbf,in,rad] );
   displacement = ColumnUnits( Zero([ total_node,6]), [in]);
   rotation     = ColumnUnits( Zero([ total_node,6]), [rad]);
   curvature    = ColumnUnits( Zero([ total_elmt,6]), [rad/in]);

/* Store the initial tip response */

   load_dir = 2;
   dof     = GetDof([total_node]);
   max_dof = dof[1][load_dir];
   tip_response[1][1] = P_old[max_dof][1];
   tip_response[1][2] = displ[max_dof][1];
   tip_response[1][3] = displ[max_dof+1][1];

/* Add nodal load at each step */

   flag = -10; index = 1;
   for( step=1 ; step <= total_step ; step=step+1 ) {

      dPk = [ 0 lbf, 1 lbf, 0 lbf*in ];
      NodeLoad( total_node, dPk ); 
      P_new = ExternalLoad();
      dP    = P_new - P_old;
      P_old = P_new;

      /* Newton-Raphson Iteration */

      while( L2Norm(dP) > 0.001 ) {

         dp    = Solve( Ks, dP );
         displ = displ + dp;

         ElmtStateDet( dp );

         Ks = Stiff();               /* Compute new stiffness    */
         PR = InternalLoad( displ ); /* Compute internal load    */
         dP = P_new - PR;            /* Compute unbalanced force */

      }  /* end of while loop, dP converges */

      /* Update element history data and save tip response */

      UpdateResponse();
      tip_response[step+1][1] = P_new[max_dof][1];
      tip_response[step+1][2] = displ[max_dof][1];
      tip_response[step+1][3] = displ[max_dof+1][1];

      /* store all element results every 10 steps */

      flag = flag+1;
      if( flag == 0 ) { 
         for( node=2 ; node<=total_node ; node=node+1 ) {
             node_displacement         = GetDispl([node],displ);
             displacement[node][index] = node_displacement[1][2];
             rotation[node][index]     = node_displacement[1][3];
         }
         for(elmt_no=1 ; elmt_no<=total_elmt ; elmt_no=elmt_no+1){
             curvature[elmt_no][index] = (rotation[elmt_no+1][index]
                                        - rotation[elmt_no][index])/(5 in);
         }
         flag  = -10;
         index = index+1;
      }
   } /* end of load step */

/* Print out the result matrices */

   PrintMatrix(tip_response);
   PrintMatrix(displacement,rotation);
   PrintMatrix(curvature);

   print "\n";
   print "===========================\n";
   print "Nonlinear Analysis Complete\n";
   print "===========================\n";

   quit;