Cantilever beam using beam elements

Previously, a thick cantilever beam was modeled with volume elements. In the present section quadratic beam elements are used for a similar exercise (Section 6.2.28). Beam elements are easy to define: they consist of three nodes on a line. Internally, they are expanded into volumetric elements. There are two types of beam elements: B32 elements, which are expanded into C3D20 elements, and B32R (reduced integration) elements, which are expanded into C3D20R elements. Based on the results in the present section, the B32R element is highly recommended. The B32 element, on the other hand, should be avoided especially if section forces are needed.

The first cantilever beam which is looked at is 100 mm long and has a square cross section of 2 x 2 mm$^2$. The axis of the beam is along the global z-direction. This beam is modeled with just one element and loaded at its end by a unit force in x-direction, Figure 44. We are interested in the stresses at integration point a and at node b, the section forces at the beam's fixed end, and the displacement in x at the free end. The location of the integration point a is at $x=-1/\sqrt{3}$, $y=1/\sqrt{3}$ and $z=50(1+1/\sqrt{3})$, the nodal coordinates of b are $x=-1$, $y=1$ and $z=100$ [17]. The material is isotropic linear elastic with a Young's modulus of 100,000 MPa and a Poisson's ratio of 0.3.

Figure 44: Geometry of the beam

The input deck for this example is very similar to the simplebeam.inp example in the test suite (Figure 45).

Figure 45: Input deck for the beam

The stresses at the integration points are obtained by a *EL PRINT card, the stresses at the nodes by the OUTPUT=3D option (default) on the *EL FILE card, whereas for the section forces the SECTION FORCES option on the same card is used (this option is mutually exclusive with the OUTPUT=3D option). The displacements are best obtained in the non-expanded view, i.e. using the OUTPUT=2D option. This means that for the present results the example had to be run twice: once with the OUTPUT=3D option and once with the SECTION FORCES option.

The results are summarized in Table 3. The $\{$mm$,$   N$,$   s$,$   K$\}$ system is used. The reference results are analytical results using simple beam theory [55]. The agreement is overwhelming. The stresses at the integration points match exactly, so do the extrapolated normal stresses to the nodes. The shear stresses need special attention. For a beam the shear stress varies parabolically across the section. A quadratic volumetric element can simulate only a linear stress variation across the section. Therefore, the parabolic variation is approximated by a constant shear stress across the section. Since the reduced integration points (at $\pm 1/\sqrt{3}$) happen to be points at which the parabolic stress variation attains its mean value the values at the integration points are exact! The extrapolated values to the nodes take the same constant value and are naturally wrong since the exact value at the corners is zero.

The section forces are obtained by

1. calculating the stresses at the integration points (inside the element, such as integration point a)
2. extrapolating those stresses to the corner nodes (such as node b)
3. calculating the stresses at the middle nodes by interpolation between the adjacent corner nodes
4. interpolating the stresses at all nodes within a section face onto the reduced integration points within the face (such as integration point c, using the shape functions of the face)
5. integrating these stresses numerically.

As shown by Table 3 this procedure yields the correct section forces for the square beam.

The displacements at the beam tip are off by 10 %. The deformation of a beam subject to a shear force at its end is third order, however, the C3D20R element can only simulate a quadratic behavior. The deviation is reduced to 2.4 % by using 5 elements (Table 4). Notice that integration point a is now closer to the fixation (same position is before but in the element adjacent to the fixation).

Table 3: Results for the square section beam subject to bending (1 element).

result	value	reference
$\sigma_{zz}(a)$	34.151	34.151
$\sigma_{xz}(a)$	-0.25	-0.25
$F_{xx}$	-1.	-1.
$M_{yy}$	100.	100.
$\sigma_{zz}(b)$	75.	75.
$\sigma_{xz}(b)$	-0.25	0.
$u_{x}$	2.25	2.50

Table 4: Results for the square section beam subject to bending (5 elements).

result	value	reference
$\sigma_{zz}(a)$	41.471	41.471
$\sigma_{xz}(a)$	-0.25	-0.25
$F_{xx}$	-1.	-1.
$M_{yy}$	100.	100.
$\sigma_{zz}(b)$	75.	75.
$\sigma_{xz}(b)$	-0.25	0.
$u_{x}$	2.44	2.50

The same beam was now subjected to a torque of 1 Nmm at its free end. The results are summarized in Table 5.

Table 5: Results for the square section beam subject to torsion (1 element).

