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;