input-psps-cantilever

有限元程序

/*
 *  =====================================================================
 *  Plane stress analysis of rectangular cantilever beam subject to an
 *  end moment.
 *
 *  Written By : Clara Popescu                               Spring, 1997
 *  =====================================================================
 */ 

print "*** DEFINE PROBLEM SPECIFIC PARAMETERS \n\n";

NDimension         = 2;
NDofPerNode        = 2;
MaxNodesPerElement = 4;

StartMesh();

   node = 1;
   AddNode(node, [  0.0 cm ,  0.0 cm ]);
   node = 2;
   AddNode(node, [ 10.0 cm ,  0.0 cm ]);
   node = 3;
   AddNode(node, [ 10.0 cm , 10.0 cm ]);
   node = 4;
   AddNode(node, [  0.0 cm , 10.0 cm ]);

/* Attach elements to grid of nodes */

   elmtno = 1;
   AddElmt( elmtno, [ 1, 2, 3, 4 ], "name_of_elmt_attr");
	
/* Define element section and material properties */

   ElementAttr("name_of_elmt_attr") { type     = "PLANE_STRESS";
                                      section  = "mysection";
                                      material = "mymaterial"; }

   SectionAttr("mysection") { depth = 10.0cm;
                              width = 10.0cm; }

   MaterialAttr("mymaterial") { density = 10000 kg/m^3;  
                                poisson = 1.0/3.0;   
                                yield   = 36000;  
                                E       = 20000000 kN/m^2; }

/* Setup boundary conditions */ 

   FixNode( 1, [1,1] );
   FixNode( 4, [1,0] );

/* Add point nodal loads to end of cantilever */

   Fx = 10.0 kN; Fy = 0.0 kN;

   NodeLoad( 2, [  Fx,  Fy]);
   NodeLoad( 3, [  Fx,  Fy]);

/* Compile and print finite element mesh */

   EndMesh();
   PrintMesh();

/* Compute Mass and Stiffness Matrices */

   stiff = Stiff();
   eload = ExternalLoad();

/* Solve static analysis problem */

   print "\n*** STATIC ANALYSIS PROBLEM \n\n";

   displ  = Solve(stiff, eload);
   iload  = InternalLoad(displ);

   PrintMatrix(iload);
   PrintDispl(displ);

/* Systematically retrieve stresses from individual elements */

   for( ii = 1; ii <= 1 ; ii = ii + 1 ) {
        print "Element : ", ii, "\n";
        actions = GetStress( [ii], displ );
        PrintMatrix( actions );
   }

/* Setup matrix for exptrapolation of stresses */

   M = Zero([4,4]);

   M[1][1] = (sqrt(3) + 1)^2;
   M[1][2] = 2;
   M[1][3] = (sqrt(3) - 1)^2;
   M[1][4] = 2;
   M[2][1] = 2;
   M[2][2] = (sqrt(3) + 1)^2;
   M[2][3] = 2;
   M[2][4] = (sqrt(3) - 1)^2;
   M[3][1] = (sqrt(3) - 1)^2;
   M[3][2] = 2;
   M[3][3] = (sqrt(3) + 1)^2;
   M[3][4] = 2;
   M[4][1] = 2;
   M[4][2] = (sqrt(3) - 1)^2;
   M[4][3] = 2;
   M[4][4] = (sqrt(3) + 1)^2;

   PrintMatrix(M);
   T = Inverse((1/12)*M);
   PrintMatrix(T*[ 1.5 Pa;
                   1.5 Pa;
                   2.0 Pa;
                   2.0 Pa ] );

   quit;
  
   PrintMatrix(T*[ actions[1][3];
                   actions[2][3];
                   actions[3][3];
                   actions[4][3] ]);

   PrintMatrix(T*[ actions[1][4];
                   actions[2][4];
                   actions[3][4];
                   actions[4][4] ]);

   quit;