result	value	reference
$\sigma_{xz}(a)$	-0.21651	-
$\sigma_{yz}(a)$	-0.21651	-
$M_{zz}$	1.	1.
$\sigma_{xz}(b)$	-0.375	0
$\sigma_{yz}(b)$	-0.375	0
$u_{y}$	$9.75 \cdot 10^{-4}$	$1.1525 \cdot 10^{-3}$

The torque is matched perfectly, the torsion at the end of the beam ($u_y$ is the displacement in y-direction at the corresponding node of node b) is off by 15 % [55]. The shear stresses at node b are definitely not correct (there is no shear stress at a corner node), however, the integration of the values interpolated from the nodes at the facial integration points yields the exact torque! Using more elements does not change the values in Table 5.

The same exercise is now repeated for a circular cross section (radius = 1 mm, same length, boundary conditions and material data as for the rectangular cross section). For such a cross section the vertex nodes of the element lie at $x,y = \pm 0.7071, \pm 0.7071$, whereas the middle nodes lie at $x,y = 0,\pm 1$ and $x,y = \pm 1,0$. The integration points are located at $x,y = \pm 0.5210$. The results for bending with just one element are shown in Table 6 and with 5 elements in Table 7.

Table 6: Results for the circular section beam subject to bending (1 element).

result	value	reference
$\sigma_{zz}(a)$	34.00	52.26
$\sigma_{xz}(a)$	-0.322	-0.318
$F_{xx}$	-0.99996	-1.
$M_{yy}$	58.7	100.
$\sigma_{zz}(b)$	62.8	90.03
$\sigma_{xz}(b)$	-0.322	-0.318
$u_{x}$	2.91	4.24

Table 7: Results for the circular section beam subject to bending (5 elements).

result	value	reference
$\sigma_{zz}(a)$	59.77	63.41
$\sigma_{xz}(a)$	-0.322	-0.318
$F_{xx}$	-0.99996	-1.
$M_{yy}$	102.	100.
$\sigma_{zz}(b)$	109.	90.03
$\sigma_{xz}(b)$	-0.322	-0.318
$u_{x}$	3.86	4.24

For just one element the shear stress is quite close to the analytical value, leading to a even better match of the shear force. This is remarkable an can only be explained by the fact that the cross area of the piecewise quadratic approximation of the circular circumference is smaller and exactly compensates the slightly higher shear stress. A similar effect will be noticed for the torque. The normal stress, however, is far off at the integration points as well as at the nodes leading to a bending moment which is way too small. The same applies to the deformation in x-direction. Using five elements leads to a significant improvement: the bending moment is only 2 % off, the deformation at the free end 9 %. Here again one can argue that the deformation is of cubic order, whereas a quadratic element can only simulate a quadratic change. Using more elements consequently improves the results.

The results for a torque applied to a circular cross section beam is shown in Table 8 (1 element; the results for 5 elements are identical).

Table 8: Results for the circular section beam subject to torsion (1 element).

result	value	reference
$\sigma_{xz}(a)$	-0.309	-0.331
$\sigma_{yz}(a)$	-0.309	-0.331
$M_{zz}$	0.999994	1.
$\sigma_{xz}(b)$	-0.535	-0.450
$\sigma_{yz}(b)$	-0.535	-0.450
$u_{y}$	$1.54 \cdot 10^{-3}$	$1.66 \cdot 10^{-3}$

Again, it is remarkable that the torque is perfectly matched, although the shear stress at the integration points is 6 % off. This leads to shear values at the vertex nodes which are 19 % off. Interpolation to the facial integration points yields shear stresses of -0.305 MPa. Integration of these stresses finally leads to the perfect torque values. The torsion angle at the end of the beam is 7 %off.

Summarizing, one can state that the use of C3D20R elements leads to quite remarkable results:

• For a rectangular cross section:
  • the section forces are correct
  • the stresses at the integration points are correct
  • the displacements for bending are correct, provided enough elements are used
  • the torsion angle is somewhat off (15 %).
• For a circular cross section:
  • the shear force and torque section forces are correct
  • the bending moment is correct if enough elements are used
  • the displacements for bending are correct, provided enough elements are used
  • the torsion angle is somewhat off (7 %).

It is generally recommended to calculate the stresses from the section forces. The only drawback is the C3D20R element may lead to hourglassing, leading to weird displacements. However, the mean of the displacements across the cross section is usually fine. An additional problem which can arise is that nonlinear geometric calculations may not converge due to this hourglassing. This is remedied in CalculiX by slightly perturbing the coordinates of the expanded nodes (by about 0.1 %).

A similar exercise was performed for the B32 element, however, the results were quite discouraging. The section forces were, especially for bending, way off.