Solid and Structural Mechanics Demo 4: A planar beam
This code solves an elastic beam problem in structural mechanics theory. The libray OFELI
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
IPF data("beam - 1.0",argv);
int verbose = data.getVerbose();
int output_flag = data.getOutput();
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:
Vect<double> b(ms), bc(ms), bf(ms);
- The global matrix and right-hand side are, as usual, filled
by looping over elements. For each element, an instance of class
is considered. This class deals with quadrilateral elements. The member function
Stiffness adds the contribution of the element stiffness
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
b = bf;
- We can now prescribe boundary conditions by penalty method
and solve the resulting linear system by factorization and backsubstitution.
- 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.
cout << "Element Axial Force Bending Moment Twisting Moment Shear Force" << endl;
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
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.
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" ?>
<data_file value="beam.dat" />
<verbose value="1" />
<Mesh dim="3" nb_dof="6">
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
<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
<PointForce x="1" dof="2">-0.01</PointForce>
<PointForce x="1" dof="3"> 0.01</PointForce>
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.