Demo: A planar beam

This code solves an elastic beam problem in structural mechanics theory. The libray contains the class Beam3DL2 that handles elastic beams in the space with 6 degrees of freedom for each node. These degrees of freedom are supported by the vertices of each element line.

  • We start by defining an instance of class IPF to deal with project data. We extract useful information from this instance: A verbosity and an output flag to deal with output levels. We next construct the Mesh and the Prescription instances. 

       IPF data("beam - 1.0",argv[1]);
   int verbose = data.getVerbose();
   int output_flag = data.getOutput();
   Mesh ms(data.getMeshFile());
   Prescription p(ms,data.getDataFile());

  • The linear system uses a symmetric matrix. For this, we we use an instance of class SkSMatrix<double>. Then we declare a vector that will contain the right-hand side and then the solution. Next, we read boundary conditions and prescribed concentrated loads from the prescription file. These are stored in vectors bc and bf respectively: 

       SkSMatrix<double> A(ms);
   Vect<double> b(ms), bc(ms), bf(ms);
   p.get(BOUNDARY_CONDITION,bc,0);
   p.get(POINT_FORCE,bf);

  • The global matrix and right-hand side are, as usual, filled by looping over elements. For each element, an instance of class Beam3DL2 is considered. This class deals with quadrilateral elements. The member function Stiffness adds the contribution of the element stiffness matrix. The other member functions have already been considered in previous lessons. Note that the load vector is a copy of the concentrated load vector. No need for assembly here. 

       MeshElements(ms) {
      Beam3DL2 eq(theElement,0.1,0.1,0.1);
      eq.Stiffness();
      eq.ElementAssembly(A);
      eq.ElementAssembly(b);
   }
   b = bf;

  • We can now prescribe boundary conditions by penalty method and solve the resulting linear system by factorization and backsubstitution. 

       A.Prescribe(b,bc);
   A.solve(b);

  • Let us now extract some additional information from the class. We want to output, for each element, the axial force, the bending moment, the twisting moment and the shear force. This is done simply by invoking corresponding member functions in the class Beam3DL2. Note that we need, before that, to construct the vector of node displacements. A dedicated constructor is made for that. Then we recall the member functions we need for post-processing. 

       Vect<double> u(ms,3);
   Beam3DL2 eq(ms,uf,u);
   cout << "Element     Axial Force       Bending Moment    Twisting Moment   Shear Force" << endl;
   MeshElements(ms) {
      Beam3DL2 eq(theElement,0.1,0.1,0.1,b);
      cout << setw(8) << theElementLabel;
      cout << " " << eq.AxialForce();
      cout << " " << eq.BendingMoment();
      cout << " " << eq.TwistingMoment();
      cout << " " << eq.ShearForce() << endl;
   }

  • We finally show here how to create a deformed mesh for a graphical presentation. We have for this the function DeformMesh that modifies node coordinates by adding the displacement vector to the coordinates. The function can also multiply the displacement by an amplification factor which is equal here to 1 by default. Then we store the resulting mesh into a new file. 

       DeformMesh(ms,u);
   ms.put("deformed_beam.m");

  

An example

Let us run this program with an example of a straight bar [0,1]. We define the project file: 

<?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
   <info>
      <title />
      <date />
      <author />
   </info>
                              
   <Project name="beam">
      <mesh_file value="beam.dat"/>
      <data_file value="beam.dat" />
      <verbose value="1" />
      <output value="0">
      <save value="1"/>
      <plot value="1"/>
      <init value="1"/>
      <plot_file value="beam.d"/>
   </Project>
                              
   <Mesh dim="3" nb_dof="6">
      <Nodes>
         0.00000000e+00      0.00000000e+00      0.00000000e+00   111111
         1.25000000e-01      0.00000000e+00      0.00000000e+00        0
         2.50000000e-01      0.00000000e+00      0.00000000e+00        0
         3.75000000e-01      0.00000000e+00      0.00000000e+00        0
         5.00000000e-01      0.00000000e+00      0.00000000e+00        0
         6.25000000e-01      0.00000000e+00      0.00000000e+00        0
         7.50000000e-01      0.00000000e+00      0.00000000e+00        0
         8.75000000e-01      0.00000000e+00      0.00000000e+00        0
         1.00000000e+00      0.00000000e+00      0.00000000e+00        0
      </Nodes>
      <Elements shape="line"  nodes="2">
         1       2       1          2       3       1          
         3       4       1          4       5       1          
         5       6       1          6       7       1
         7       8       1          8       9       1   
      </Elements>
   </Mesh>
   <Prescription>
      <PointForce x="1" dof="2">-0.01</PointForce>
      <PointForce x="1" dof="3"> 0.01</PointForce>
   </Prescription>
</OFELI_File>

Let us give some explanations about this project file:

  • The project file contains all data necessary for the code: General problem parameters, mesh data and prescriptions.
  • The finite element mesh consists in 8 beam elements, where the first node (x=0) has all its degrees of freedom constrained.
  • We impose a load concentrated at the end point (x=1): The degrees of freedom 2 and 3 are loaded. Note that since we don't prescribe any boundary condition, then by default, the first end (x=0) is fixed.