# Beams - structural beams

## Contents

# Beams - structural beams#

This notebook shows how to use Beam elements model cantilever beams, how to get the results and how those results compare to theory.

A “Beam” in DAVE consists of a number of discrete segments and has a different formulation than the beams typically found in FEM packages.

```
from DAVE import *
s = Scene();
```

```
Equilibrium-core version = 2.48 from c:\python\miniconda3\envs\book\lib\site-packages\pyo3d.cp39-win_amd64.pyd
default resource folders:
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\resources
C:\Users\beneden\DAVE_models
C:\data\Dave\Book\DAVE-book\DAVE-notebooks
Blender found at: C:\Program Files\Blender Foundation\Blender 3.3\blender.exe
```

Beams are elastic nodes that can be created between axis systems.

```
left = s.new_axis("left")
right = s.new_axis("right")
beam = s.new_beam("beam", nodeA=left, nodeB = right,
EA = 1000,
EIy = 1e5, EIz = 1e5, GIp = 1e3,
L=10,
mass = 0,
n_segments=20)
right.set_free()
s.solve_statics()
```

```
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\scene.py:1668: UserWarning: new_axis is deprecated, use new_frame instead
warnings.warn("new_axis is deprecated, use new_frame instead")
```

```
True
```

```
from DAVE.jupyter import *;
```

```
show(s, camera_pos=(5,-10,6), lookat = (5,0,2), width=1000, height = 200)
```

```
Warning! Please use "settings.use_parallel_projection" instead!
```

```
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\jupyter\jupyter.py:194: UserWarning: VTK/Vedo issue: plotter is None
warnings.warn("VTK/Vedo issue: plotter is None")
```

A Beam node is rigidly fixed to the axis systems at its ends. It is fixed in the positive X-direction, so it departs from node A along the X-axis and arrives at ndoe B from the negative X-direction.

At the momement we have a massless beam with its left side fixed (axis system “left” is fixed) and its right size free (the right end of the beam is fixed to axis system “right”, but “right” is not connected to anything and is free to move).

This is one of those cases that frequently occur in textbooks and for which we can easily check the results.

First place a point on the right side so that we can apply a force:

```
s.new_point("point", parent = right)
force = s.new_force("force", parent = "point")
```

## Cantilever beam with a force#

```
force.force = (0,0,-100)
```

```
s.solve_statics()
```

```
True
```

```
show(s, camera_pos=(5,-10,0), lookat = (5,0,0), width=1000, height = 300)
```

```
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\visual_helpers\vtkHelpers.py:533: RuntimeWarning: invalid value encountered in divide
axis = axis / length
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\visual_helpers\vtkHelpers.py:594: RuntimeWarning: invalid value encountered in divide
axis = axis / length
```

```
Warning! Please use "settings.use_parallel_projection" instead!
```

Theory tells us that the deflection and rotation of the right end should be:

\(\delta = -PL^3 / 3EI\)

and

\(\Theta = PL^2 / 2EI\)

Get the values from the model:

```
L = beam.L
EI = beam.EIy
P = -force.force[2]
```

And calculate their values

The deflection is the vertical displacement of the axis system on the right side of the beam:

```
right.z
```

```
-0.33604653085104813
```

```
-P*L**3 / (3*EI)
```

```
-0.3333333333333333
```

Theta is the slope. The slope can be calculated from the rotation of the axis system, but it is also available as “tilt”, which is a percentage.

```
right.tilt_y / 100
```

```
0.05014161613258149
```

```
P*L**2 / (2*EI)
```

```
0.05
```

It is also possible to obtain the forces in the beam.

The moments are available at the nodes.

```
import matplotlib.pyplot as plt
```

```
plt.plot(beam.X_nodes, beam.bending[:,1])
plt.xlabel('Distance along the beam [m]')
plt.ylabel('Moment about y-axis [kN*m]');
```

## Cantilever with moment#

We can do the same with a moment

```
force.force = (0,0,0)
force.moment = (0,1000,0)
```

```
s.solve_statics()
```

```
True
```

```
show(s, camera_pos=(5,-10,3), lookat = (5,0,-1), width=1000, height = 300)
```

```
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\visual_helpers\vtkHelpers.py:533: RuntimeWarning: invalid value encountered in divide
axis = axis / length
```

```
Warning! Please use "settings.use_parallel_projection" instead!
```

```
c:\python\miniconda3\envs\book\lib\site-packages\DAVE\jupyter\jupyter.py:194: UserWarning: VTK/Vedo issue: plotter is None
warnings.warn("VTK/Vedo issue: plotter is None")
```

Theory tells us that the deflection and rotation of the right end should be:

\(\delta = ML^2 / 2EI\)

and

\(\Theta = ML / EI\)

Get the values from the model:

```
M = force.moment[1] # moment about Y-axis
```

```
M*L**2 / (2*EI)
```

```
0.5
```

```
-right.z
```

```
0.4995839925972531
```

```
M*L / EI
```

```
0.1
```

```
right.tilt_y / 100
```

```
0.09983341664682813
```

## Tension and torsion#

Beam can also take tension and torsion.

```
force.force = (100,0,0)
force.moment = (100,0,0)
```

```
s.solve_statics()
```

```
True
```

```
show(s, camera_pos=(5,-10,3), lookat = (5,0,-1), width=1000, height = 300)
```

```
Warning! Please use "settings.use_parallel_projection" instead!
```

Tension and torsion forces in the beam are available not at the nodes but at the centers for the beam segments.

```
plt.plot(beam.X_midpoints, beam.torsion);
```

```
plt.plot(beam.X_midpoints, beam.tension);
```

**Now this looks horrible!**

Fortunately this is just the way matplotlib displays data which is almost constant.

Looking at the data directly it appears to be comfortingly constant:

```
beam.torsion
```

```
[99.99999999999781,
99.99999999999814,
99.99999999999709,
99.99999999999842,
99.99999999999537,
99.99999999999955,
99.9999999999967,
99.99999999999459,
100.00000000000048,
99.99999999999915,
99.99999999999892,
100.00000000000448,
100.00000000000048,
99.9999999999956,
100.00000000000381,
99.99999999998425,
100.00000000000159,
100.00000000000026,
100.00000000000136,
99.99999999998536]
```

```
import matplotlib
matplotlib.rcParams['axes.formatter.useoffset'] = False
```

```
plt.plot(beam.X_midpoints, beam.torsion);
plt.ylim((99.9, 100.1));
```

```
plt.plot(beam.X_midpoints, beam.tension);
plt.ylim((99.9, 100.1));
```

## Distributed load#

Modelling a distributed load is not possible in this way. Loads can only be added at points.

To model a distributed load we would have to manually discretise the model and add discrete loads at each node.

But these is a shortcut that we can take: The beams can have a mass. We can use the mass to model a distributed load as gravity.

```
s.delete('force')
```

```
q = 10 # kN/m
```

```
beam.mass = beam.L * q / 9.81
```

```
s.solve_statics()
```

```
True
```

```
show(s, camera_pos=(5,-10,0), lookat = (5,0,0), width=1000, height = 300)
```

```
Warning! Please use "settings.use_parallel_projection" instead!
```

Theory tells us that the deflection and rotation of the right end should be:

\(\delta = qL^4 / 8EI\)

and

\(\Theta = qL^3 / 6EI\)

```
-q*L**4 / (8*EI)
```

```
-0.125
```

```
right.z
```

```
-0.12539340701469195
```

```
q*L**3 / (6*EI)
```

```
0.016666666666666666
```

```
right.tilt_y / 100
```

```
0.01668937486888241
```

## Large deflections#

All the tests done so far were performed on small deflections. That is for a good reason. The Euler/Bernouilli beam theory is only applicable to small displacements. This is because in this theory the moment is derived from the deflection:

\(d^2/dx^2(EI d^2w / dx^2) = q\)

so the moment in the beam is proprtional to the change in slope. This is only valid for small changes.

The implementation in DAVE **is valid for large deflections**, it is just that the analytical formula used above are not.

## Number of segments#

Beams use a discrete model (see theory https://davedocs.online/beams.html).

When a higher number of segments is used:

The model become more accurate

The solver becomes slower

(The numerical errors build up)

So the number of nodes should be selected with some caution. In general 20 seems to be a good number for beams although in some cases (for example pure tension or torsion) a single segment is just as good.

It is easy to use the “plot_effect” function in scene to check the effect of the number of segments.

Here we calculate the effect of the number of segments on the vertical position of the second end of the beam:

```
s.plot_effect(evaluate="s['point'].gz",
change_property="n_segments",
change_node="beam",
start=1,
to=25.0,
steps=25);
```

```
setting 1.0 results in -0.25043539538010456
setting 2.0 results in -0.15639089359229383
setting 3.0 results in -0.13899103714473474
```

```
setting 4.0 results in -0.1329040421012639
setting 5.0 results in -0.1300872247358817
```

```
setting 6.0 results in -0.1285572544156294
setting 7.0 results in -0.12763482040356924
```

```
setting 8.0 results in -0.1270361175408637
```

```
setting 9.0 results in -0.12662568342623784
```

```
setting 10.0 results in -0.12633209478116098
```

```
setting 11.0 results in -0.12611496195310684
```

```
setting 12.0 results in -0.12594964178054247
```

```
setting 13.0 results in -0.12582112540100712
```

```
setting 14.0 results in -0.12571911178709774
```

```
setting 15.0 results in -0.12563692195711026
```

```
setting 16.0 results in -0.12556940415251053
```

```
setting 17.0 results in -0.12551358019856465
```

```
setting 18.0 results in -0.12546679931017593
```

```
setting 19.0 results in -0.12542720873744273
```

```
setting 20.0 results in -0.1253934070146907
```

```
setting 21.0 results in -0.12536431820383623
```

```
setting 22.0 results in -0.12533910520139316
```

```
setting 23.0 results in -0.12531710883253766
```

```
setting 24.0 results in -0.1252978043245916
```

```
setting 25.0 results in -0.12528076971477245
```

– end –